Preconditioning the Restarted and Shifted Block FOM Algorithm for Matrix Exponential Computation

The approximation of $e^{tA}B$ where $A$ is a large sparse matrix and $B$ a rectangular matrix is the key ingredient in many scientific and engineering computations. A powerful tool to manage the matrix exponential function is to resort to a suitable rational approximation such as the Carath$\acute{\rm e}$odory-Fej$\acute{\rm e}$r approximation, whose core reduces to solve shifted linear systems with multiple right-hand sides. The restarted and shifted block FOM algorithm is a commonly used technique for this problem. However, determining good preconditioners for shifted systems that preserve the original structure is a difficult task. In this paper, we propose a new preconditioner for the restarted and shifted block FOM algorithm. The key is that the absolute values of the poles of the Carath$\acute{\rm e}$odory-Fej$\acute{\rm e}$r approximation are medium sized and can be much smaller than the norm of the matrix in question. The advantages of the proposed strategy are that we can precondition all the shifted linear systems simultaneously, and preserve the original structure of the shifted linear systems after restarting. Theoretical results are provided to show the rationality of our preconditioning strategy. Applications of the new approach to Toeplitz matrix exponential problem are also discussed. Numerical experiments illustrate the superiority of the new algorithm over many state-of-the-art algorithms for matrix exponential.

Related papers:

title: Computing the Action of the Matrix Exponential, with an Application to Exponential Integrators
abstract: A new algorithm is developed for computing $e^{tA}B$, where $A$ is an $n\times n$ matrix and $B$ is $n\times n_0$ with $n_0 \ll n$. The algorithm works for any $A$, its computational cost is dominated by the formation of products of $A$ with $n\times n_0$ matrices, and the only input parameter is a backward error tolerance. The algorithm can return a single matrix $e^{tA}B$ or a sequence $e^{t_kA}B$ on an equally spaced grid of points $t_k$. It uses the scaling part of the scaling and squaring method together with a truncated Taylor series approximation to the exponential. It determines the amount of scaling and the Taylor degree using the recent analysis of Al-Mohy and Higham [ (2009), pp. 970-989], which provides sharp truncation error bounds expressed in terms of the quantities $\|A^k\|^{1/k}$ for a few values of $k$, where the norms are estimated using a matrix norm estimator. Shifting and balancing are used as preprocessing steps to reduce the cost of the algorithm. Numerical experiments show that the algorithm performs in a numerically stable fashion across a wide range of problems, and analysis of rounding errors and of the conditioning of the problem provides theoretical support. Experimental comparisons with MATLAB codes based on Krylov subspace, Chebyshev polynomial, and Laguerre polynomial methods show the new algorithm to be sometimes much superior in terms of computational cost and accuracy. An important application of the algorithm is to exponential integrators for ordinary differential equations. It is shown that the sums of the form $\sum_{k=0}^p \varphi_k(A)u_k$ that arise in exponential integrators, where the $\varphi_k$ are related to the exponential function, can be expressed in terms of a single exponential of a matrix of dimension $n+p$ built by augmenting $A$ with additional rows and columns, and the algorithm of this paper can therefore be employed.

title: A Krylov projection method for systems of ODEs
abstract: In this paper we consider approximations of solutions of IVPs obtained through projections into Krylov subspaces. Numerical experiments on parabolic equations illustrate the performance of the method.

title: Approximation of Large-Scale Dynamical Systems
abstract: Preface Part I.  1.  Motivating examples Part II. Preliminaries: 3. Tools from matrix theory 4. Linear dynamical systems, Part 1 5. Linear dynamical systems, Part 2 6. Sylvester and Lyapunov equations Part III. SVD-based Approximation  7. Balancing and balanced approximations 8. Hankel-norm approximation 9. Special topics in SVD-based approximation methods Part IV. Krylov-based Approximation . Eigenvalue computations 11. Model reduction using Krylov methods Part V. SVD-Krylov Methods and Case Studies: 12. SVD-Krylov methods 13. Case studies 14. Epilogue 15. Problems Bibliography Index.

title: Approximation of the matrix exponential operator by a structure-preserving block Arnoldi-type method
abstract: 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) , Hochbruck and Lubich (1997) , Saad (1992) . 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) . A symplectic Krylov method to approximate exp(A)V was also proposed in Agoujil et al. (2012)  with [email protected]?R^2^n^x^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 ([email protected]?n). Our approach is based on block Krylov subspace methods that preserve Hamiltonian and skew-Hamiltonian structures.

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

abstract: 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.

title: Exponential time integration and Chebychev discretisation schemes for fast pricing of options
abstract: We consider exponential time differencing (ETD) schemes for the numerical pricing of options. Special treatments for the implementation of the boundary conditions that arise in finance are described. We show that only one explicit time step computation gives unconditional second order accuracy for European, Barrier and Butterfly spread options under both Black-Scholes geometric Brownian motion model and Merton's jump diffusion model with constant coefficients. In comparison, the commonly used Crank-Nicolson scheme is shown to be only conditionally stable due to lack of L"0-stability. Finally, we describe how the use of spectral spatial discretisation based on a Chebychev grid point concentration strategy gives fourth order accurate option prices for both the Black-Scholes and Merton's jump-diffusion model.

