A study of defect-based error estimates for the Krylov approximation of $\varphi$-functions

Prior recent work, devoted to the study of polynomial Krylov techniques for the approximation of the action of the matrix exponential ${\rm e}^{tA}v$, is extended to the case of associated $\varphi$-functions (which occur within the class of exponential integrators). In particular, a~posteriori error bounds and estimates, based on the notion of the defect (residual) of the Krylov approximation are considered. Computable error bounds and estimates are discussed and analyzed. This includes a new error bound which favorably compares to existing error bounds in specific cases. The accuracy of various error bounds is characterized in relation to corresponding Ritz values of $A$. Ritz values yield properties of the spectrum of $A$ (specific properties are known a~priori, e.g. for Hermitian or skew-Hermitian matrices) in relation to the actual starting vector $v$ and can be computed. This gives theoretical results together with criteria to quantify the achieved accuracy on the run. For other existing error estimates the reliability and performance is studied by similar techniques. Effects of finite precision (floating point arithmetic) are also taken into account.


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 e tA v for a time step t, for a given matrix A and a given vector v, 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 [27].
If the matrix A is sparse and large, approximation of the action of these matrix functions in the class of Krylov subspaces is a general and well-established technique. For the matrix exponential and ϕ-functions this goes back to early works in the field of chemical physics [39,44], parabolic problems [20], some nonlinear problems [18], etc. The case of a symmetric or skew-Hermitian matrix A is the most prominent one. Krylov approximations of the matrix exponential were early studied for the symmetric case in [12,13,46], and together with ϕ-functions in a more general setting [28,26].
Concerning different approaches for the numerical approximation of the matrix exponential see [36]. In [46] 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 near-best approximation among other choices of interpolation nodes, see also [12,52] and future works [3] with similar results for the non-symmetric case. For other polynomial approaches approximating the matrix exponential we mention truncated Taylor series [2] (and many works well in advance), Chebychev polynomial interpolation [54], or the Leja method [8].
In general, Krylov approximation (or other polynomial approximations) result in an accurate approximation if the time step t in e tA v 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 [26]. 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 t 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 [44]. In particular we refer to the EXPOKIT package [49]. Similar ideas can be applied for the evaluation of ϕ-functions [28,41].
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 [46,49] or use a formulation via the defect (also called residual) of the Krylov approximation [11,9,5]. Further a priori as well as a posteriori error estimates are given in [37,34,10,3,31,56,30]. Restarting via substeps based on different choices of error estimates is further discussed in [30]. A restart with substeps together with a strategy to choose the Krylov dimension in terms of computational cost was presented in [41,6]. For various other approaches for restarting (without adapting the time step) we refer to [9,15,53,40,1,16,5,48].
The influence of round-off errors on the construction of the Krylov basis in floating point arithmetic was early studied for the symmetric case in [43,45]. The orthogonalization procedure can behave numerically unstable, typically due to a loss of orthogonality. Nevertheless, the near-best approximation property and related a priori convergence results are not critically affected [13,11]. Following [11], in the symmetric case the defect obtained in floating point arithmetic results in numerically stable error estimates.
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 [57] and is based on an inhomogeneous differential equation for the approximation error. This is used in Theorem 1 to obtain an integral representation of the error, also taking effects of floating point arithmetic into account. 1 In contrast to previous works ( [11,30]), this result is extended to ϕ-functions here. Theorem 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). We also study the accuracy of our and other existing defect-based bounds [30] with respect to spectral properties of the Krylov Hessenberg matrix (the representation of A 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 [28]. 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 so-called 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.

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 t > 0 and matrix A ∈ C n×n applied to an initial vector v ∈ C n . Here, The matrix exponential u(t) = e tA v is the solution of the differential equation The associated ϕ-functions are given by This includes the case ϕ 0 = exp. The matrix functions (2.1) and (2.2) are defined according to their scalar counterparts. The following definitions of ϕ p are equivalent to (2.2): For z ∈ C we have ϕ 0 (z) = e z , and (See also [24,Subsection 10.7.4].) The function w p (t) = t p ϕ p (tA)v (p ∈ N) is the solution of an inhomogeneous differential equation of the form see for instance [41]. This follows from (2.2), The ϕ-functions appear for instance in the field of exponential integrators, see for instance [27]. For the case of A 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 We briefly recapitulate the usual notation and properties of standard Krylov subspaces, see for instance [47]. For a given matrix A ∈ C n×n , a starting vector v ∈ C n and Krylov dimension 0 < m ≤ n, the Krylov subspace is given by Let V m ∈ C n×m represent the orthonormal basis of K m (A, v) with respect to the Hermitian inner product, constructed by the Arnoldi method and satisfying V * m V m = I m×m . Its first column is given by is upper Hessenberg. We further use the notation h m+1,m = (H m+1 ) m+1,m ∈ R, and v m+1 ∈ C n for the (m + 1)-th column of V m+1 , with V * m v m+1 = 0 and v m+1 2 = 1. The Arnoldi decomposition (in exact arithmetic) can be expressed in matrix form, (2.5)

Remark 1
The numerical range W(A) = {y * Ay/y * y, 0 = y ∈ C n } plays a role in our analysis. Note that W(H m ) ⊆ W(A) (see (A.1)).

Remark 2
The case (H m ) k+1,k = 0 occurs if K k (A, v) is an invariant subspace of A, 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 H m are real and positive, 0 < (H m ) j+1, j for j = 1, . . . , m − 1, and 0 < h m+1,m ∈ R.
For the special case of a Hermitian or skew-Hermitian matrix A the Arnoldi iteration simplifies to a three-term recurrence, the so-called 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 [23]. 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 Gram-Schmidt, modified Gram-Schmidt, the Householder algorithm, the Givens algorithm, or variants of Gram-Schmidt with reorthogonalization (see also [47, Algorithm 6.1-6.3] and others). We refer to [7] and references therein for an overview on the stability properties of these different variants.
In the sequel the notation V m , H m , 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, H m remains to be upper Hessenberg with positive lower subdiagonal entries. Assuming floating point arithmetic we use the notation U m ∈ C n×m for a perturbation of the Arnoldi decomposition (2.5) caused by round-off, i.e., An upper norm bound for U m was first introduced in [43] for the Lanczos iteration in real arithmetic. For different variants of the Arnoldi or Lanczos iteration this is discussed in [58] and others. We assume U m 2 is bounded by a constant C 1 which can depend on m and n in a moderate way and is sufficiently small in a typical setting, We further assume that the normalization of the columns of V m is accurate, in particular that the (m + 1)-th basis vector v m+1 is normalized correctly up round-off with a sufficiently small constant C 2 (see e.g. [43, eq. (14)]), Concerning V m+1 which represents an orthogonal basis in exact arithmetic, numerical loss of orthogonality has been well-studied. Loss of orthogonality can be significant (see for instance [45,7] and others), depending on the starting vector v. Reorthogonalization schemes or orthogonalization via Householder or Givens algorithm can be used to obtain orthogonality of V m+1 on a sufficiently accurate level. The numerical range of H m obtained in floating point arithmetic (see (2.6)) can be characterized as with U C 3 ε (W(A)) being the neighborhood of W(A) in C with a distance C 3 ε. With the assumption that V m+1 is sufficiently close to orthogonal (e.g., semiorthogonal [50]), the constant C 3 in (2.7c) (which also depends on C 1 and problem sizes) can be shown to be moderate-sized. Further details on this aspect are given in Appendix A.
Krylov approximation of ϕ-functions. 3 Let V m ∈ C n×m , H m ∈ C m×m and β ∈ R be the result of the Arnoldi method in floating point arithmetic for K m (A, v) as described above. For a time-step 0 < t ∈ R and p ≥ 0 the vector ϕ p (tA)v can be approximated in the Krylov subspace K m (A, v) by the Krylov propagator The special case p = 0 reads u 0,m (t) = βV m e tH m e 1 . (2.8b) We remark that the small-dimensional problem ϕ p (tH m )e 1 ∈ C m , typically with m n, can be evaluated cheaply by standard methods. In the sequel we denote y p,m (t) = β ϕ p (tH m )e 1 ∈ C m , i.e., u p,m (t) = V m y p,m (t). (2.9) For p = 0 the small dimensional problem y 0,m (t) = β e tH m e 1 solves the differential equation For later use we introduce the notation y p,m (t) = t p y p,m (t), (2.11a) which for p ∈ N and according to (2.4) satisfies the differential equation Remark 3 Although we take rounding effects in the Arnoldi decomposition into account, we do not give a full study of round-off errors at this point. Round-off errors in substeps such as the evaluation of y p,m (t) or the matrix-vector multiplication V m y p,m (t) will be ignored. We refer to [23] for a more general study of these effects.

Remark 4
In the special cases A = B or A = iB for a Hermitian matrix B ∈ C n×n (with A being skew-Hermitian in the latter case) the orthogonalization of the Krylov basis of K m (B, v) simplifies to a three-term recursion, the so-called Lanczos method. In the skew-Hermitian case (A = iB) the Krylov propagator (2.8a) can be evaluated by βV m ϕ p (itH m )e 1 , i.e., we approximate the function λ → ϕ p (itλ ) in the Krylov subspace K m (B, v). 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 e σtB as introduced in [30], with σ = ±i and a Hermitian matrix B for the skew-Hermitian case.
The error of the Krylov propagator. We denote the error of the Krylov propagator given in (2.9) by We are further interested in computable a posteriori estimates for the error norm ζ p,m (t) ≈ l p,m (t) 2 , which in the best case can be proven to be upper bounds on the error norm l p,m (t) 2 ≤ ζ p,m (t).
Norm estimates of the error (2.12) can be used in practice to stop the Krylov iteration after k steps if l p,k (t) 2 satisfies (2.13) below, or to restrict the time-step t 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 [49,30].
A prominent task is to test if the error norm per unit step is bounded by a tolerance tol, ζ p,m (t) ≤ t · tol, which should entail l p,m (t) 2 ≤ t · tol. (2.13) In case of ζ p,m (t) being an upper bound on the error norm, this results in a reliable bound on the error norm (2.13).

An integral representation for the error of the Krylov propagator
We proceed with discussing the error l p,m of the Krlyov propagator. To this end we first define its scalar defect by δ p,m (t) = β e * m t p ϕ p (tH m )e 1 = t p y p,m (t) m ∈ C, (3.1a) and the defect integral by 4 Theorem 1 Let δ p,m (t) ∈ C be the defect defined in (3.1a). For y p,m (t) ∈ C m defined in (2.9) and a numerical perturbation U m ∈ C n×m of the Arnoldi decomposition (see (2.6)), we have: (a) The error l p,m (t) of the Krylov propagator (see (2.12)) enjoys the integral representation For given machine precision ε and constants C 1 , C 2 representing round-off effects (see (2.7a),(2.7b)), and with κ 1 = max s∈[0,t] e sA 2 and κ 2 = max s∈[0,t] e sH m 2 the error norm is bounded by with the defect integral L p,m (t) defined in (3.1b).

Proof
(a) For the exact matrix function we use the notation u p (t) = ϕ p (tA)v, and w p (t) = t p u p (t).
For the Krylov propagator we denote u p,m (t) = V m y p,m (t) with y p,m (t) = β ϕ p (tH m )e 1 (see (2.9)), and we also define w p,m (t) = t p u p,m (t) = V m y p,m (t), with y p,m (t) = t p y p,m (t) defined in (2.11a).
Local error representation in terms of the defect. We defined the scaled error -For p ∈ N the scaled error satisfies with the defect of w p,m (t) with respect to the differential equation (3.3), Together with (2.6) and using of βV m e 1 = v the defect can be written as -For p = 0, in an analogous way we obtain We conclude with the scalar defect defined in (3.1a). Due to (3.4) we have and for l p,m (t) = t −p l p,m (t) together with (3.5) this implies (3.2a).
With the defect integral L p,m (t) defined in (3.1b) we obtain the first term in (3.2b). For the second integral term (with y p,m (t) = β ϕ p (tH m )e 1 ) we use the upper bound -For p ∈ N we apply the integral representation due to (2.3) for ϕ p (tH m )e 1 to obtain the norm bound -For p = 0 we obtain (3.8) in a direct way. Combining (3.7) with (3.8) and denoting κ 2 = max s∈[0,t] e sH m 2 we obtain Combining these estimates with (3.6) we conclude (3.2b).

Remark 5
The error norm of the Krylov propagator scales with κ 1 = max s∈[0,t] e sA 2 and κ 2 = max s∈[0,t] e sH m 2 in a natural way. 5 . It is well known that e tA 2 ≤ e tµ 2 (A) with the logarithmic norm µ 2 (A) = max{Re(W(A))} = max{spec(A + A * )/2}, see for instance [24,Theorem 10.11]. Problems with µ 2 (A) > 0 can be arbitrary ill-conditioned and difficult to solve with proper accuracy. (For further results on the stability of the matrix exponential see also [36,55].) We will not further discuss problems with µ 2 (A) > 0 and assume µ 2 (A) ≤ 0. We refer to the case µ 2 (A) ≤ 0 as the dissipative case, with κ 1 = 1.
For the dissipative case with µ 2 (A) ≤ 0 the error bound (3.2b) from Theorem 1 reads The dissipative behavior of e tA carries over to the Krylov propagator up to a perturbation which depends on round-off errors, including the loss of orthogonality of V m . In terms of the numerical range Our aim is to construct an upper norm bound for the error per unit step (2.13) via (3.9). Let the tolerance tol be given and t be a respective time step for (2.13). Then the round-off error terms in (3.9) are negligible if Concerning the constants C 1 , C 2 and C 3 see (2.7). We recapitulate that C 1 and C 2 given in (2.7a) and (2.7b) can be considered to be small enough in a standard Krylov setting. The constant C 3 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 C 3 only appears as an exponential prefactor for the round-off term in (3.10) and is less critical compared to C 1 and C 2 .
With the previous observation on the round-off 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 tol, see (3.10).
Corollary 1 For the case µ 2 (A) ≤ 0 and with the assumption that round-off error is negligible, the error of the Krylov propagator is bounded by the defect integral L p,m (t), Note that the defect norm |δ p,m (s)| 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 f over the nodes λ 1 , . . . , λ m . (This is to be understood in the confluent sense for the case of multiple nodes λ j , see for instance [24,Section B.16].) Theorem 2 (see for instance [9]) Let H m ∈ C m×m be an upper Hessenberg matrix with positive secondary diagonal entries, 0 < (H m ) j+1, j ∈ R for j = 1, . . . , m − 1, and eigenvalues λ 1 , . . . , λ m . Let f be an analytic function for which f (H m ) is well defined. Then, For f = (ϕ p ) t : λ → ϕ p (tλ ) we will also make use of the following result. 6 Theorem 3 (Corollary 1 in [49]) (Expressing ϕ-functions via dilated exp-functions.) For t ∈ R, The matrix H p,m in Theorem 3 is block triangular with eigenvalues equal to those of H m and J p×p . Therefore, spec( H m ) = {λ 1 , . . . , λ m , 0, . . . , 0}, with 0 as an eigenvalue of multiplicity p (at least). In our context, H m is upper Hessenberg with a positive lower secondary diagonal and ( H m ) j+1, j . In accordance with Theorem 2 the result of Theorem 3 holds for divided differences in a similar manner, With Theorem 2 and 3 the following equivalent formulations can be used the rewrite the scalar defect δ p,m (t) defined in (3.1a).
Let H p,m ∈ C m+p be given as in Theorem 3. For the scalar defect we obtain the following equivalent formulations: We remark that the eigenvalues λ 1 , . . . , λ m of the Krylov Hessenberg matrix H m are also referred to as Ritz values (of A) in the literature.

Computable a posteriori error bounds for the Krylov propagator
The following two propositions are used for the proof of Theorem 4 below. 8 Proposition 1 For arbitrary nodes λ j ∈ C and p ∈ N 0 , Proposition 2 (Lemma including eq. (5.1.1) in [35]) For arbitrary nodes λ j = ξ j + iη j ∈ C, Proof See Appendix B.
We now derive upper bounds for the error via its representation by the defect integral (3.1b). 6 Theorem 3 can be generalized to the case t p e * m ϕ k+p (tH m )e 1 = e * m+p ϕ k (t H p,m )e 1 with k ∈ N, see [2, Theorem 2.1]. The case k = 0 is sufficient for our purpose. 7 Here we introduce the notation (λ 1 , . . . , λ m , 0 p ) = (λ 1 , . . . , λ m , 0, . . . , 0) ∈ C m+p for p ∈ N 0 . 8 We use the notation introduced in the previous sections.
For the case of H m having real eigenvalues, the assertion of Theorem 4 can be reformulated in the following way (see [30,Proposition 6]).
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 but may be less tight. Using the Mean Value Theorem, [24, eq. (B.26)] or [4, eq. (44)], for the divided differences in Theorem 4, eq. (4.1) we obtain the following result which corresponds to [30,Theorem 1 and 2]. For the exponential of a skew-Hermitian matrix a similar error estimate has been used in [33] and is based on ideas of [44] with some lack of theory.
Corollary 4 Let p ∈ N 0 , µ 2 (A) ≤ 0, and assume that round-off errors are sufficiently small (see Corollary 1). Let ξ max = 0 for p ∈ N and ξ max = max j=1,...,m ξ j ≤ 0 for p = 0 and eigenvalues λ j = ξ j + iµ j ∈ C of H m . An upper bound on the error norm is given by For the case of H m having purely imaginary eigenvalues, the divided differences in Theorem 4 (see (4.1)) can be evaluated directly via [24, eq. (B.27)], hence the assertions of Theorem 4 and Corollary 4 coincide in this case.
Accuracy of the previously specified upper bounds on the error norm. In the following we again denote λ 1 , . . . , λ m ∈ C for the eigenvalues of H m , with λ j = ξ j + iη j . For the scalar defect δ p,m (t) (see (3.1a)) we recapitulate Corollary 2, in particular Theorem 4 and its corollaries make use of the error bound given in Corollary 1 and computable upper bounds on the defect integral L p,m (t). A refinement of the upper bound from Corollary 1 would require further applications of the large-dimensional matrix-vector product with A ∈ C n×n and has been shown to be inefficient in terms of computational cost, see also [30,Remark 7]. The computable upper bounds on the defect integral L p,m (t) will be further discussed. We recapitulate the upper bound of the divided differences given in Proposition 2, Thus, in the case of H m having eigenvalues with a sufficiently small imaginary part, the upper bound in Proposition 2, is tight. In the following proposition this statement is made more precise.
Proposition 3 (Part of a proof in [35], eq. (5.2.3)) For nodes λ j = ξ j + iη j ∈ C and t ≥ 0 with max j t|η j | ≤ η t < π/2, Under the assumptions of Proposition 3 we conclude Summarizing, we see that the defect integral can be computed exactly for the case of H m having real eigenvalues (Corollary 3), and a computable upper bound can be given which is tight for the case of H m having eigenvalues sufficiently close to the real axis (Theorem 4 and eq. (4.7)).

Remark 6
For H m with purely imaginary eigenvalues (λ j ∈ iR), e.g. in the skew-Hermitian case, the following asymptotic expansion for the defect is obtained from (4.8), 9 We use the expansion from (4.8) for | exp t [λ 1 , . . . , λ m , 0 p ]| and exp t [ξ 1 , . . . , ξ m , 0 p ] to obtain (4.10) Termwise integration of (4.10) and the proper prefactor gives an asymptotic expansion for the defect integral L p,m (t), similar to (4.7), Omitting further details we state that (4.11) is to be understood in an asymptotic sense with an remainder of O(t 3 |ξ ||η| 2 +t 4 |η| 4 ). In contrast to (4.7) the remainder is depending on ξ terms but (4.11) reveals further constants which can be relevant for practical applications. A similar criterion can be given for the accuracy of the upper bound, which appears in Corollary 4 (with ξ max = 0) and [30, Theorem 1 and 2]. With (4.8), and ρ 1 and ρ 2 given therein, the defect integral can be written as for t → 0. In contrast to the error bound in Corollary 4, the formulas for ρ 1 and ρ 2 in (4.8) require the evaluation of the eigenvalues of H m . The following Proposition gives a formula for ρ 1 and ρ 2 which does not require computation of the eigenvalues of H m and can be evaluated on the fly.
Proposition 4 (Evaluation of ρ 1 and ρ 2 in terms of entries of H m ) The coefficients ρ 1 and ρ 2 in (4.8) can be rewritten as Proof For the coefficients ρ 1 and ρ 2 we use (C.17) with m ← m + p and S 1 and S 2 from (C.3). For the nodes λ 1 , . . . , λ m , 0 p (with λ 1 , . . . , λ m eigenvalues of H m ) we obtain (4.14) The identity for Trace(H 2 m ) in (4.14) holds true due to the upper Hessenberg structure of H m .
Following the proof of Theorem 5 we observe that the case ρ 1 = 0 is possible but results in ρ 2 = 0.  .6)). For the case of H m having eigenvalues with a significant imaginary part, tight estimates are more difficult to obtain. It can be favorable to approximate the defect integral (3.1b) by quadrature to obtain an error estimate via Corollary 1. The aim of using quadrature is to obtain an error estimate which is tighter compared to previous upper norm bounds on the error. In contrast to the proven upper error bounds given in Theorem 4, Corollary 3 and 4 the following quadrature estimates do not result in upper error bounds in general. However, in many practical cases such quadrature estimates turn out to be still reliable.
Here, some remarks on the defect are in order to explain some subtleties with quadrature estimates for the defect integral L p,m (t). We discuss a test problem with a skew-Hermitian matrix A ∈ C n×n . Following Remark 4 we choose A = iB with a Hermitian matrix B, in particulary, B = tridiag(1, −2, 1) ∈ R n×n with n = 10 000. The matrix B is related to a finite difference discretization of the one-dimensional Laplacian operator and A corresponds to a free Schrödinger type problem. The eigenvalues σ j , for j = 1, . . . , n, of B are well studied, we obtain σ j = 4 sin( jπ/(2(n + 1))) 2 with respective eigenvector ψ j ∈ R n . (4.15) Here, µ 2 (A) = 0, and the conditions of Corollary 1 hold. For a given starting vector v ∈ C n the time propagation for the discretized free Schrödinger equation is given by exp(tA)v and can be approximated by the Krylov propagator with p = 0. The following different cases for the starting vector v will be discussed.
In addition to the setting from (a)-(c) we normalize v, v 2 = 1. The defect δ p,m (t) for p = 0 is computed in MATLAB, using expm to evaluate the matrix exponential of H m and divided differences for a fixed Krylov dimension m = 20.
In Fig. 1 we observe |δ p,m (t)| = O(t m−1 ) (for t → 0) up to t ≈ 10 1 for the case (a)-(c). The values of |δ p,m (t)| in this time regime vary strongly between among these cases. We further remark that in the case (b) for t ≥ 4 · 10 1 the defect |δ p,m (t)| behaves similar to the divided differences of the exponential over the first eigenvalues λ ] is illustrated by ('+'). The asymptotic expansion of the divided differences for t → 0 given in (4.9) is illustrated using dashed lines. The dash-dotted line is O(t 6 ).
As a conclusion from the example in Fig.3 1, we observe that quadrature of the defect can be relevant up to a time t for which the quadrature based estimate of l p,m (t) 2 (via the defect integral) is equal to a given tolerance, see (2.13). This regime of t would depend on the choice of tol and additional factors such as β , h m+1,m etc. which appear in the error bound from Corollary 1. Depending on parameters and the starting vector v the defect can be highly oscillatory for relevant times t and, respectively, a quadrature estimate of the defect integral can be difficult to obtain. Such effects seem to be relevant for special choices of starting vectors v, for example case (b) and (c). The effect of H m having clustered eigenvalues and the prefactor used in Fig. 1 ('+') are explained in the following model problem, see Fig. 2.
Divided differences with clustered nodes: an example. Choose m = 3 with nodes a 1 = 1.123, a 2 = 1.231, a 3 = 5.43. With this choice we obtain cluster of nodes, a 1 ≈ a 2 . For the given example we obtain | exp t [ia 2 , ia 3 ]| | exp t [ia 1 , ia 2 ]| for t large enough, hence, using the recursive definition of the divided differences (see [24, eq. (B.24)] or others) we obtain This example is illustrated in Fig. 2. This behavior can be generalized for a larger number of nodes and is also observed in Fig. 1.
Quadrature estimates for the defect integral. With the previous observations on the defect we now discuss different quadrature-based estimates. We can extend the result of the generalized residual estimate, which was introduced in [28] and appeared in a similar manner in [11,46,34,5], to ϕ-functions using the defect integral according to Corollary 1.
Remark 9 (Generalized residual estimate, see also [28]) Applying the right-endpoint rectangle rule we have and with Corollary 1 (and δ p,m (t) given in (3.1a)) we obtain the error estimate Assume that max s∈[0,t] |δ p,m (t)| = |δ p,m (t)|, e.g. |δ p,m (t)| is monotonically increasing in t. Then,  In this case the generalized residual estimate from Remark 9 results in an upper bound on the error norm.
In the most general case the defect is of a high order for t → 0 and in a relevant time regime, see also Fig. 1 case (a) and previous remarks. Then the defect is a higher order function and the right-endpoint quadrature does result in an upper bound but is not tight. In this case we can improve the estimate by a prefactor depending on the effective order defined in Appendix C. If the defect is sufficiently smooth in a relevant time regime this results in a tight upper bound on the error norm.
Remark 10 (Effective order estimate, see also [30]) Denote f (t) = | exp t [λ 1 , . . . , λ m , 0 p ]| for the timedependent part of the defect with eigenvalues λ 1 , . . . , λ m of H m . Assume f (t) > 0 for a sufficiently small time regime t > 0. We consider the effective order ρ(t) to be defined for the divided differences f (t) as given in (C.4a). With the following estimate for the integral of the defect, and from Corollary 1 (with δ p,m (t) given in (3.1a)) we obtain In [30] the effective order is defined for |e * m e tH m e 1 | (p = 0) which is equivalent to the definition via the divided differences of f (t). (This follows from Corollary 2 and the definition of the effective order which is independent of a constant prefactor.) Some of the following observations already appeared in [30]. The quadrature scheme in Remark 10 is motivated by the following relation of the effective order and the integral of the divided differences f (t). From eq. (C.4a), Integration and application of the mean value theorem shows the existence of and integration by parts gives .
Combining (4.17) and Corollary 1 we obtain the upper bound A computable expression for the effective order was given in [30, eq. (6.10)]. This result can be generalized to the case p ∈ N 0 , ρ(t) = t Re (H m ) m,m + (H m ) m,m−1 (y p,m (t)) m−1 /(y p,m (t)) m for p = 0, and Re((y p−1,m (t)) m /(y p,m (t)) m ) for p ∈ N, with y p,m (t) ∈ C m from (2.9). The expression for the case p ∈ N can be obtained by [30, eq. (6.10)] applied on the representation |e * m+p e t H m e 1 | for the defect (iii. in Corollary 2) and making use of the special structure of H m , β e * m+p e t H m e 1 = t p (y p,m (t)) m (see Corollary 2) and β e * m+p−1 e t H m e 1 = t p−1 (y p−1,m (t)) m (see [49,Corollary 1]).
As illustrated in Fig. 1 the defect can be highly oscillatory in a relevant time regime, especially for specific starting vectors, and in this case the quadrature estimates should be handled with care.

