Block Krylov Subspace Research Articles

A review of block Krylov subspace methods for multisource electromagnetic modelling

Practical applications of controlled-source electromagnetic (EM) modelling require solutions for multiple sources at several frequencies, thus leading to a dramatic increase of the computational cost. In this paper, we present an approach using block Krylov subspace solvers that are iterative methods especially designed for problems with multiple right-hand sides (RHS). Their main advantage is the shared subspace for approximate solutions, hence, these methods are expected to converge in less iterations than the corresponding standard solver applied to each linear system. Block solvers also share the same preconditioner, which is constructed only once. Simultaneously computed block operations have better utilization of cache due to the less frequent access to the system matrix. In this paper, we implement two different block solvers for sparse matrices resulting from the finite-difference and the finite-element discretizations, discuss the computational cost of the algorithms and study their dependence on the number of RHS given at once. The effectiveness of the proposed methods is demonstrated on two EM survey scenarios, including a large marine model. As the results of the simulations show, when a powerful preconditioning is employed, block methods are faster than standard iterative techniques in terms of both iterations and time.

A block MINRES algorithm based on the band Lanczos method

We develop a block minimum residual (MINRES) algorithm for symmetric indefinite matrices. This version is built upon the band Lanczos method that generates one basis vector of the block Krylov subspace per iteration rather than a whole block as in the block Lanczos process. However, we modify the method such that the most expensive operations are still performed in a block fashion. The benefit of using the band Lanczos method is that one can detect breakdowns from scalar values arising in the computation, allowing for a handling of breakdown which is straightforward to implement. We derive a progressive formulation of the MINRES method based on the band Lanczos process and give some implementation details. Specifically, a simple reordering of the steps allows us to perform many of the operations at the block level in order to take advantage of communication efficiencies offered by the block Lanczos process. This is an important concern in the context of next-generation super computing applications. We also present a technique allowing us to maintain the block size by replacing dependent Lanczos vectors with pregenerated random vectors whose orthogonality against all Lanczos vectors is maintained. Numerical results illustrate the performance on some sample problems. We present experiments that show how the relationship between right-hand sides can effect the performance of the method.

Block Krylov subspace methods for large-scale matrix computations in control

In this paper we show how to improve the approximate solution of the large Sylvester equation obtained by an arbitrary method. Such problems appear in many areas of control theory such as the computation of Hankel singular values, model reduction algorithms and others. Moreover, we propose a new method based on refinement process and weighted block Arnoldi algorithm for solving large Sylvester matrix equation. The numerical tests report the effectiveness of these methods.

Iterative across-time solution of linear differential equations: Krylov subspace versus waveform relaxation

The aim of this paper is two-fold. First, we propose an efficient implementation of the continuous time waveform relaxation (WR) method based on block Krylov subspaces. Second, we compare this new WR–Krylov implementation against Krylov subspace methods combined with the shift and invert (SAI) technique. Some analysis and numerical experiments are presented. Since the WR–Krylov and SAI–Krylov methods build up the solution simultaneously for the whole time interval and there is no time stepping involved, both methods can be seen as iterative across-time methods. The key difference between these methods and standard time integration methods is that their accuracy is not directly related to the time step size.

A deflated conjugate gradient method for multiple right hand sides and multiple shifts

We consider the task of computing solutions of linear systems that only differ by a shift with the identity matrix as well as linear systems with several different right-hand sides. In the past, Krylov subspace methods have been developed which exploit either the need for solutions to multiple right-hand sides (e.g. deflation type methods and block methods) or multiple shifts (e.g. shifted CG) with some success. In this paper we present a block Krylov subspace method which, based on a block Lanczos process, exploits both features--shifts and multiple right-hand sides--at once. Such situations arise, for example, in lattice quantum chromodynamics (QCD) simulations within the Rational Hybrid Monte Carlo (RHMC) algorithm. We present numerical evidence that our method is superior compared to applying other iterative methods to each of the systems individually as well as, in typical situations, to shifted or block Krylov subspace methods.

A block Krylov subspace time‐exact solution method for linear ordinary differential equation systems

SUMMARYWe propose a time‐exact Krylov‐subspace‐based method for solving linear ordinary differential equation systems of the form y ′ = − Ay + g(t) and y ′ ′ = − Ay + g(t), where y(t) is the unknown function. The method consists of two stages. The first stage is an accurate piecewise polynomial approximation of the source term g(t), constructed with the help of the truncated singular value decomposition. The second stage is a special residual‐based block Krylov subspace method. The accuracy of the method is only restricted by the accuracy of the piecewise polynomial approximation and by the error of the block Krylov process. Because both errors can, in principle, be made arbitrarily small, this yields, at some costs, a time‐exact method. Numerical experiments are presented to demonstrate efficiency of the proposed method. Copyright © 2013 John Wiley & Sons, Ltd.