title: Acceleration Techniques for Approximating the Matrix Exponential Operator
abstract: In this paper we investigate some well-established and more recent methods that aim at approximating the vector $\exp(A)v$ when $A$ is a large symmetric negative semidefinite matrix, by efficiently combining subspace projections and spectral transformations. We show that some recently developed acceleration procedures may be restated as preconditioning techniques for the partial fraction expansion form of an approximating rational function. These new results allow us to devise a priori strategies to select the associated acceleration parameters; theoretical and numerical results are shown to justify these choices. Moreover, we provide a performance evaluation among several numerical approaches to approximate the action of the exponential of large matrices. Our numerical experiments provide a new, and in some cases, unexpected picture of the actual behavior of the discussed methods.

title: Residual, restarting and Richardson iteration for the matrix exponential
abstract: A well-known problem in computing some matrix functions iteratively is a lack of a clear, commonly accepted residual notion. An important matrix function for which this is the case is the matrix exponential. Assume, the matrix exponential of a given matrix times a given vector has to be computed. We interpret the sought after vector as a value of a vector function satisfying the linear system of ordinary differential equations (ODE), whose coefficients form the given matrix. The residual is then defined with respect to the initial-value problem for this ODE system. The residual introduced in this way can be seen as a backward error. We show how the residual can efficiently be computed within several iterative methods for the matrix exponential. This completely resolves the question of reliable stopping criteria for these methods. Furthermore, we show that the residual concept can be used to construct new residual-based iterative methods. In particular, a variant of the Richardson method for the new residual appears to provide an efficient way to restart Krylov subspace methods for evaluating the matrix exponential.

title: Krylov subspace approximation of eigenpairs and matrix functions in exact and computer arithmetic
abstract: Many researchers are now working on computing the product of a matrix function and a vector,using approximations in a Krylov subspace. We review our results on the analysis of one implemen-tation of that approach for symmetric matrices, which we call the Spectral Lanczos DecompositionMethod (SLDM).We have proved a general convergence estimate, relating SLDM error bounds to those obtainedthrough approximation of the matrix function by a part of its Chebyshev series. Thus, we arrivedat e ective estimates for matrix functions arising when solving parabolic, hyperbolic and ellipticpartial di erential equations. We concentrate on the parabolic case, where we obtain estimatesthat indicate superconvergence of SLDM. For this case a combination of SLDM and splittingmethods is also considered and some numerical results are presented.We implement our general estimates to obtain convergence bounds of Lanczos approximationsto eigenvalues in the internal part of the spectrum. Unlike Kaniel-Saad estimates, our estimatesare independent of the set of eigenvalues between the required one and the nearest spectrumbound.We consider an extension of our general estimate to the case of the simple Lanczos method(without reorthogonalization) in nite computer arithmetic which shows that for a moderatedimension of the Krylov subspace the results, proved for the exact arithmetic, are stable up toroundo .

title: Extended Krylov Subspaces: Approximation of the Matrix Square Root and Related Functions
abstract: We introduce an economical Gram--Schmidt orthogonalization on the extended Krylov subspace originated by actions of a symmetric matrix and its inverse. An error bound for a family of problems arising from the elliptic method of lines is derived. The bound shows that, for the same approximation quality, the diagonal variant of the extended subspaces requires about the square root of the dimension of the standard Krylov subspaces using only positive or negative matrix powers. An example of an application to the solution of a 2.5-D elliptic problem attests to the computational efficiency of the method for large-scale problems.

title: A Restarted Krylov Subspace Method for the Evaluation of Matrix Functions
abstract: We show how the Arnoldi algorithm for approximating a function of a matrix times a vector can be restarted in a manner analogous to restarted Krylov subspace methods for solving linear systems of equations. The resulting restarted algorithm reduces to other known algorithms for the reciprocal and the exponential functions. We further show that the restarted algorithm inherits the superlinear convergence property of its unrestarted counterpart for entire functions and present the results of numerical experiments.

title: RD-Rational Approximations of the Matrix Exponential
abstract: Restricted Denominator (RD) rational approximations to the matrix exponential operator are constructed by interpolation in points related to Krylov subspaces associated to a rational transform of the particular matrix considered. Convergence analysis are provided. Numerical experiments show the effectiveness of the proposed methods in applications involving discretizations of differential operators.

title: Analysis of some Krylov subspace approximations to the matrix exponential operator
abstract: In this note a theoretical analysis of some Krylov subspace approximations to the matrix exponential operation $\exp (A)v$ is presented, and a priori and a posteriors error estimates are established. Several such approximations are considered. The main idea of these techniquesis to approximately project the exponential operator onto a small Krylov subspace and to carry out the resulting small exponential matrix computation accurately. This general approach, which has been used with success in several applications, provides a systematic way of defining high-order explicit-type schemes for solving systems of ordinary differential equations or time-dependent partial differential equations.