A stopping criterion for lucky breakdown.
The special case h k+1,k = 0 during the construction of the Krylov subspace is considered to be a lucky breakdown, a breakdown of the Arnoldi or Lanczos iteration with the benefit of an exact approximation of ϕ p (tA)v for any t > 0 via the Krylov subspace K k (A, v). In floating point arithmetic the lucky breakdown results in h k+1,k ≈ 0 and can lead to stability issues if the Arnoldi or Lanczos method is not stopped properly. The condition that the Krylov propagator is exact is not exactly determinable in floating point arithmetic but can be weakened to the error condition in (2.13) for a given tolerance tol per unit step. With this approach we introduce a stopping criterion which can be applied on the fly to detect a lucky breakdown and satisfies an error bound. This does not depend on any a priori information as long the tolerance tol is chosen properly so that round-off errors can be neglected, see remarks before Corollary 1.
Proposition 5 Let µ 2 (A) ≤ 0 and assume that round-off errors are sufficiently small, see Corollary 1. Let tol be a given tolerance and β h k+1,k (p + 1)! ≤ tol (4.18) be satisfied at the k-th step of the Arnoldi or Lanczos iteration. Then the iteration can be stopped and the Krylov subspace K k (A, v) can be used to approximate the vector ϕ p (tA)v with a respective error per unit step l p,k (t) 2 ≤ t · tol.
Proof We use the upper bound on the error norm from Corollary 1, To obtain a uniform bound on the defect integral we use Together with (4.19) and (4.18) we conclude l p,k (t) 2 ≤ t · tol.