Approximation of the matrix exponential operator by a structure-preserving block Arnoldi-type method

The approximation of exp(A)V, where A is a real matrix and V a rectangular matrix, is the key ingredient of many exponential integrators for solving systems of ordinary differential equations. The use of Krylov subspace techniques in this context has been actively investigated; see Calledoni and Moret (1997) [10], Hochbruck and Lubich (1997) [17], Saad (1992) [20]. An appropriate structure preserving block method for approximating exp(A)V, where A is a large square real matrix and V a rectangular matrix, is given in Lopez and Simoncini (2006) [18]. A symplectic Krylov method to approximate exp(A)V was also proposed in Agoujil et al. (2012) [2] with V∈R2n×2. The purpose of this work is to describe a structure preserving block Krylov method for approximating exp(A)V when A is a Hamiltonian or skew-Hamiltonian 2n-by-2n real matrix and V is a 2n-by-2s matrix (s≪n). Our approach is based on block Krylov subspace methods that preserve Hamiltonian and skew-Hamiltonian structures.

Block Arnoldi-based methods for large scale discrete-time algebraic Riccati equations

In the present paper, we present block Arnoldi-based methods for the computation of low rank approximate solutions of large discrete-time algebraic Riccati equations (DARE). The proposed methods are projection methods onto block or extended block Krylov subspaces. We give new upper bounds for the norm of the error obtained by applying these block Arnoldi-based processes. We also introduce the Newton method combined with the block Arnoldi algorithm and present some numerical experiments with comparisons between these methods.

Block Krylov subspace methods for the computation of structural response to turbulent wind

In this paper the numerical computation of the dynamic response to turbulent wind excitations of slender structures is addressed. A numerical procedure capable to effectively estimate the three-dimensional structural behavior is proposed, based on a direct frequency domain approach. A probabilistic description of the wind velocity field, accounting for the correlation between the turbulence components, is combined to a linearized fluid–structure interaction model, under the quasi-steady hypothesis. We propose robust implementations of multiple right-hand side and multiple shift Krylov subspace methods with deflation of basis vectors, which allow us to efficiently analyze the dynamic response for a wide range of frequency values and wind time histories. Numerical experiments are reported with data stemming from real structure modeling.

Model-order reduction of kth order MIMO dynamical systems using block kth order Krylov subspaces

In this article, we study numerical methods for model-order reduction of large-scale kth order multi-input multi-output dynamical systems. We propose a new structure-preserving projection method, of which the projection subspace is a block kth order Krylov subspace based on a square matrix sequence and an initial rectangle matrix. A procedure, named as block kth order Arnoldi process, is presented for establishing an orthonormal basis of the projection subspace. Moreover, we show that the reduced system constructed by the new method has the same order of approximation as the standard block Krylov subspace method via linearization. Numerical experiments report the effectiveness of this method.

Stein implicit Runge–Kutta methods with high stage order for large-scale ordinary differential equations

The Runge–Kutta method is one of the most popular implicit methods for the solution of stiff ordinary differential equations. For large problems, the main drawback of such methods is the cost required at each integration step for computing the solution of a nonlinear system of equations. In this paper, we propose to reduce the cost of the computation by transforming the linear systems arising in the application of Newton's method to Stein matrix equations. We propose an iterative projection method onto block Krylov subspaces for solving numerically such Stein matrix equations. Numerical examples are given to illustrate the performance of our proposed method.

Deflated block Krylov subspace methods for large scale eigenvalue problems

We discuss a class of deflated block Krylov subspace methods for solving large scale matrix eigenvalue problems. The efficiency of an Arnoldi-type method is examined in computing partial or closely clustered eigenvalues of large matrices. As an improvement, we also propose a refined variant of the Arnoldi-type method. Comparisons show that the refined variant can further improve the Arnoldi-type method and both methods exhibit very regular convergence behavior.

Application of block Krylov subspace algorithms to the Wilson–Dirac equation with multiple right-hand sides in lattice QCD

It is well known that the block Krylov subspace solvers work efficiently for some cases of the solution of differential equations with multiple right-hand sides. In lattice QCD calculation of physical quantities on a given configuration demands us to solve the Dirac equation with multiple sources. We show that a new block Krylov subspace algorithm recently proposed by the authors reduces the computational cost significantly without losing numerical accuracy for the solution of the O ( a ) -improved Wilson–Dirac equation.