title: Error Bounds for the Lanczos Methods for Approximating Matrix Exponentials
abstract: In this paper, we present new error bounds for the Lanczos method and the shift-and-invert Lanczos method for computing $e^{-\tau A} v$ for a large sparse symmetric positive semidefinite matrix $A$. Compared with the existing error analysis for these methods, our bounds relate the convergence to the condition numbers of the matrix that generates the Krylov subspace. In particular, we show that the Lanczos method will converge rapidly if the matrix $A$ is well-conditioned, regardless of what the norm of $\tau A$ is. Numerical examples are given to demonstrate the theoretical bounds.

title: 2-Norm Error Bounds and Estimates for Lanczos Approximations to Linear Systems and Rational Matrix Functions
abstract: The Lanczos process constructs a sequence of orthonormal vectors $v_m$ spanning a nested sequence of Krylov subspaces generated by a hermitian matrix $A$ and some starting vector $b$. In this paper we show how to cheaply recover a secondary Lanczos process starting at an arbitrary Lanczos vector $v_m$. This secondary process is then used to efficiently obtain computable error estimates and error bounds for the Lanczos approximations to the action of a rational matrix function on a vector. This includes, as a special case, the Lanczos approximation to the solution of a linear system $Ax = b$. Our approach uses the relation between the Lanczos process and quadrature as developed by Golub and Meurant. It is different from methods known so far because of its use of the secondary Lanczos process. With our approach, it is now possible in particular to efficiently obtain upper bounds for the error in the 2-norm, provided a lower bound on the smallest eigenvalue of $A$ is known. This holds in particular for a lar...

title: Stopping Criteria for Rational Matrix Functions of Hermitian and Symmetric Matrices
abstract: Building upon earlier work by Golub, Meurant, Strakos, and Tichý, we derive new a posteriori error bounds for Krylov subspace approximations to $f(A)b$, the action of a function $f$ of a real symmetric or complex Hermitian matrix $A$ on a vector $b$. To this purpose we assume that a rational function in partial fraction expansion form is used to approximate $f$, and the Krylov subspace approximations are obtained as linear combinations of Galerkin approximations to the individual terms in the partial fraction expansion. Our error estimates come at very low computational cost. In certain important special situations they can be shown to actually be lower bounds of the error. Our numerical results include experiments with the matrix exponential, as used in exponential integrators, and with the matrix sign function, as used in lattice quantum chromodynamics simulations, and demonstrate the accuracy of the estimates. The use of our error estimates within acceleration procedures is also discussed.

title: Analysis of Projection Methods for Rational Function Approximation to the Matrix Exponential
abstract: Krylov subspace methods for approximating the action of the matrix exponential exp(A) on a vector v are analyzed with A Hermitian and negative semidefinite. Our approach is based on approximating the exponential with the commonly employed diagonal Pade and Chebyshev rational functions, which yield a system of equations with a polynomial coefficient matrix. We derive optimality properties and error bounds for the convergence of a Galerkin-type approximation and of a computationally feasible and extensively used alternative. As complementary results, we theoretically justify the use of a popular a posteriori error estimate, and we provide upper bounds for the components of the solution vector. Our theoretical and numerical results show that this methodology may provide an appropriate framework to devise new strategies such as more powerful acceleration schemes.

title: The Carathéodory–Fejér Method for Real Rational Approximation
abstract: A “Caratheodory–Fejer method” is presented for near-best real rational approximation on intervals, based on the eigenvalue (or singular value) analysis of a Hankel matrix of Chebyshev coefficients. In approximation of a smooth function F, the CF approximant $R^{cf} $ frequently differs from the best approximation $R^ * $ by only one part in millions or billions. To account for this we show here under weak assumptions that if F is approximated on $[ - \varepsilon ,\varepsilon ]$, then as $\varepsilon \to 0$, $||F - R^ * || = O(\varepsilon ^{m + n + 1} )$ while $||R^{cf} - R^ * || = O(\varepsilon ^{3m + 2n + 3} )$. In contrast, the latter figure would be $O(\varepsilon ^{m + n + 2} )$ for the Chebyshev economization approximant of Maehly or the Chebyshev–Pade approximant of Gragg. It follows that as $\varepsilon \to 0$, best approximation error curves approach the real parts of $m + n + 1$-winding rational functions of constant modulus to within $O(\varepsilon ^{3m + 2n + 3} )$. Numerical examples are given...

title: Talbot quadratures and rational approximations
 
 
 
 

abstract: Many computational problems can be solved with the aid of contour integrals containing ez in the integrand: examples include inverse Laplace transforms, special functions, functions of matrices and operators, parabolic PDEs, and reaction-diffusion equations. One approach to the numerical quadrature of such integrals is to apply the trapezoid rule on a Hankel contour defined by a suitable change of variables. Optimal parameters for three classes of such contours have recently been derived: (a) parabolas, (b) hyperbolas, and (c) cotangent contours, following Talbot in 1979. The convergence rates for these optimized quadrature formulas are very fast: roughly O(3-N), where N is the number of sample points or function evaluations. On the other hand, convergence at a rate apparently about twice as fast, O(9.28903-N), can be achieved by using a different approach: best supremum-norm rational approximants to ez for z∈(–∞,0], following Cody, Meinardus and Varga in 1969. (All these rates are doubled in the case of self-adjoint operators or real integrands.) It is shown that the quadrature formulas can be interpreted as rational approximations and the rational approximations as quadrature formulas, and the strengths and weaknesses of the different approaches are discussed in the light of these connections. A MATLAB function is provided for computing Cody–Meinardus–Varga approximants by the method of Carathéodory–Fejér approximation.