Numerical experiments
The notation for the error l p,m (t), the estimate of the error norm ζ p,m (t) and the tolerance tol have been introduced in (2.12) and (2.13). The notation ζ p,m will be used for different choices of error estimates discussed in the previous section. Theorem 4 and Corollary 4 result in upper bounds on the error norm, l p,m (t) 2 ≤ ζ p,m (t) . The quadrature-based error estimates given in Remark 9 and 10 result in estimates for the error norm, l p,m (t) 2 ≈ ζ p,m (t), and with additional conditions also give upper bounds. For a fixed tolerance tol we use the notation t(m) for the smallest time t with ζ p,m (t) = t · tol, see (2.13). This choice of t(m) helps us to verify the tested error estimates for a time t which is of the most practical interest. With the help of a reference solution the true error norm per unit step can be tested by l p,m (t(m)) 2 /t(m). For the numerical experiments we focus on the most prominent case p = 0 and simplify the notation by writing δ m (t) = δ 0,m (t), l m (t) = l 0,m (t) and ζ m (t) = ζ 0,m (t).

Convection-diffusion equation
Consider the following two-dimensional convection-diffusion equation with t ≥ 0 and x ∈ [0, 1] 2 , Let A ∈ R n×n be obtained by the two-dimensional finite difference discretization of the operator L in (5.1) with zero dirichlet boundary conditions and N = 500 inner mesh points in each spatial direction, hence, n = N 2 . This test problem is similar to other convection-diffusion equations appearing in the study of Krylov subspace methods, see also [30,15,19,6] and others. The symmetric case with ν = 0 results in the Heat equation and has already been discussed in [30]. For the convection parameter we choose ν = 100, 500 which results in a non-normal matrix A. Considering the spectrum of A the case ν = 100 is closer to the Hermitian case and ν = 500 is closer to the skew-Hermitian case. We remark that for the real matrix A ∈ R n×n the terms Hermitian and symmetric can be used in an equivalent manner. In both cases the numerical range of A is in the left complex plane, µ 2 (A) ≤ 0.
We discuss error estimates for the case p = 0, hence, e tA v is approximated in the Krylov subspace K m (A, v), see (2.8b). As a starting vector we choose the normalized vector v = (1, . . . , 1) * ∈ R n . The error estimates given in Theorem 4, Corollary 4 and Remark 9 and 10 are compared for this problem in Fig. 3.
For the case ν = 100 the eigenvalues of H m have a negligible imaginary part and the upper bound given in Theorem 4 shows tight results. For ν = 500 and larger choices of m this bound is less tight. The criterion ac.est.1(t) given in Remark 7 is evaluated for ν = 100, 500 with t(m) corresponding to Theorem 4 (see caption of Fig. 3). For ν = 100 we obtain ac.est.1(t(m)) < 0.1 for any m tested and for ν = 500 the smallest m with ac.est.1(t(m)) > 0.1 is m = 45. The upper bound of Corollary 4 is applied with ξ max = 0 (the effect of ξ max is negligible in this case). Similar to the criterion ac.est.1(t), we test ac.est.2(t) given in Remark 8 for t(m) corresponding Corollary 4. The smallest m with ac.est.2(t(m)) > 0.1 is m = 9 and m = 10 for ν = 100 and ν = 500, respectively. The quadrature-based error estimates given in Remark 9 and 10 and both result in upper bounds on the error norm for this example, whereas the effective order estimate in Remark 10 results in a tighter bound.

