1 Introduction
Overview on prior work.
The matrix exponential and associated functions play a crucial role in some numerical methods for solving systems of differential equations. In practice this means that the vector for a given matrix and a given vector , representing the time propagation for a linear initial value problem, is to be approximated. Similarly, the associated functions (see (2.2) below) conform to solutions of certain inhomogeneous differential equations. In particular, evaluation of functions is used in exponential integrators HO10 .
If the matrix is sparse and large, approximation of the action of these matrix functions in the class of Krylov subspaces is a general and wellestablished technique. For the matrix exponential and functions this goes back to early works in the field of chemical physics NR83 ; PL86 , parabolic problems GS92 , some nonlinear problems FTDR89 , etc. The case of a symmetric or skewHermitian matrix is the most prominent one. Krylov approximations of the matrix exponential were early studied for the symmetric case in DK89 ; DK91 ; Sa92 , and together with functions in a more general setting HLS98 ; HL97 .
Concerning different approaches for the numerical approximation of the matrix exponential see ML03 . In Sa92
it is shown for the symmetric case that the Krylov approximation is equivalent to interpolation of the exponential function at associated Ritz values. This automatically results in a nearbest approximation among other choices of interpolation nodes, see also
DK89 ; SL96 and future works BR09 with similar results for the nonsymmetric case. For other polynomial approaches approximating the matrix exponential we mention truncated Taylor series MH11 (and many works well in advance), Chebychev polynomial interpolation TK84 , or the Leja method CKOR16 .In general, Krylov approximation (or other polynomial approximations) result in an accurate approximation if the time step in is sufficiently small or the dimension of the Krylov subspace (i.e., the degree of the approximating matrix polynomial) is sufficiently large, see for instance HL97 . The dimension of the Krylov subspace is limited in practice, and large time steps require a restart of the iteration generating the Krylov basis. A larger time step can be split into smaller substeps for which the Krylov approximation can be applied in a nested way. Such a restart strategy in the sense of a time integrator was already exploited in PL86 . In particular we refer to the EXPOKIT package Si98 . Similar ideas can be applied for the evaluation of functions HLS98 ; NW12 .
In practice, a posteriori error estimates are used to choose a proper Krylov dimension or proper (adaptive) substeps if the method is restarted as a time integrator. Different approaches for a posteriori error estimation make use of a series expansion for the error given Sa92 ; Si98 or use a formulation via the defect (also called residual) of the Krylov approximation DGK98 ; CM97 ; BGH13 . Further a priori as well as a posteriori error estimates are given in MN01 ; Lu08 ; DMR09 ; BR09 ; JL15 ; WY17 ; JAK19 . Restarting via substeps based on different choices of error estimates is further discussed in JAK19 . A restart with substeps together with a strategy to choose the Krylov dimension in terms of computational cost was presented in NW12 ; BK19 . For various other approaches for restarting (without adapting the time step) we refer to CM97 ; EE06 ; Ta07 ; Ni07 ; AEEG08 ; EEG11 ; BGH13 ; Schwe15 .
The influence of roundoff errors on the construction of the Krylov basis in floating point arithmetic was early studied for the symmetric case in Pa76 ; Par98 . The orthogonalization procedure can behave numerically unstable, typically due to a loss of orthogonality. Nevertheless, the nearbest approximation property and related a priori convergence results are not critically affected DK91 ; DGK98 . Following DGK98 , in the symmetric case the defect obtained in floating point arithmetic results in numerically stable error estimates.
Beside the polynomial Krylov method, further studies are devoted to the approximation of matrix functions using so called extended Krylov subspaces DK98 ; KS10 ; GG13 , rational Krylov subspaces MN04 ; EH06 ; Gu10 , or polynomial Krylov subspaces with a harmonic Ritz approach HH05 ; Schwe15 ; WZX16 .
Overview on results presented here.
In Section 2 we introduce the problem setting and recapitulate basic properties of Krylov subspaces.
In Section 3 we introduce the defect associated with Krylov approximations to functions, including the exponential function as the basic case. Our approach for the defect is different from WZX16 and is based on an inhomogeneous differential equation for the approximation error. This is used in Theorem 3.1 to obtain an integral representation of the error, also taking effects of floating point arithmetic into account. ^{1}^{1}1Cf. DGK98 for the case of the matrix exponential. In contrast to previous works (DGK98 ; JAK19 ), this result is extended to functions here. Theorem 3.1 also includes an a priori upper bound on the error norm based on an integral of the defect norm.
This upper bound is further analyzed in Section 4 to obtain computable a posteriori bounds, in particular a new a posteriori bound (Theorem 4.1). We also study the accuracy of our and other existing defectbased bounds JAK19 with respect to spectral properties of the Krylov Hessenberg matrix (the representation of in the orthogonal Krylov basis). To this end we use properties of divided differences including a new asymptotic expansion for such given in Appendix C. In Subsection 4.1 we recapitulate error estimates based on a quadrature estimate of the defect norm integral, e.g. the generalized residual estimate HLS98 . We also discuss cases for which the defect norm behaves oscillatory and reliable quadrature estimates may be difficult to obtain. In Subsection 4.2 we specify a stopping criterion for the socalled lucky breakdown in floating point arithmetic which is justified by our a posteriori error bounds.
In Section 5 we illustrate our results via numerical experiments. Here we confine ourselves to the exponential function. Particular application problems where functions are used will be documented elsewhere.
2 Problem statement and Krylov approximation
We discuss the approximation via Krylov techniques for evaluation of the matrix exponential, and in particular of the associated functions, for a step size and matrix applied to an initial vector . Here,
(2.1) 
The matrix exponential is the solution of the differential equation
The associated functions are given by
(2.2) 
This includes the case . The matrix functions (2.1) and (2.2) are defined according to their scalar counterparts. The following definitions of are equivalent to (2.2): For we have , and
(2.3) 
(See also (Hi08, , Subsection 10.7.4).) The function () is the solution of an inhomogeneous differential equation of the form
(2.4) 
see for instance NW12 . This follows from (2.2),
The functions appear for instance in the field of exponential integrators, see for instance HO10 .
For the case of being a large and sparse matrix, e.g., the spatial discretization of a partial differential operator using a localized basis, Krylov subspace techniques are commonly used to approximate (2.2) in an efficient way.
Notation and properties of Krylov subspaces.
^{2}^{2}2In the sequel, denotes the th unit vector in or , respectively.We briefly recapitulate the usual notation and properties of standard Krylov subspaces, see for instance Sa03 . For a given matrix , a starting vector and Krylov dimension , the Krylov subspace is given by
Let represent the orthonormal basis of with respect to the Hermitian inner product, constructed by the Arnoldi method and satisfying . Its first column is given by with . Here, the matrix
is upper Hessenberg. We further use the notation , and for the th column of , with and .
The Arnoldi decomposition (in exact arithmetic) can be expressed in matrix form,
(2.5) 
Remark 1
The numerical range plays a role in our analysis. Note that (see (A.1)).
Remark 2
The case occurs if is an invariant subspace of , whence the Krylov approximation given in (2.9) below is exact. This exceptional case is referred to as a lucky breakdown. In general we assume that no lucky breakdown occurs, whence the lower subdiagonal entries of are real and positive, for , and .
For the special case of a Hermitian or skewHermitian matrix the Arnoldi iteration simplifies to a threeterm recurrence, the socalled Lanczos iteration. This case will be addressed in Remark 4 below.
Krylov subspaces in floating point arithmetic.
We proceed with some results for the Arnoldi decomposition in computer arithmetic, assuming complex floating point arithmetic with a relative machine precision , see also Hi02 . For practical implementation different variants of the Arnoldi procedure exist, using different ways for the orthogonalization of the Krylov basis. These are based on classical GramSchmidt, modified GramSchmidt, the Householder algorithm, the Givens algorithm, or variants of GramSchmidt with reorthogonalization (see also (Sa03, , Algorithm 6.1–6.3) and others). We refer to BLR00 and references therein for an overview on the stability properties of these different variants.
In the sequel the notation , , etc., will again be used for the result of the Arnoldi method in floating point arithmetic. We now accordingly adapt some statements formulated in the previous paragraph. By construction, remains to be upper Hessenberg with positive lower subdiagonal entries. Assuming floating point arithmetic we use the notation for a perturbation of the Arnoldi decomposition (2.5) caused by roundoff, i.e.,
(2.6) 
An upper norm bound for was first introduced in Pa76 for the Lanczos iteration in real arithmetic. For different variants of the Arnoldi or Lanczos iteration this is discussed in Ze03 and others. We assume is bounded by a constant which can depend on and in a moderate way and is sufficiently small in a typical setting,
(2.7a)  
We further assume that the normalization of the columns of is accurate, in particular that the th basis vector is normalized correctly up roundoff with a sufficiently small constant (see e.g. (Pa76, , eq. (14))),  
(2.7b)  
Concerning which represents an orthogonal basis in exact arithmetic, numerical loss of orthogonality has been wellstudied. Loss of orthogonality can be significant (see for instance Par98 ; BLR00 and others), depending on the starting vector . Reorthogonalization schemes or orthogonalization via Householder or Givens algorithm can be used to obtain orthogonality of on a sufficiently accurate level. 
The numerical range of obtained in floating point arithmetic (see (2.6)) can be characterized as
(2.7c) 
with being the neighborhood of in with a distance . With the assumption that is sufficiently close to orthogonal (e.g., semiorthogonal Si84 ), the constant in (2.7c) (which also depends on and problem sizes) can be shown to be moderatesized. Further details on this aspect are given in Appendix A.
Krylov approximation of functions.
^{3}^{3}3Remark concerning notation: ’’ objects live in , and ’’ objects live in .Let , and be the result of the Arnoldi method in floating point arithmetic for as described above. For a timestep and the vector can be approximated in the Krylov subspace by the Krylov propagator
(2.8a)  
The special case reads  
(2.8b) 
We remark that the smalldimensional problem , typically with , can be evaluated cheaply by standard methods. In the sequel we denote
(2.9) 
For the small dimensional problem solves the differential equation
(2.10) 
For later use we introduce the notation
(2.11a)  
which for and according to (2.4) satisfies the differential equation  
(2.11b) 
Remark 3
Although we take rounding effects in the Arnoldi decomposition into account, we do not give a full study of roundoff errors at this point. Roundoff errors in substeps such as the evaluation of or the matrixvector multiplication will be ignored. We refer to Hi02 for a more general study of these effects.
Remark 4
In the special cases or for a Hermitian matrix (with being skewHermitian in the latter case) the orthogonalization of the Krylov basis of simplifies to a threeterm recursion, the socalled Lanczos method. In the skewHermitian case () the Krylov propagator (2.8a) can be evaluated by , i.e., we approximate the function in the Krylov subspace . The advantage is a cheaper computation of the Krylov subspace in terms of computational cost and better conservation of geometric properties. For details we refer to the notation as introduced in JAK19 , with and a Hermitian matrix for the skewHermitian case.
The error of the Krylov propagator.
We denote the error of the Krylov propagator given in (2.9) by
(2.12) 
We are further interested in computable a posteriori estimates for the error norm , which in the best case can be proven to be upper bounds on the error norm . Norm estimates of the error (2.12) can be used in practice to stop the Krylov iteration after steps if satisfies (2.13) below, or to restrict the timestep to obtain an accurate approximation and restart the method with the remaining time. For details on the total error with this restarting approach see also Si98 ; JAK19 .
A prominent task is to test if the error norm per unit step is bounded by a tolerance ,
(2.13) 
In case of being an upper bound on the error norm, this results in a reliable bound on the error norm (2.13).
3 An integral representation for the error of the Krylov propagator
We proceed with discussing the error of the Krlyov propagator. To this end we first define its scalar defect by
(3.1a)  
and the defect integral by^{4}^{4}4 This and the result of Theorem 3.1 remain valid for the case .  
(3.1b) 
Theorem 3.1
Proof