title: On Analysis of Projection Methods for Rational Function Approximation to the Matrix Exponential
abstract: We revisit the analysis of projection methods for rational function approximation to the matrix exponential, and present a new result for establishing convergence of the Pade-plus-Krylov approximation when $||A||>1$.

title: Recent computational developments in Krylov subspace methods for linear systems
abstract: Many advances in the development of Krylov subspace methods for the iterative solution of linear systems during the last decade and a half are reviewed. These new developments include different versions of restarted, augmented, deflated, flexible, nested, and inexact methods. Also reviewed are methods specifically tailored to systems with special properties such as special forms of symmetry and those depending on one or more parameters.

title: BiCGStab(ℓ) for Families of Shifted Linear Systems
abstract: We consider a seed system Ax = b together with a shifted linear system of the formWe develop modifications of the BiCGStab(ℓ) method which allow to solve the seed and the shifted system at the expense of just the matrix-vector multiplications needed to solve Ax = b via BiCGStab(ℓ). On the shifted system, these modifications do not perform the corresponding BiCGStab(ℓ)-method, but we show, that in the case that A is positive real and σ ≥ 0, the resulting method is still a well-smoothed variant of BiCG. Numerical examples from an application arising in quantum chromodynamics are given to illustrate the efficiency of the method developed.

title: Restarted GMRES for Shifted Linear Systems
abstract: Shifted matrices, which differ by a multiple of the identity only, generate the same Krylov subspaces with respect to any fixed vector. This fact has been exploited in Lanczos-based methods like CG, QMR, and BiCG to simultaneously solve several shifted linear systems at the expense of only one matrix--vector multiplication per iteration. Here, we develop a variant of the restarted GMRES method exhibiting the same advantage and we investigate its convergence for positive real matrices in some detail. We apply our method to speed up "multiple masses" calculations arising in lattice gauge computations in quantum chromodynamics, one of the most time-consuming supercomputer applications.

title: A FLEXIBLE PRECONDITIONED ARNOLDI METHOD FOR SHIFTED LINEAR SYSTEMS
abstract: We are interested in the numerical solution of the large nonsymmetric shifted linear system, (A + αI)x = b, for many different values of the shift α in a wide range. We apply the Saad’s flexible preconditioning technique [14] to the solution of the shifted systems. Such flexible preconditioning with a few parameters could probably cover all the shifted systems with the shift in a wide range. Numerical experiments report the effectiveness of our approach on some problems.

title: Restarted Full Orthogonalization Method for Shifted Linear Systems
abstract: Restarted GMRES is known to be inefficient for solving shifted systems when the shifts are handled simultaneously. Variants have been proposed to enhance its performance. We show that another restarted method, restarted Full Orthogonalization Method (FOM), can effectively be employed. The total number of iterations of restarted FOM applied to all shifted systems simultaneously is the same as that obtained by applying restarted FOM to the shifted system with slowest convergence rate, while the computational cost grows only sub-linearly with the number of shifts. Numerical experiments are reported.

title: Krylov Subspace Recycling for Sequences of Shifted Linear Systems
abstract: We study the use of Krylov subspace recycling for the solution of a sequence of slowly-changing families of linear systems, where each family consists of shifted linear systems that differ in the coefficient matrix only by multiples of the identity. Our aim is to explore the simultaneous solution of each family of shifted systems within the framework of subspace recycling, using one augmented subspace to extract candidate solutions for all the shifted systems. The ideal method would use the same augmented subspace for all systems and have fixed storage requirements, independent of the number of shifted systems per family. We show that a method satisfying both requirements cannot exist in this framework. As an alternative, we introduce two schemes. One constructs a separate deflation space for each shifted system but solves each family of shifted systems simultaneously. The other builds only one recycled subspace and constructs approximate corrections to the solutions of the shifted systems at each cycle of the iterative linear solver while only minimizing the base system residual. At convergence of the base system solution, we apply the method recursively to the remaining unconverged systems. We present numerical examples involving systems arising in lattice quantum chromodynamics.

title: Deflated GMRES for systems with multiple shifts and multiple right-hand sides☆
abstract: We consider solution of multiply shifted systems of nonsymmetric linear equations, possibly also with multiple right-hand sides. First, for a single right-hand side, the matrix is shifted by several multiples of the identity. Such problems arise in a number of applications, including lattice quantum chromodynamics where the matrices are complex and non-Hermitian. Some Krylov iterative methods such as GMRES and BiCGStab have been used to solve multiply shifted systems for about the cost of solving just one system. Restarted GMRES can be improved by deflating eigenvalues for matrices that have a few small eigenvalues. We show that a particular deflated method, GMRES-DR, can be applied to multiply shifted systems. In quantum chromodynamics, it is common to have multiple right-hand sides with multiple shifts for each right-hand side. We develop a method that efficiently solves the multiple right-hand sides by using a deflated version of GMRES and yet keeps costs for all of the multiply shifted systems close to those for one shift. An example is given showing this can be extremely effective with a quantum chromodynamics matrix.