Free Schrödinger equation with a double well potential
Consider the one-dimensional free Schrödinger equation with a double well potential. which is evaluated on the mesh and normalized to obtain a discrete starting vector v ∈ R n . This problem also appears in [29,51]. Similar to the previous subsection we discuss error estimates for the case p = 0, hence, the Krylov approximation of e −itB v. The implementation of the skew-Hermitian problem is described in Remark 4. In Fig. 4 the upper bound given in Theorem 4 and Corollary 4, which coincidence in the skew-Hermitian case, and the error estimates given in Remark 9 and 10 are compared. With our choice starting vector the matrix H m does have clustered eigenvalues, see also Subsection 4.1. The defect δ m (t), which is presented in the lower right corner of Fig. 4, does have an oscillatory behavior which is luckily not in the relevant time regime. Therefore, with our choice of tolerance tol = 10 −8 the quadrature-based error estimates are still valid. In other cases this oscillatory behavior of the defect can lead to failure of the error estimates given in Remark 9 and 10. The upper bound given in Corollary 4 is reliable but not tight for this example which can also be explained by the loss of order of the defect caused by the starting vector.

Conclusions and outlook
In this work various a posteriori bounds and estimates on the error norm, which have their origin in an integral representation of the error using the defect (residual), are studied. We have characterized the accuracy of these error bounds by the positioning of Ritz values (i.e., eigenvalues of H m ) on the complex plane. The case of real Ritz values is the most favorable one to obtain a tight error bound via an integral on the defect norm (Corollary 3). A new error bound (Theorem 4) has shown to be tight if Ritz values are close to the real axis and in this case favorably compares with existing error bounds. We further recapitulate an existing error bound (Corollary 4) which remains relevant, especially for the case of Ritz values with a significant imaginary part. In addition for the error bound in Theorem 4 and Corollary 4, we have provided a criterion to quantify the achieved accuracy on the run. For an illustration of the claims concerning the new error bound we primary refer to the numerical example given in Subsection 5.1. The quadrature-based error estimates in Subsection 4.1 (e.g. the generalized residual estimate) do not yield proven upper bounds on the error norm and we adressed special cases (e.g. the numerical example in Subsection 5.2) for which the reliability of these estimates can be problematic. These cases are also analyzed in terms of Ritz values in Subsection 4.1 and this relation can be of further interested for a numerical implementation. Nevertheless, in most cases the quadrature-based estimates remain valid, whereat the effective order quadrature stands out in terms of performance.
We also remark that the theory provided in our work gives the possibility to adapt the choice of the error estimate on the fly to obtain an estimate which is as reliable, accurate and economic as possible. A numerical implementation is the topic of further work.
A Properties of the Krylov subspace in exact and floating point arithmetic Similar results hold in floating point arithmetic with relative machine precision ε and certain additional assumptions. Assume there exists an orthonormal basis V m ∈ C n×m and a perturbation U m ∈ C n×m , which is sufficiently small in norm (i.e., there exists a moderate constant C 3 with U m 2 ≤ C 3 ε), with With assumption (A.2) and basic properties of the numerical range we obtain Then we combine (A.3) and (A.4) and make use of U m 2 ≤ C 3 ε to obtain In [50,Theorem 5] the existence of the representation (A.2) is proven for the Lanczos method with a sufficiently small constant C 3 and the assumption that the Krylov basis is semiorthogonal.
For the general case of the Arnoldi method the representation (A.2) can be derived using (2.6), (2.7a) and an additional condition on the level of orthogonality of the Krylov basis, e.g., assuming that an orthonormal basis V m exists for which V m − V m 2 is small enough (see also [7, Theorem 2.1] and references therein).