For the exact matrix function we use the notation
For the Krylov propagator we denote
(see (2.9)), and we also define

For , i.e., and , according to (2.10) we have
Local error representation in terms of the defect.

With , and , the representation (3.2a) implies the upper bound
(3.6) With the defect integral defined in (3.1b) we obtain the first term in (3.2b). For the second integral term (with ) we use the upper bound
(3.7) 
For we apply the integral representation due to (2.3) for to obtain the norm bound
(3.8) 
For we obtain (3.8) in a direct way.
Combining (3.7) with (3.8) and denoting we obtain

Remark 5
The error norm of the Krylov propagator scales with and in a natural way. ^{5}^{5}5Taking the maximum in the definition of and is necessary to cover the case . For the special case the upper norm bound given in Theorem 3.1 can be adapted to scale with . It is well known that
see for instance (Hi08, , Theorem 10.11). Problems with can be arbitrary illconditioned and difficult to solve with proper accuracy. (For further results on the stability of the matrix exponential see also ML03 ; Lo77 .) We will not further discuss problems with and assume . We refer to the case as the dissipative case, with .
For the dissipative case with the error bound (3.2b) from Theorem 3.1 reads
(3.9) 
The dissipative behavior of carries over to the Krylov propagator up to a perturbation which depends on roundoff errors, including the loss of orthogonality of . In terms of the numerical range , with we have , for a constant depending on roundoff effects (2.7c). Thus, and .
Our aim is to construct an upper norm bound for the error per unit step (2.13) via (3.9). Let the tolerance be given and be a respective time step for (2.13). Then the roundoff error terms in (3.9) are negligible if
(3.10) 
Concerning the constants , and see (2.7). We recapitulate that and given in (2.7a) and (2.7b) can be considered to be small enough in a standard Krylov setting. The constant can be larger in the case of a loss of orthogonality of the Krylov subspace, which can however be avoided at the cost of additional computational effort. The constant only appears as an exponential prefactor for the roundoff term in (3.10) and is less critical compared to and .
With the previous observation on the roundoff errors taken into account in (3.9) we consider the following upper bound to be stable in computer arithmetic in accordance to a proper value of , see (3.10).
Corollary 1
For the case and with the assumption that roundoff error is negligible, the error of the Krylov propagator is bounded by the defect integral ,
Note that the defect norm cannot be integrated exactly in general. This point will further be studied in the sequel.
Representing the defect in terms of divided differences.
Divided differences play an essential role in this work. We use the notation
for the divided differences of a function over the nodes . (This is to be understood in the confluent sense for the case of multiple nodes , see for instance (Hi08, , Section B.16).)
Theorem 3.2 (see for instance Cm97 )
Let be an upper Hessenberg matrix with positive secondary diagonal entries, for
, and eigenvalues
. Let be an analytic function for which is well defined. Then,with .
For we will also make use of the following result. ^{6}^{6}6Theorem 3.3 can be generalized to the case with , see (MH11, , Theorem 2.1). The case is sufficient for our purpose.
Theorem 3.3 (Corollary 1 in Si98 )
(Expressing functions via dilated functions.) For ,
with
The matrix in Theorem 3.3 is block triangular with eigenvalues equal to those of and . Therefore, , with as an eigenvalue of multiplicity (at least). In our context, is upper Hessenberg with a positive lower secondary diagonal and . In accordance with Theorem 3.2 the result of Theorem 3.3 holds for divided differences in a similar manner,
With Theorem 3.2 and 3.3 the following equivalent formulations can be used the rewrite the scalar defect defined in (3.1a).
Corollary 2
We remark that the eigenvalues of the Krylov Hessenberg matrix are also referred to as Ritz values (of ) in the literature.
4 Computable a posteriori error bounds for the Krylov propagator
The following two propositions are used for the proof of Theorem 4.1 below.^{8}^{8}8 We use the notation introduced in the previous sections.
Proposition 1
For arbitrary nodes and ,
Proof
See Appendix B.
Proposition 2 (Lemma including eq. (5.1.1) in Mnp84 )
For arbitrary nodes ,
Proof
See Appendix B.
We now derive upper bounds for the error via its representation by the defect integral (3.1b).
Theorem 4.1
Let , , and assume that roundoff errors are sufficiently small (see Corollary 1). For the eigenvalues of we write , . An upper bound on the error norm is given by
(4.1) 
Proof
For the case of having real eigenvalues, the assertion of Theorem 4.1 can be reformulated in the following way (see (JAK19, , Proposition 6)).
Corollary 3
As a further corollary we formulate an upper bound on the error norm which is cheaper to evaluate compared to the bound from Theorem 4.1 but may be less tight. Using the Mean Value Theorem, (Hi08, , eq. (B.26)) or (Bo05, , eq. (44)), for the divided differences in Theorem 4.1, eq. (4.1) we obtain the following result which corresponds to (JAK19, , Theorem 1 and 2). For the exponential of a skewHermitian matrix a similar error estimate has been used in KBC05 and is based on ideas of PL86 with some lack of theory.
Corollary 4
Let , , and assume that roundoff errors are sufficiently small (see Corollary 1). Let for and for and eigenvalues of . An upper bound on the error norm is given by
Comments
There are no comments yet.