title: Restarted block-GMRES with deflation of eigenvalues
abstract: Block-GMRES is an iterative method for solving nonsymmetric systems of linear equations with multiple right-hand sides. Restarting may be needed, due to orthogonalization expense or limited storage. We discuss how restarting affects convergence and the role small eigenvalues play. Then a version of restarted block-GMRES that deflates eigenvalues is presented. It is demonstrated that deflation can be particularly important for block methods.

title: A Preconditioned and Shifted GMRES Algorithm for the PageRank Problem with Multiple Damping Factors
abstract: Google has become one of the most popular and successful search engines in recent years. Google's success can be attributed to its simple and elegant algorithm: PageRank. In practice, one often needs to solve the PageRank problem with multiple damping factors or with multiple damping factors and multiple personalization vectors. The conventional PageRank algorithm has to solve these problems one by one. The shifted GMRES$(m)$ algorithm can be used to solve them in the same search subspace. However, there are two disadvantages to this algorithm. The first is “near singularity,” and the second is “stagnation.” In this paper, we first present a modified and shifted GMRES$(m)$ algorithm to deal with the problem of near singularity. In order to overcome the drawback of stagnation and to improve convergence, we propose a polynomial preconditioner for the modified algorithm. We show that the resulting algorithm can circumvent the drawbacks of near singularity and stagnation that occur in its original counterpart...

title: Approximate Inverse Preconditioning for Shifted Linear Systems
abstract: In this paper we consider the problem of preconditioning symmetric positive definite matrices of the form Aα=A+αI where α>0. We discuss how to cheaply modify an existing sparse approximate inverse preconditioner for A in order to obtain a preconditioner for Aα. Numerical experiments illustrating the performance of the proposed approaches are presented.

title: Efficient Preconditioner Updates for Shifted Linear Systems
abstract: We present a technique for building effective and low cost preconditioners for sequences of shifted linear systems $(A + \alpha I) x_\alpha = b$, where $A$ is symmetric positive definite and $\alpha > 0$. This technique updates a preconditioner for $A$, available in the form of an $LDL^T$ factorization, by modifying only the nonzero entries of the $L$ factor in such a way that the resulting preconditioner mimics the diagonal of the shifted matrix and reproduces its overall behavior. This approach is supported by a theoretical analysis as well as by numerical experiments, showing that it works efficiently for a broad range of values of $\alpha$.