B Some properties of divided differences
Proof (of Proposition 1) For p ∈ N 0 and any A ∈ C m×m , w ∈ C m , from the series representation (2.2) we obtain This identity carries over to divided differences in the following way. Let As a consequence of the Opitz formula, see [42] and remarks in [4, Proposition 25], we have which completes the proof.
Here, in the last step we have used the Mean Value Theorem for the integral. In this way we end up with the estimate With |tx|, |ty| ≤ η t < π/2 we obtain cos( η t ) ≤ cos(tx) ≤ | cos(tx) + i sin(ty)|, which completes the proof.
With the notation (C.1) and the series representation of the exponential function we obtain We also introduce the notation For κ 0 , κ 1 and κ 2 we obtain the following formula.
The following asymptotic expansion of f (t) for t → 0 is motivated by the concept of effective order. The effective order of the function f (t) can be understood as the slope of the double-logarithmic function ln( f (e τ )) with τ = lnt, and with derivative f (e τ ) e τ f (e τ ) .
We denote the effective order by satisfying ρ(t)/t = log( f (t)) . (C.4b) We now analyze the divided differences close to an asymptotic regime under the assumption f (t) > 0, which holds for sufficiently small t > 0. The effective order ρ(t) is then well-defined by (C.4a). The following expansion (C.5) for ρ(t) is be considered in an asymptotic sense for t → 0; convergence of the series is not an issue here.

(C.17)
To study the influence of the real and imaginary parts of the nodes λ j = ξ j + iη j we introduce the notation After all these technicalities we arrive at the following asymptotic expansion.