Order reduction of bilinear MIMO dynamical systems using new block Krylov subspaces

In this paper we study numerical methods for the model-order reduction of large-scale bilinear multi-input multi-output systems. A new projection method is proposed. The projection subspace is the union of some new block Krylov subspaces. We show that the reduced-order bilinear system constructed by the new method can match a desired number of moments of multivariable transfer functions corresponding to the kernels of Volterra series representation of the original system. Some numerical examples are presented to illustrate the effectiveness of the proposed method.

Block BiCGGR: a new Block Krylov subspace method for computing high accuracy solutions

Block BiCGGR: a new Block Krylov subspace method for computing high accuracy solutions

A Spectral Time-Domain Method for Computational Electrodynamics

Block Krylov subspace spectral (KSS) methods have previously been applied to the variable-coefficient heat equation and wave equation, and have demonstrated high-order accuracy, as well as stability characteristic of implicit time-stepping schemes, even though KSS methods are explicit. KSS methods for scalar equations compute each Fourier coefficient of the solution using techniques developed by Gene Golub and Gerard Meurant for approximating elements of functions of matrices by Gaussian quadrature in the spectral, rather than physical, domain. We show how they can be generalized to non-self-adjoint systems of coupled equations, such as Maxwell’s equations.

The block grade of a block Krylov space

The aim of the paper is to compile and compare basic theoretical facts on Krylov subspaces and block Krylov subspaces. Many Krylov (sub)space methods for solving a linear system Ax = b have the property that in exact computer arithmetic the true solution is found after ν iterations, where ν is the dimension of the largest Krylov subspace generated by A from r 0 , the residual of the initial approximation x 0 . This dimension is called the grade of r 0 with respect to A. Though the structure of block Krylov subspaces is more complicated than that of ordinary Krylov subspaces, we introduce here a block grade for which an analogous statement holds when block Krylov space methods are applied to linear systems with multiple, say s , right-hand sides. In this case, the s initial residuals are bundled into a matrix R 0 with s columns. The possibility of linear dependence among columns of the block Krylov matrix ( R 0 AR 0 ⋯ A ν - 1 R 0 ) , which in practical algorithms calls for the deletion (or, deflation) of some columns, requires extra care. Relations between grade and block grade are also established, as well as relations to the corresponding notions of a minimal polynomial and its companion matrix.

On Padé-type model order reduction of J-Hermitian linear dynamical systems

A simple, yet powerful approach to model order reduction of large-scale linear dynamical systems is to employ projection onto block Krylov subspaces. The transfer functions of the resulting reduced-order models of such projection methods can be characterized as Padé-type approximants of the transfer function of the original large-scale system. If the original system exhibits certain symmetries, then the reduced-order models are considerably more accurate than the theory for general systems predicts. In this paper, the framework of J-Hermitian linear dynamical systems is used to establish a general result about this higher accuracy. In particular, it is shown that in the case of J-Hermitian linear dynamical systems, the reduced-order transfer functions match twice as many Taylor coefficients of the original transfer function as in the general case. An application to the SPRIM algorithm for order reduction of general RCL electrical networks is discussed.

Restarted block Lanczos bidiagonalization methods

The problem of computing a few of the largest or smallest singular values and associated singular vectors of a large matrix arises in many applications. This paper describes restarted block Lanczos bidiagonalization methods based on augmentation of Ritz vectors or harmonic Ritz vectors by block Krylov subspaces.

Preserving geometric properties of the exponential matrix by block Krylov subspace methods

Given a large square real matrix A and a rectangular tall matrix Q, many application problems require the approximation of the operation \(\exp(A)Q\). Under certain hypotheses on A, the matrix \(\exp(A)Q\) preserves the orthogonality characteristics of Q; this property is particularly attractive when the associated application problem requires some geometric constraints to be satisfied. For small size problems numerical methods have been devised to approximate \(\exp(A)Q\) while maintaining the structure properties. On the other hand, no algorithm for large A has been derived with similar preservation properties. In this paper we show that an appropriate use of the block Lanczos method allows one to obtain a structure preserving approximation to \(\exp(A)Q\) when A is skew-symmetric or skew-symmetric and Hamiltonian. Moreover, for A Hamiltonian we derive a new variant of the block Lanczos method that again preserves the geometric properties of the exact scheme. Numerical results are reported to support our theoretical findings, with particular attention to the numerical solution of linear dynamical systems by means of structure preserving integrators.