title: EFFICIENT PRECONDITIONING FOR SEQUENCES OF PARAMETRIC COMPLEX SYMMETRIC LINEAR SYSTEMS
abstract: Solution of sequences of complex symmetric linear systems of the form Ajxj = bj, j = 0; :::; s, Aj = A + jEj, A Hermitian, E0; :::; Es complex diagonal matrices and 0; s scalar complex parameters arise in a variety of challenging problems. This is the case of time dependent PDEs; lattice gauge computations in quantum chromodynamics; the Helmholtz equation; shift-and- invert and Jacobi{Davidson algorithms for large-scale eigenvalue calculations; problems in control theory and many others. If A is symmetric and has real entries then Aj is complex symmetric. The case A Hermitian positive semidenite, Re( j) 0 and such that the diagonal entries of Ej, j = 0; :::; s have nonnegative real part is considered here. Some strategies based on the update of incomplete factorizations of the matrix A and A 1 are introduced and analyzed. The numerical solution of sequences of algebraic linear systems from the discretization of the real and complex Helmholtz equation and of the diusion equation in a rectangle illustrate the performance of the proposed approaches.

title: On conjugate gradient type methods and polynomial preconditioners for a class of complex non-hermitian matrices
abstract: SummaryWe consider conjugate gradient type methods for the solution of large linear systemsA x=b with complex coefficient matrices of the typeA=T+iσI whereT is Hermitian and σ a real scalar. Three different conjugate gradient type approaches with iterates defined by a minimal residual property, a Galerkin type condition, and an Euclidean error minimization, respectively, are investigated. In particular, we propose numerically stable implementations based on the ideas behind Paige and Saunder's SYMMLQ and MINRES for real symmetric matrices and derive error bounds for all three methods. It is shown how the special shift structure ofA can be preserved by using polynomial preconditioning, and results on the optimal choice of the polynomial preconditioner are given. Also, we report on some numerical experiments for matrices arising from finite difference approximations to the complex Helmholtz equation.

title: Updating incomplete factorization preconditioners for shifted linear systems arising in a wind model
abstract: The efficiency of a finite element mass consistent model for wind field adjustment depends on the stability parameter @a which allows adjustment from a strictly horizontal wind to a pure vertical one. Each simulation with the wind model leads to the resolution of a linear system of equations, the matrix of which depends on a function @e(@a), i.e., (M+@eN)x"@e=b"@e, where M and N are constant, symmetric and positive definite matrices with the same sparsity pattern for a given level of discretization. The estimation of this parameter may be carried out by using genetic algorithms. This procedure requires the evaluation of a fitness function for each individual of the population defined in the searching space of @a, that is, the resolution of one linear system of equations for each value of @a.  Conjugate Gradient algorithm (PCG) is usually applied for the resolution of these types of linear systems due to its good convergence results. In order to solve this set of linear systems, we could either construct a different preconditioner for each of them or use a single preconditioner constructed from the first value of @e to solve all the systems. In this paper, an intermediate approach is proposed. An incomplete Cholesky factorization of matrix A"@e is constructed for the first linear system and it is updated for each @e at a low computational cost. Numerical experiments related to realistic wind field are presented in order to show the performance of the proposed preconditioning strategy.

title: A Flexible Krylov Solver for Shifted Systems with Application to Oscillatory Hydraulic Tomography
abstract: We discuss efficient solutions to systems of shifted linear systems arising in computations for oscillatory hydraulic tomography. The reconstruction of hydrogeological parameters such as hydraulic conductivity and specific storage using limited discrete measurements of pressure (head) obtained from sequential oscillatory pumping tests, leads to a nonlinear inverse problem. We tackle this using the quasi-linear geostatistical approach [Kitanidis, Water Resources Res., 31 (1995), pp. 2411--2419]. This method requires repeated solution of the forward (and adjoint) problem for multiple frequencies, for which we use flexible preconditioned Krylov subspace solvers specifically designed for shifted systems based on ideas in [Gu, Zhou, and Lin,  (2007), pp. 522--530]. The solvers allow the preconditioner to change at each iteration. We analyze the convergence of the solver and perform an error analysis when an iterative solver is used for inverting the preconditioner matrices. Finally, we appl...

title: Shift-Invert Arnoldi Approximation to the Toeplitz Matrix Exponential
abstract: The shift-invert Arnoldi method is employed to generate an orthonormal basis from the Krylov subspace corresponding to a real Toeplitz matrix and an initial vector. The vectors and recurrence coefficients produced by this method are exploited to approximate the Toeplitz matrix exponential. Toeplitz matrix inversion formula and rapid Toeplitz matrix-vector multiplications are utilized to lower the computational costs. For convergence analysis, a sufficient condition is established to guarantee that the error bound is independent of the norm of the matrix. Numerical results are given to demonstrate the efficiency of the method.

title: Evaluating matrix functions for exponential integrators via Carathéodory-Fejér approximation and contour integrals
abstract: Among the fastest methods for solving stiff PDE are exponential integrators, which require the evaluation of , where is a negative semidenite matrix and is the exponential function or one of the related i functionsi such as . Building on previous work by Trefethen and Gutknecht, Minchev, and Lu, we propose two methods for the fast evaluation of that are especially useful when shifted systems can be solved efciently , e.g. by a sparse direct solver. The rst method is based on best rational approximations to on the negative real axis computed via the CarathOodory-FejOr procedure. Rather than using optimal poles we approximate the functions in a set of common poles, which speeds up typical computations by a factor of to . The second method is an application of the trapezoid rule on a Talbot-type contour.

title: On the use of matrix functions for fractional partial differential equations
abstract: The main focus of this paper is the solution of some partial differential equations of fractional order. Promising methods based on matrix functions are taken in consideration. The features of different approaches are discussed and compared with results provided by classical convolution quadrature rules. By means of numerical experiments accuracy and performance are examined.

title: Extended numerical computations on the “1/9” conjecture in rational approximation theory
abstract: The behavior of the constants λn,n(e−x), denoting the errors of best uniform approximation to e−z on the interval [0,+∞) by real rational functions having numerator and denominator polynomials of degree at most n, has generated much recent interest in the approximation theory literature. Based on high-precision calculations, we present here the table of constants {λn,n (e−x)} n 30 =0, rounded to forty significant digits, and we discuss their significance to related conjectures in this area.

title: Asymptotics and Super Asymptotics for Best Rational Approximation Error Norms to the Exponential Function (The ‘1/9’ Problem) by the Carathéodory-Fejér Method
abstract: Let E nbe the error norm of the best L ∞ rational approximation of degree n to the exponential function exp(-t) on [0,∞). Grounds are given for setting the conjectured limit E n / q n →2q 1/2 when n →∞, where q is the known constant ‘l/9’= 1/9.2890254919208189187554494359517450610316948677…, based on the singular values and functions of the relevant Henkel operator (Caratheodory-Fejer’s method).

title: A thick-restarted block Arnoldi algorithm with modified Ritz vectors for large eigenproblems
abstract: The block Arnoldi method is one of the most commonly used techniques for large eigenproblems. In this paper, we exploit certain modified Ritz vectors to take the place of Ritz vectors in the thick-restarted block Arnoldi algorithm, and propose a modified thick-restarted block Arnoldi algorithm for large eigenproblems. We then consider how to periodically combine the refined subspace iterative method with the modified thick-restarting block Arnoldi algorithm for computing a few dominant eigenpairs of a large matrix. The resulting algorithm is called a Subspace-Block Arnoldi algorithm. Numerical experiments show the efficiency of our new algorithms.

title: Numerical Approximation of Partial Differential Equations
abstract: This is the softcover reprint of the very popular hardcover edition. This book deals with the numerical approximation of partial differential equations. Its scope is to provide a thorough illustration of numerical methods, carry out their stability and convergence analysis, derive error bounds, and discuss the algorithmic aspects relative to their implementation. A sound balancing of theoretical analysis, description of algorithms and discussion of applications is one of its main features. Many kinds of problems are addressed. A comprehensive theory of Galerkin method and its variants, as well as that of collocation methods, are developed for the spatial discretization. These theories are then specified to two numerical subspace realizations of remarkable interest: the finite element method and the spectral method.

title: Extended Krylov subspace for parameter dependent systems
abstract: The Extended Krylov subspace has recently received considerable attention as a powerful tool for matrix function evaluations and other problems involving large matrices. In this paper we show that this space has a great potential within projection-type methods for effectively solving several other important large-scale algebraic problems: we focus on the solution of shifted systems and of more general parameter-dependent matrix equations, and on the approximation of the transfer function by projection. Numerical experiments stemming from real applications show the effectiveness of the approach.

title: Shift-invert Lanczos method for the symmetric positive semidefinite Toeplitz matrix exponential
abstract: The Lanczos method with shift-invert technique is exploited to approximate the symmetric positive semidefinite Toeplitz matrix exponential. The complexity is lowered by the Gohberg–Semencul formula and the fast Fourier transform. Application to the numerical solution of an integral equation is studied. Numerical experiments are carried out to demonstrate the effectiveness of the proposed method.

title: Inexact Matrix-Vector Products in Krylov Methods for Solving Linear Systems: A Relaxation Strategy
abstract: Embedded iterative linear solvers are being used more and more often in linear algebra. An important issue is how to tune the level of accuracy of the inner solver to guarantee the convergence of the outer solver at the best global cost. As a first step towards the challenging goal of controlling embedded linear solvers, inexact Krylov methods are used as a model of inner-outer iterations with external Krylov scheme. This paper experimentally shows that Krylov methods for solving linear systems can still perform very well in the presence of carefully monitored inexact matrix-vector products. This surprising behavior of inexact Krylov methods, as opposed to Newton-like methods, is investigated in detail, and potentially important applications are mentioned. A new relaxation strategy for the inner accuracy is proposed for Krylov methods with inexact matrix-vector products; its efficiency is supported by a wide range of numerical experiments on different algorithms and contrasted against other potential approaches.

title: Theory of Inexact Krylov Subspace Methods and Applications to Scientific Computing
abstract: We provide a general framework for the understanding of inexact Krylov subspace methods for the solution of symmetric and nonsymmetric linear systems of equations, as well as for certain eigenvalue calculations. This framework allows us to explain the empirical results reported in a series of CERFACS technical reports by Bouras, Fraysse, and Giraud in 2000. Furthermore, assuming exact arithmetic, our analysis can be used to produce computable criteria to bound the inexactness of the matrix-vector multiplication in such a way as to maintain the convergence of the Krylov subspace method. The theory developed is applied to several problems including the solution of Schur complement systems, linear systems which depend on a parameter, and eigenvalue problems. Numerical experiments for some of these scientific applications are reported.

title: The Scaling and Squaring Method for the Matrix Exponential Revisited
abstract: The scaling and squaring method is the most widely used method for computing the matrix exponential, not least because it is the method implemented in the MATLAB function expm. The method scales the matrix by a power of 2 to reduce the norm to order 1, computes a Pade approximant to the matrix exponential, and then repeatedly squares to undo the effect of the scaling. We give a new backward error analysis of the method (in exact arithmetic) that employs sharp bounds for the truncation errors and leads to an implementation of essentially optimal efficiency. We also give a new rounding error analysis that shows the computed Pade approximant of the scaled matrix to be highly accurate. For IEEE double precision arithmetic the best choice of degree of Pade approximant turns out to be 13, rather than the 6 or 8 used by previous authors. Our implementation of the scaling and squaring method always requires at least two fewer matrix multiplications than the expm function in MATLAB 7.0 when the matrix norm exceeds 1, which can amount to a 37% saving in the number of multiplications, and it is typically more accurate, owing to the fewer required squarings. We also investigate a different scaling and squaring algorithm proposed by Najfeld and Havel that employs a Pade approximation to the function $x \coth(x)$. This method is found to be essentially a variation of the standard one with weaker supporting error analysis.

title: Expokit: a software package for computing matrix exponentials
abstract: Expokit provides a set of routines aimed at computing matrix exponentials. More precisely, it computes either a small matrix exponential in full, the action of a large sparse matrix exponential on an operand vector, or the solution of a system of linear OBEs with constant inhomogeneity. The backbone of the sparse routines consists of matrix-free Krylov subspace projection methods (Arnoldi and Lanczos processes), and that is why the toolkit is capable of coping with sparse matrices of large dimension. The software handles real and complex matrices and provides specific routines for symmetric and Hermitian matrices. The computation of matrix exponentials is a numerical issue of critical importance in the area of Markov chains and furthermore, the computed solution is subject to probabilistic constraints. In addition to addressing general matrix exponentials, a distinct attention is assigned to the computation of transient states of Markov chains.

title: Algorithm 919: A Krylov Subspace Algorithm for Evaluating the ϕ-Functions Appearing in Exponential Integrators
abstract: We develop an algorithm for computing the solution of a large system of linear ordinary differential equations (ODEs) with polynomial inhomogeneity. This is equivalent to computing the action of a certain matrix function on the vector representing the initial condition. The matrix function is a linear combination of the matrix exponential and other functions related to the exponential (the so-called ϕ-functions). Such computations are the major computational burden in the implementation of exponential integrators, which can solve general ODEs. Our approach is to compute the action of the matrix function by constructing a Krylov subspace using Arnoldi or Lanczos iteration and projecting the function on this subspace. This is combined with time-stepping to prevent the Krylov subspace from growing too large. The algorithm is fully adaptive: it varies both the size of the time steps and the dimension of the Krylov subspace to reach the required accuracy. We implement this algorithm in the matlab function phipm and we give instructions on how to obtain and use this function. Various numerical experiments show that the phipm function is often significantly more efficient than the state-of-the-art.

title: Deflated Restarting for Matrix Functions
abstract: We investigate an acceleration technique for restarted Krylov subspace methods for computing the action of a function of a large sparse matrix on a vector. Its effect is to ultimately deflate a specific invariant subspace of the matrix which most impedes the convergence of the restarted approximation process. An approximation to the subspace to be deflated is successively refined in the course of the underlying restarted Arnoldi process by extracting Ritz vectors and using those closest to the spectral region of interest as exact shifts. The approximation is constructed with the help of a generalization of Krylov decompositions to linearly dependent vectors. A description of the restarted process as a successive interpolation scheme at Ritz values is given in which the exact shifts are replaced with improved approximations of eigenvalues in each restart cycle. Numerical experiments demonstrate the efficacy of the approach.

title: A second-order accurate numerical approximation for the fractional diffusion equation
abstract: Fractional order diffusion equations are generalizations of classical diffusion equations, treating super-diffusive flow processes. In this paper, we examine a practical numerical method which is second-order accurate in time and in space to solve a class of initial-boundary value fractional diffusive equations with variable coefficients on a finite domain. An approach based on the classical Crank-Nicholson method combined with spatial extrapolation is used to obtain temporally and spatially second-order accurate numerical estimates. Stability, consistency, and (therefore) convergence of the method are examined. It is shown that the fractional Crank-Nicholson method based on the shifted Grunwald formula is unconditionally stable. A numerical example is presented and compared with the exact analytical solution for its order of convergence.

title: Finite di$erence approximations for fractional advection-dispersion &ow equations
abstract: Fractional advection–dispersion equations are used in groundwater hydrology to model the transport of passive tracers carried by &uid &ow in a porous medium.In this paper we develop practical numerical methods to solve one dimensional fractional advection–dispersion equations with variable coe4cients on a 5nite domain.The practical application of these results is illustrated by modeling a radial &ow problem.Use of the fractional derivative allows the model equations to capture the early arrival of tracer observed at a 5eld site.

title: A direct O(N log2 N) finite difference method for fractional diffusion equations
abstract: Fractional diffusion equations model phenomena exhibiting anomalous diffusion that can not be modeled accurately by the second-order diffusion equations. Because of the nonlocal property of fractional differential operators, the numerical methods have full coefficient matrices which require storage of O(N^2) and computational cost of O(N^3) where N is the number of grid points. In this paper we develop a fast finite difference method for fractional diffusion equations, which only requires storage of O(N) and computational cost of O(Nlog^2N) while retaining the same accuracy and approximation property as the regular finite difference method. Numerical experiments are presented to show the utility of the method.

title: The Australian National University TR-CS-95-07 Conjugate Gradient Methods for Toeplitz Systems
abstract: In this expository paper, we survey some of the latest developments on using preconditioned conjugate gradient methods for solving Toeplitz systems. One of the main results is that the complexity of solving a large class of n-by-n Toeplitz systems is reduced to O(n logn) operations as compared to O(n log n) operations required by fast direct Toeplitz solvers. Di erent preconditioners proposed for Toeplitz systems are reviewed. Applications to Toeplitz-related systems arising from partial di erential equations, queueing networks, signal and image processing, integral equations, and time series analysis are given.

title: An Introduction to Iterative Toeplitz Solvers
abstract: Preface 1.  Circulant preconditioners 3. Unified treatment from kernels 4. Ill-conditioned Toeplitz systems 5. Block Toeplitz systems A. M-files used in the book Bibliography.

title: A CG Method for Multiple Right Hand Sides and Multiple Shifts in Lattice QCD Calculations
abstract: 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 QCD simulations within the Rational Hybrid Monte Carlo algorithm. We give numerical evidence that our method is superior to applying other iterative methods to each of the systems individually as well as, in some cases, to shifted or block Krylov subspace methods.

title: Application of block Krylov subspace algorithms to the Wilson-Dirac equation with multiple right-hand sides in lattice QCD
abstract: 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.