1. Introduction. Consider a matrix A ∈ R n×n and a function f that is analytic in a neighborhood of the numerical range of A. This paper studies efficient iterative methods for approximating a matrix function f (A) multiplied by a vector b ∈ R n . For large-scale problems, we approximate f (A)b by restricting A to a subspace of dimension m (m n) and obtain

where V m ∈ R n×m contains orthonormal basis vectors of the subspace and A m = V * m AV m is the restriction of A to this subspace. Matrix function problems arise in the numerical solution of differential equations [14,37,47], matrix functional integrators [44,45], model order reduction [2,31], optimization problems [9,64], and others.

One of the classical methods of subspace projection for matrix function approximations is the standard Krylov subspace method that generates the subspaces K m (A, b) = span b, Ab, . . . , A m-1 b ; see, e.g., [42,55]. A few restarted variants were proposed in [1,26,27,32,33]. Methods based on rational approximations have also been studied, such as the extended Krylov subspace method (EKSM) [22,46] and the adaptive rational Krylov subspace method (RKSM) [23,24,38,39,50]. In this paper, we consider a generic RKSM that generates subspaces of the form

where q m-1 (A) is a polynomial of degree not larger than m -1 with respect to A.

However, it is impossible to compute the residual norm directly because f (A)b is unknown. A practical stopping criterion is to monitor the difference between the approximations obtained in two successive iterations. The RKSM can be terminated if such a difference becomes sufficiently small, but this criterion may lead to premature termination if the approximation stagnates without actual convergence to f (A)b. A more reliable alternative stopping criterion is to evaluate upper bounds for the residual norm [21,43,63], especially for exponential-type functions. There are also some results on a posteriori error bounds [35]. In addition, we may compute |e * m f (A m )e 1 | to evaluate the accuracy of the approximation; see, e.g., [11,21,46,59]. In this paper, our stopping criterion is based on the norm of (AV m -V m A m )f (A m )V * m b, which can be interpreted as the residual norm of an associated differential equation; see, e.g., [11,21,56].

Given a square band matrix B and a sufficiently regular function f , the magnitude of the entries of f (B) below the main diagonal can be characterized by a decaying behavior that depends on the row index relative to the diagonal [4,6,7,52]. A priori estimates of the decay rate have been discussed in [5,15,18,34]. For RKSMs, the restricted matrix A m is not banded, but upper bounds for the entries of f (A m ) have been derived, which also exhibit a decaying behavior below the main diagonal [57]. In this paper, we further show a similar decaying behavior for the entries of K -1 m f (A m ) below the diagonal, where K m is an upper Hessenberg matrix generated by the RKSM. The matrix K -1 m f (A m ) is directly related to the residual of the RKSM, and it can be used to determine an a priori tolerance relaxation for the linear solve at each RKSM step to enable an inexact RKSM for approximating f (A)b.

Specifically, at each step of the RKSM, we compute a shift-invert matrix-vector product of the form (A -sI) -1 (A -σI)u, which is equivalent to the solution of the linear system (A -sI)x = (A -σI)u. For large-scale problems, these linear systems are solved approximately by iterative methods. Earlier studies on the inexact Krylov methods based on inexact matrix vector products can be found in [13,61]. Errors are introduced at each RKSM step, and they accumulate in the rational Krylov subspace. The motivation of this paper is to find a strategy to relax the accuracy of the linear solves without negatively impacting the convergence of the RKSM to f (A)b. This motivation is the same as that for the study of inexact standard [20,56] and rational [8,40] Krylov methods for approximating f (A)b for Hermitian matrices. In particular, in [40] a strategy of relaxing the inner tolerance of the shift-invert Lanczos method or the EKSM is proposed that applies only one fixed pole repeatedly; in [8] an effective preconditioner construction for the iterative solution of linear systems with different shifts in the RKSM are considered, but a relaxation of the tolerance of the inner linear solves is not discussed. Inexact RKSMs has also been used in evolution equations [41], Lyapunov equations [48], eigenvalue problems [49,66], and model reductions [65]. In this study, we consider RKSMs for matrices regardless of their symmetry and focus on how the tolerances of the inner linear systems with different shifts can be relaxed without delaying the convergence of the RKSM.

The inexact RKSM in our problem setting relies on the association between the upper bounds for the residual norm for approximating f (A)b and the decay bounds for the entries below the main diagonal of K -1 m f (A m ). Such associations can be established for both analytic functions and Markov functions. These results are largely consistent with the relationships between the poles and the convergence of the RKSM for approximating the actions of the exponential function [53] and Markov functions [3]. A tolerance relaxation strategy for the iterative linear solve at each step of the inexact RKSM is derived based on the decay bounds for the entries of K -1 m f (A m ). Compared with the decay bounds for the entries of f (A m ) in [57], our bounds for the entries of K -1 m f (A m ) are sharper as they keep the original integrals, which can be evaluated efficiently by appropriate quadrature rules; more importantly, they directly inform and enable an inexact RKSM in our problem setting.

The rest of the paper is organized as follows. In Section 2, we review the RKSM for approximating f (A)b and derive a sparsity pattern of certain rational functions of matrices restricted to the RKSM subspaces. In Section 3, we introduce the theorems for the decay bounds for the entries of K -1 m f (A m ) below the diagonal. A tolerance relaxation strategy is derived for the inexact RKSM in Section 4 to ensure that the difference between the true and the derived residuals of the inexact method is bounded by a given tolerance. A heuristic tolerance relaxation strategy is proposed for the inexact linear solves arising at each RKSM step. In Section 5, we show numerical results to support the theorems of the decay bounds for the RKSM residual norms and also show the advantage of the inexact method with a heuristic tolerance relaxation strategy over the exact method. In Section 6, our main theorems are proved using the Faber-Dzhrbashyan rational approximations. Conclusions of this paper are given in Section 7.

In this section, we give some preliminary results to facilitate our later discussion of the inexact RKSM for approximating f (A)b.

A brief review of the RKSM. The RKSM starts with a vector b ∈ R n \ {0} to construct rational Krylov subspaces Q m (A, v 1 ), where A ∈ R n×n and v 1 = b/ b 2 . At step k, the RKSM chooses a pole s k = 0 and a zero σ k = s k , and applies the linear operator (I -A/s k ) -1 (A -σ k I) to the vector v k , which is the last vector of the orthonormal basis vectors {v 1 , v 2 , . . . , v k } that span the current subspace Q k (A, v 1 ). To build an orthonormal basis of the enlarged subspace Q k+1 (A, v 1 ), we adopt the modified Gram-Schmidt orthogonalization and obtain (2.1)

Repeat the above operation for each index value k = 1, 2, . . . , m, assuming that there is no breakdown. It is not difficult to get the rational Arnoldi relation:

or equivalently, (2.2)

where

] contains the orthonormal basis vectors of the rational Krylov subspace

/s m ) and P m = diag(σ 1 , . . . , σ m ), and K m , G m ∈ R (m+1)×m are both upper Hessenberg matrices:

(2.

4) ETNA Kent State University and Johann Radon Institute (RICAM) INEXACT RKSM FOR APPROXIMATING THE ACTION OF FUNCTIONS OF MATRICES 541 From (2.2), it follows that (2.5)

An alternative approach to enlarge the rational Krylov subspace is to apply the operator (I -A/s k ) -1 to the vector v k ; see, e.g., [24,57]. It generates exactly the same subspace as in (2.3) if all poles s k (1 ≤ k ≤ m) remain the same for both approaches. In this paper, we follow (2.1) to enlarge the subspace because our numerical experience suggests that applying the operator (I -A/s k ) -1 (A -σ k I) tends to achieve a smaller final residual norm. This can probably be attributed to an improvement in floating-point accuracy, though we have no additional insight here. The choice of σ k = s k does not impact the generated rational Krylov subspace, and we follow [38] to set σ k = σ = 1 for all 1 ≤ k < m.

and f is analytic in a neighborhood of the numerical range of A. Since the initial basis vector is v 1 = b/β, where β = ||b|| 2 , the RKSM approximation at step k is defined as

where

For certain functions f , we may instead define an alternative residual, associated with an ordinary differential equation that depends on f ; see, e.g., [11,56]. For the exact solution y = f (A)b, this alternative residual is zero. By the argument of continuity, if an approximate solution y m is sufficiently close to y, the alternative residual should be sufficiently small in norm. We give one example to demonstrate this point. EXAMPLE 2.1. Assume that the matrix A has no negative real eigenvalues. Consider the elliptic Dirichlet problem:

The exact solution is y(t) = exp(-t √ A)b. If we denote the RKSM approximation as

Since the residual of the exact solution is R(t) = Ay(t) -y (t) = 0 for all t ≥ 0, the residual R m (t) = Ay m (t) -y m (t) should have a small norm if y m (t) ≈ y(t). In particular, at t = 1, the residual of y m (t) is

where

There are more examples based on differential equations showing that the residual R m = (AV m -V m A m )f (A m )βe 1 can be used to determine the accuracy of the RKSM approximation y m ≈ f (A)b for other functions; see, e.g., [10,11,12,21].

From the Arnoldi relation in (2.2) and (2.5), we get

Following the implementation by Güttel in [38], at step k (k ≤ m), we can first temporarily choose the infinite pole s k = +∞, so that the Arnoldi relation in (2.2) becomes

We can easily obtain the restricted matrix

with the temporary G k and K k and compute the residual of

not sufficiently small, then we choose a finalized pole s k = 0 and form the finalized G k , K k by updating the last column of the temporary G k and K k and then proceed to the next RKSM step. The description of this method is shown in Algorithm 1.

Input:

. . , v k , and normalize into (a temporary) v k+1 .

Compute the restricted matrix

Compute the approximate solution

Return y k as the approximation to f (A)b.

end if

Determine the finalized pole s k = σ k .

10:

Recompute

. . , v k , and normalize into (the finalized) v k+1 .

11:

Update the last columns of G k and K k . 12: end for Based on the Arnoldi relation in (2.8), for k = m and s m = ∞, we get

Note that, since we usually have h m+1,m = O(1), the residual norm is directly associated with the (m, 1)-entry of the matrix K -1 m f (A m ). 2.3. A sparsity pattern of functions of restricted matrices for the RKSM. In this section, we show that the entries of certain rational functions of restricted matrices constructed by the RKSM have a sparsity pattern. This observation will be used to prove two main theorems in Section 3 on the decay bounds for the entries of K -1 m f (A m ) below the diagonal

INEXACT RKSM FOR APPROXIMATING THE ACTION OF FUNCTIONS OF MATRICES 543 for analytic functions and Markov functions. To this end, we first propose a lemma that states several properties of A m = V * m AV m , the restriction of A to the rational Krylov subspace (2.3). LEMMA 2.2. Let V m ∈ R n×m contain an orthonormal basis of Q m

Define q = q m-1 (A) -1 v 1 , where q m-1 (z) = m-1 j=1 (1 -z/s j ). Define P m as the set of all polynomials of degree less than or equal to m. The following statements hold:

(i) For any matrix X ∈ R m×m and 0 ≤ j ≤ m, which is sufficient to complete the proof. LEMMA 2.3. Suppose that m -1 steps of the RKSM are performed without breakdown as in (2.1), with s i = 0, s i = σ i , for 1 ≤ i < m, which leads to the Arnoldi relation in (2.2). Assume that K m is nonsingular. Define two new vector spaces

The quality of a candidate approximation to f (A)b from the RK subspace

follows the quality of a corresponding rational function

is used in most RKSM implementations, where the continuation vector is chosen as the last vector of the basis that has been generated.

where t ≥ 1 and p j-1 (z) ∈ P j-1 . For any indices k, ℓ

Proof. The ℓ-th orthonormal basis vector v ℓ of Q m (A, v 1 ) can be written as r ℓ-1 (A)v 1 , where r ℓ-1 (z) ∈ P ℓ-1 /q ℓ-1 , such that r

m in (ii), and hence

Left multiplying V * m on both sides of (ii) and letting X = I, we obtain the relation

Note from (ii) that this equality also holds if r m is replaced with r ℓ-1 because r ℓ-1 ∈ P m /q m-1 . Therefore,

Combining (2.13) and (2.14), we get

Left multiplying v * k on both sides, we have

By Lemma 2.3 and the fact that r (t)

By (2.1), we have

.17) ETNA Kent State University and Johann Radon Institute (RICAM) INEXACT RKSM FOR APPROXIMATING THE ACTION OF FUNCTIONS OF MATRICES 545 Left multiplying K -1 m V * m on both sides of (2.17), we get (2.18)

Combining (2.16) and (2.18), we have

for all k ≥ j + t. Note that since ℓ ≤ t and t ≤ k -j, we have ℓ ≤ k -j with j ≥ 1 and hence ℓ < k. This means that the above sparsity pattern holds in the strictly lower triangular portion of K -1 m r (t) j (A m ). Lemma 2.4 shows that for the RKSM, there exists a sparsity pattern for the entries of

involving the class of rational functions (2.12) and the restricted matrices A m obtained by the RKSM, and Figure 2.1 illustrates two examples of the sparsity patterns for certain t-and j-values. In [57], a corresponding result has been derived for the entries of

m is directly associated with the residual of the RKSM approximations for f (A)b as given in (2.9). We will show that this sparsity property helps to establish the two main theorems in Section 3, derive the convergence of the RKSM, and develop the tolerance relaxation for the iterative linear solve at each step of the inexact RKSM.

In this section, we investigate the upper bounds for the entries in the first column of

The core point is to show how quickly these entries decay with the row index k. In Section 4, we shall explore the connections between the decay bounds and the residual of the RKSM for approximating f (A)b.

k K -1 m f (A m )βe 1 for analytic functions and Markov functions. Several estimates for decay bounds for the entries of functions of matrices have been proposed; see, e.g., [4,Theorem 10], [6, Theorem 3.7], [52,Theorem 2.6], and [56,Theorem 2.3]. In this paper, we use the Faber-Dzhrbashyan (FD) rational functions [62, Ch. XIII, Section 3] and [51] to find upper bounds for e * k K -1 m f (A m )βe 1 . We begin with some definitions. For any matrix A ∈ R n×n , we let The matrix function f (A) can be defined by Cauchy's integral formula as follows [42]

where Γ E is a closed contour in E that encloses the spectrum of A.

In [57, Theorem 4.2], decay bounds for the entries of f (A m ) are derived. Since K m in (2.4) is an upper Hessenberg matrix, K -1 m is the inverse function of a band matrix, and decay bounds for the entries of K -1 m can be derived [19]. One can combine the results of the decay bounds for both f (A m ) and K -1 m to get those for K -1 m f (A m ). However, the decay bounds for K -1 m require spectral information of K m defined in (2.4), which has no straightforward connections to W (A). In this paper, we follow the work in [57] to directly derive decay bounds for K -1 m f (A m ) by using the Faber-Dzhrbashyan (FD) rational functions. Our first main theorem reads as follows.

THEOREM 3.2. Assume that K m defined in (2.4) is nonsingular and

Suppose that all poles s 1 , . . . , s m of the RKSM are located in the exterior of the set E, where

where

and c j is independent of the value of τ . Moreover, a simplified bound holds in the form

The proof of Theorem 3.2 is given in Section 6.

Next, we study the bounds for the entries of K -1 m f (A m ) for an important class of nonanalytic functions, namely, Markov (Cauchy-Stieltjes) functions, defined as [3] (3.4)

where µ is a positive measure with supp(µ) ⊂ (-∞, 0]. Markov functions are not analytic in the entire set E τ for any τ > 1 that is not sufficiently small.

Below are two examples of Markov functions. ETNA Kent State University and Johann Radon Institute (RICAM) INEXACT RKSM FOR APPROXIMATING THE ACTION OF FUNCTIONS OF MATRICES 547 EXAMPLE 3.3. For f

The two Markov functions above are not analytic in (-∞, 0]. If we use Theorem 3.2 to determine an upper bound, we should choose τ such that 1 < τ < |φ(0)|. Such a bound is usually a significant overestimate of the actual rate of decay. Instead, we present a similar theorem for Markov functions with decay bounds that are much sharper. Note that [7,34] have established decay bounds for Markov functions of matrices with banded or Kronecker structure, whereas our results hold without assumptions on the matrix structure. THEOREM 3.5. With the same setting as Theorem 3.2, except that f is a Markov function defined in (3.4), and assuming that E lies strictly in the right half complex plane, it holds that

where

Moreover, a simplified bound holds in the form

The proof of Theorem 3.5 is also given in Section 6. [57], we may further relax the bounds for e * k K -1 m f (A m )e ℓ in Theorem 3.2 and Theorem 3.5 by replacing the integrals in the formulas with rough bounds in terms of elementary functions. However, we point out that keeping the integrals and efficiently approximating them by customized quadrature rules can achieve significantly sharper bounds at a cost not much higher than that needed to evaluate the rough bounds without integrals. In fact, based on the numerical tests in [57], theoretical bounds roughly give overestimates that are four orders of magnitude larger than the actual values. We shall present customized numerical quadrature rules to efficiently approximate our bounds with integrals. This will be discussed in detail in Section 5.1.

Poles and the rate of convergence of the RKSM. Though the main goal of this paper is to study the mechanism to enable the inexact RKSM for approximating f (A)b, in this section we give a brief discussion about the implications of the bounds (3.1) and (3.6) to explore numerical or theoretical relationships between the poles and the asymptotic convergence factor of the RKSM in this problem setting. These relationships seem consistent with those established in [53] and [3] for the exponential function and Markov functions, respectively. For the exponential function f (z) = e -hz , which is analytical in the entire complex plane, the Restricted Denominator (RD) rational approximation is a competitive method; see, e.g., [53]. RD rational approximations can be regarded as a special variant of the RKSM with a fixed repeated pole. It is shown in [53] that if W (A) is a sector in the right half complex plane with vertex at the origin and is symmetric with respect to the real axis, then the optimal pole of RD is s 0 = -m/h, where m is the maximum number of RD steps. From the definition of the residual norm in (2.9) and the upper bound (3.3) for analytical functions we get

Assume that there exists a uniform upper bound for both |h m+1,m | and e * m K -1 m independent of m. Then for the RKSM with a fixed repeated pole s ∈ R,

(1 ≤ i ≤ m -1), and therefore,

dθ.

The optimal single repeated pole s = s * ≤ 0 is defined as

We use a composite trapezoidal rule to approximate the integral and use MATLAB's fminbnd (a function which aims to find the minimum of a continuous single-variable function on a finite interval) to approximate the optimal single repeated pole s * . Assume that A is a real nonsymmetric matrix such that W (A) can be covered by an ellipse with semi-major axis of length a parallel to the real axis and semi-minor axis of length b parallel to the imaginary axis (a > b ≥ 0), centered at c ∈ C. The conformal map φ(z) is then defined by

and its inverse is defined by

where ρ = √ a 2 -b 2 and κ = (a + b)/ρ; see, e.g., [57]. For a = b, W (A) can be covered by a circle, so that

For a < b, W (A) can be covered by an ellipse with semi-major axis parallel to the imaginary axis, and we can derive the similar expressions for both φ(z) and ψ(w).

For example, assume that a matrix A is such that its numerical range W (A) can be covered by an ellipse centered at c = 101, with semi-major axis of length a = 100 lying on the 3.1 shows the comparison of the optimal single pole s * computed numerically in (3.8) and s 0 = -m h in [53] for approximating e -hA b. One can see from the table that the two poles are relatively close. The difference between these two poles might be attributed to the different shape of W (A), which is assumed to be an infinite sector in [53] but is a finite ellipse in our experiments. The result in [53] is more suitable for some problems arising from discretizing PDEs with different mesh sizes, because all of them can be fitted into identical sector with infinite radius. In this paper, following the assumptions in the literature on the convergence of the RKSM based on Riemann mappings φ(z) and the inverses ψ(w), we focus on matrices with a finite numerical range.

TABLE 3.1 Comparison of the optimal single pole s * (3.8) and s 0 = -m h [53] for approximating e -hA b, where W (A) can be covered by an ellipse centered at c = 101, with semi-major axis of length a = 100 lying on the real axis and semi-minor axis of length b = 10.

s 0 20 -174.35 -200 0.8717 -16.60 -20 0.8298 -0.92 -2 0.4578 40 -298.70 -400 0.7468 -38.08 -40 0.9521 -2.70 -4 0.6749 60 -437.64 -600 0.7294 -64.02 -60 1.0669 -4.45 -6 0.7418 80 -585.49 -800 0.7319 -93.91 -80 1.1739 -6.20 -8 0.7755 100 -738.70 -1000 0.7387 -126.80 -100 1.2680 -7.97 -10 0.7975 For Markov functions, from the definition of the residual norm (2.9) and (3.7) in Theorem 3.5, we get

Assume that there exists a uniform upper bound for

which is consistent with the findings in [3]. In particular, for a fixed repeated pole s i = s < 0, there exists

, where the last equality holds since w-φ(s) φ(s)w-1 is monotonic in w on (-∞, φ(0)]. Essentially the same result for the residual defined by f (A)b -V m f (A m )βe 1 can be found in [3, Corollary 6.4 (a)]. A similar residual bound obtained with two cyclic poles can also be derived, corresponding to the result in [3, Corollary 6.4 (b)]. The above discussion gives an alternative proof of the residual bounds for approximating the action of a Markov function f (A)b by the RKSM with a few cyclic poles.

4. An inexact RKSM for approximating f (A)b. Inexact Arnoldi algorithms have been widely used in solving numerical linear algebra problems, including approximating f (A)b. In general, these algorithms include inexact standard (polynomial) Krylov methods and inexact rational Krylov methods. They can be applied to symmetric and nonsymmetric matrices, while f can be an analytic function or a Markov function. Preliminary test results were given in [8] for the inexact standard Krylov method to approximate f (A)b, where A is symmetric and positive definite and f is analytic. Inexact standard Krylov subspace methods have also been studied in [20] and [56] for approximating f (A)b, where A is nonsymmetric and f is analytic. Several inexact rational Krylov methods, including the shift-and-invert Lanczos method and the EKSM (which rely on only one fixed pole), have been investigated in [40] for approximating the action of Markov functions of Hermitian matrices to vectors. To the best of our knowledge, no studies have been carried out to explore the inexact rational Krylov method with variable poles for approximating f (A)b involving nonsymmetric matrices for either analytic functions or Markov functions. Our goal of study is to fill this research gap.

For large-scale problems, the approximate computation of the shift-invert matrix vector product w k+1 = (I -A/s k ) -1 (A -σI)v k in (2.1) at step k of the RKSM is done by an iterative linear solver. Errors are introduced in the approximate solution and hence into the basis vectors of the rational Krylov subspaces. Let w k+1 be an approximate solution such that the residual of this linear solve is

If we choose the pole s m = ∞, then the inexact Arnoldi relation of RKSM after step m is

We are interested in exploring strategies to make the difference between the derived residual in (2.9) and the true residual in (4.2) sufficiently small, so that the error term Ξ m has little impact on the convergence of the inexact RKSM. The difference between the two residuals is

From the definition of ∆ m , we have

where the second equality holds since I -V m V * m is an orthogonal projector. In order to make ∆ m sufficiently small, either

Suppose that c m is an upper bound for K -1 m . For any vector w with unit 2-norm, we have

with semi-major axis of length a parallel to the real axis and semi-minor axis of length b (a ≥ b ≥ 0), where c 1 > a. For any point on the boundary of the ellipse that covers W (A), denoted as p = (c 1 + a cos θ, c 2 + b sin θ), consider a corresponding point p * = (c 1 + a cos θ + ε cos α, c 2 + b sin θ + ε sin α). Since the sum of the distances from p to the two foci of the ellipse is 2a, it is easy to show that the sum of the distances from p * to the two foci is less than or equal to 2a + 2ε by applying the triangle inequality involving the three triangles with vertices p, p * , and the two foci. Therefore,

Suppose that for every 1 ≤ k ≤ m, ξ k ≤ k , where

and m i=1 2 i ≤ ε 2 /c 2 m . ETNA Kent State University and Johann Radon Institute (RICAM) 552 S. XU AND F. XUE Proof. From (4.5), we conclude that k ≤ 1 m-k+1

. From (4.5), we also conclude that k ≤ tol mχ k for all 1 ≤ k ≤ m. By the expression for ∆ m in (4.3), it follows that

In Theorems 3.2 and 3.5 we derived a decaying behavior of the entries in the first column of K -1 m f (A m ). Since these entries decrease in modulus with the row index, it follows from (4.6) that the tolerance of the inexact linear solves can be relaxed with the RKSM progress. In practice, the upper bounds for e * k K -1 m f (A m )βe 1 suggested by Theorem 3.2 or Theorem 3.5 involve integrals, which may take some time to be approximated to a reasonable accuracy. Also, these bounds could give significant overestimates of the actual entries at certain RKSM steps, which may lead to an excessively conservative relaxation estimate.

Instead, we consider a heuristic estimate of e * k K -1 m f (A m )βe 1 based on R k , which usually gives a less conservative tolerance relaxation for the approximate linear solve at each RKSM step and works well in practice. To derive this heuristic, we define the actual entry of interest

, where K m , A m ∈ R m×m are obtained after applying the temporary pole s m = ∞ at step m of the RKSM. From the expression of R m in (2.9),

, where the scalar h k+1,k and the restricted matrices K k , A k ∈ R k×k are obtained after applying the finalized finite pole s k at step k. From (3.1) in Theorem 3.2 or (3.6) in Theorem 3.5, we have

and the first k -1 poles remain the same for generating K k , A k , h k+1,k , and K k , A k , h k+1,k , the definitions of c j in (4.7) and (4.8) have identical expressions, for both analytic functions and Markov functions. Suppose that e * k K -1 m and e * k K -1 k are close and that the bounds in (4.7) and (4.8) are comparably sharp. Then χ k can be approximated by R k |hk+1,k| . Based on Algorithm 1, it is possible to get both R k and h k+1,k with the temporary infinite pole at step k, before we need to use χ k ≈ R k |hk+1,k| to set the tolerance for the iterative linear solve (I -A/s k )w k+1 = (A -σ k I)v k with the finalized finite pole s k . Since this is a heuristic estimate of χ k , we may have a lower risk of applying excessive relaxation by slightly increasing this estimate so that the tolerance of the inexact linear solve can be tightened moderately, and the inexact RKSM may follow the behavior of the exact algorithm more reliably. In practice, we set χ k = 10 R k |hk+1,k| in our numerical tests. ETNA Kent State University and Johann Radon Institute (RICAM)

INEXACT RKSM FOR APPROXIMATING THE ACTION OF FUNCTIONS OF MATRICES 553 5. Numerical experiments. In this section, we first provide numerical evidence to show the sharpness of our upper bounds for the entries of K -1 m f (A m ) for both analytic functions and Markov functions and discuss efficient quadrature rules to approximate these bounds. Then we numerically show the advantage of the inexact RKSM over the exact method for approximating f (A)b. All experiments were carried out in MATLAB R2021b on a laptop running in Windows 10 with 16GB DDR4 2400 MHz memory and a 2.81 GHz Intel Dual Core CPU.

To study the decaying pattern for the residual norm R m 2 in (2.9), we are mostly interested in approximating the upper bound for e * k K -1 m f (A m )e 1 , for k ≤ m, accurately and efficiently. For analytical functions, Theorem 3.2 provides upper bounds in both (3.1) and (3.3), referred to as the original bound and the simplified bound, respectively. Although |c j | defined in (3.2) is independent of the values of τ , numerical tests show that it is unstable to approximate |c j | in (3.2) for a wide range of values of τ . A better approach is to use fminbnd in MATLAB to find the optimal τ > 1 that minimizes the partial sum of the infinite series of the |c j |'s. To compute |c j | in (3.2) with a fixed value of τ , we use a composite trapezoid rule to approximate the integral in (3.2) and then employ the fast Fourier transform (FFT) to evaluate a partial sum of the infinite series in (3.1). In a neighborhood of the optimal value of τ , numerical experiments show the efficiency of the composite trapezoid rule evaluated by the FFT.

Specifically, we evaluate the expression for c j in (3.2) by the composite trapezoid rule

where N is the number of quadrature points, θ p = 2πp N (0 ≤ p ≤ N -1) are the quadrature nodes, w p = 2π N -1 are the quadrature weights for all 0 ≤ p ≤ N -1, and

To approximate the infinite series for |c j |, we compute the first N = 2 14 = 16384 terms of c j . It follows that

where we use MATLAB's fft to compute all P j = N -1 p=0 w p Q p e -i 2πp N j (0 ≤ j ≤ N -1). Similarly, for the simplified bound in (3.3), we also use the composite trapezoid rule to approximate the integral and then call fminbnd in MATLAB to find the optimal τ > 1 .

For Markov functions, Theorem 3.5 provides the original bound in (3.6) and also the simplified bound in (3.7). Since the original bound in (3.6) involves an infinite series, we approximate it by computing the first N = 2 12 = 4096 terms. For the test function f 2 (z) = e - √ z , ETNA Kent State University and Johann Radon Institute (RICAM) 554 S. XU AND F. XUE since µ (ζ) defined in (3.5) is oscillatory, we apply integration by substitution and divide the interval of integration into several subintervals. We apply Gauss-Legendre quadrature on the first few subintervals where the quadrature values are relatively large. The remaining subintervals are approximated by a trigonometric integral. For the test function f 3 (z) = z -1/2 , we divide the interval of integration into several subintervals and apply integration by substitution and Gauss-Jacobi quadrature on appropriate subintervals. We omit these technical details but point out that the quadrature can be evaluated accurately with efficiency. Compared with MATLAB's integral with default setting, numerical experiments show that our quadrature for approximating c j is more accurate and faster. EXAMPLE 5.1. Consider a diagonal (symmetric) matrix A ∈ R 100001×100001 with diagonal entries a kk = -cos πk 100000 10 3 -10 -3 2

, for 0 ≤ k ≤ 100000. We compute several iterations of the RKSM for approximating f (A)b with 4 different functions.

We test an analytic hyperbolic sine function f 0 (z) = sin(hz) = e hz -e -hz 2 for h = 0.01, and we use one repeated single pole s = -m/h = -4000 for 40 steps of the RKSM. We also test the analytic exponential function f 1 (z) = e -z and use the repeated single pole s = m = -40 suggested in [53]. Comparisons between the actual values of e * k K -1 m f (A m )e 1 and the two upper bounds in Theorem 3.2 are reported in the upper left and upper right plots in Figure 5.1 for f 0 (z) and f 1 (z), respectively. Both upper bounds give accurate estimates of the actual values. Compared to the bounds for |e * k A m e ℓ | and |e * k f (A m )e ℓ | investigated for analytic functions in [57], our approach with composite trapezoid quadrature and FFT evaluates the upper bounds efficiently with higher accuracy.

For the Markov functions

The eigenvalues of A are located in an ellipse centered at 10 3 +10 -3

= 500.0005, with semimajor axis of length 10 3 -10 -3 2 = 499.9995 lying on the real axis and semi-minor axis of length 10. Similar to Example 5.1, we compute several iterations of the RKSM for approximating f (A)b with 4 functions. From the results shown in Figure 5.2, we can see for this non-Hermitian matrix whose eigenvalues are located in an ellipse that the upper bounds we derived have similar behavior as those for the Hermitian matrix given in Example 5.1.

Then we consider a few practical problems for which the exact method is simulated by an inexact RKSM, where the linear solve at each RKSM step is performed to a fixed high level of accuracy, to compare with the inexact RKSM with the relaxation strategy discussed in Section 4. EXAMPLE 5.3. We consider a non-Hermitian matrix A ∈ R 23560×23560 from the Mar-trixMarket problem af23560 , whose eigenvalues are located on the right half complex plane. We test the analytic function f 1 (z) = e -z and the Markov function f 2 (z) = e - √ z with the adaptive RKSM in [38,Section 4] to approximate f (A)b, where b ∈ R 23560×1 is a vector with standard normally distributed random entries. We apply both the exact RKSM and inexact RKSM to approximate f (A)b. For the exact method, the linear solves are performed by MATLAB's backslash operation, while for the inexact method, the linear systems are solved by a right-preconditioned GMRES(100) method with the relaxation strategy discussed in Theorem 4.1, where tol = 10 -9 , ξ j = tol mχj , and

The preconditioner is the incomplete LU factorization preconditioner with threshold and pivoting (ILUTP) [60, Section 10.4.4, p. 327], using a drop tolerance 0.01. Figure 5.3 illustrates that if we properly set the relaxed accuracy of the approximate linear solve at each RKSM step, we can get the desired residual norm for the inexact method, and the norm of the difference between the residuals of the exact method and the inexact method remains small through the entire process of the RKSM iterations. EXAMPLE 5.4. We test 11 nonsymmetric real matrices, all of which have the entire spectrum strictly in the right half complex plane. Two of these matrices are of the form A = M -1 K, where both M and K are sparse, but A is not formed explicitly. Specifically,

the problems obstacle and plate, which involve matrices of saddle-point structure arising from modeling incompressible fluid flows in 2D domains, are generated by the IFISS package version 3.6 [30]. The obstacle problem is generated with grid parameter 6, using the biquadratic-bilinear (Q2-Q1) element on a stretched rectangular grid, with viscosity parameter ν = 1 175 corresponding to a Reynolds number Re = 2 ν = 350. The plate problem is constructed with grid parameter 7, using the biquadratic-bilinear element on a non-stretched rectangular grid, with viscosity parameter ν = 1 500 that corresponds to a Reynolds number Re = 2 ν = 1000. Both problems give a matrix pair (K, M ), where

with F being the discrete convection-diffusion operator, B T the gradient operator for the pressure, B the divergence operator for the velocity, G the velocity mass matrix, and η = 0.01 so that the n p (the degree of freedom of the pressure space) infinite eigenvalues of

are mapped to 1 η = 100 without changing the finite eigenvalues [16]. These mapped finite eigenvalues are in the deep interior of the spectrum and should have essentially no impact on the convergence of the RKSM for approximating f (A)b with A = M -1 K. Such a matrix pair (K, M ) has been used to study the linear stability of the steady-state solution of the

Navier-Stokes equation by matrix exponentials [58]. For these two problems, where M is not the identity, we in addition let

The two larger problems LinDir2D and LinDir3D are generated from finite difference discretizations of the second-order linear PDE

], respectively. The artificial "wind" v and w are defined as

e -2xy (y 2 + 2 sin(x)) cos(4x + y)(x 3 + 3e -y ) , w 2D = erf(x -y 2 ) 2 + 2 -8 arctan(x 2 cos(y)) + π/2 , and

We use a standard second-order centered finite difference to approximate the first and second derivatives based on a uniform 2 9 × 2 10 mesh grid of Ω 2D and a 2 6 × 2 7 × 2 8 mesh grid of Ω 3D . Both problems are based on Dirichlet boundary condition but with the boundary nodes included in the matrix. This leads to matrices of order (2 9 + 1) × (2 10 + 1) = 525825 for LinDir2D and (2 6 + 1) × (2 7 + 1) × (2 8 + 1) = 2154945 for LinDir3D, respectively. The original matrices corresponding to the linear differential operator of this PDE were scaled by max{h x , h y } 2 and max{h x , h y , h z } 2 , respectively, where h x , h y , and h z are the mesh size in the x, y, and z directions, so that the scaled matrices have bounded norms independent of the mesh size (consistent with the assumption that W (A) ⊂ E and E ⊂ C is compact). An application of MATLAB's backslash to solve a shifted linear system involving the scaled matrices of LinDir2D is still faster than our iterative method, whereas for LinDir3D it used up to 16 GB memory on our machine in a few minutes and led MATLAB to crash (in fact, 32 GB memory was still not sufficient to solve the linear systems involving LinDir3D by backslash). Other matrices are selected from the SuiteSparse Matrix Collection [17].

Kent State University and Johann Radon Institute (RICAM)

To compare the behavior of the inexact RKSM with and without the relaxation strategy for the inner linear solves, we use the RKSM to approximate f (A)b for f 1 (z) = e -z , f 2 (z) = e - √ z , f 3 (z) = z -1/2 and a random vector b whose entries follow a standard normal distribution. For the inexact method, we use the same strategy as in Example 5.3, and for the exact method, we let ξ k = min 1≤i≤k ξ k to simulate the behavior of the ideal exact RKSM that performs an exact linear solve at each step. The adaptive poles of the RKSM are chosen following the strategy adopted in [38,Section 4]. We use the right-preconditioned GMRES(70) method as the inner linear solver for the RKSM. The maximum number of GMRES restart cycles is set to be J = 20. We use the ideal least-squares commutator (LSC) preconditioners [28,29] for the problems obstacle and plate from IFISS. For LinDir2D and LinDir3D, the preconditioner is one W-cycle of the geometric multigrid (GMG) mtehod with two applications of the Gauss-Seidel method as pre-and post-smoothers. For all other matrices (from SuiteSparse), we use the incomplete LU preconditioner with threshold dropping and pivoting (ILUTP) preconditioners [60,Section 10.4.4]; also see MATLAB's documentation for ilu with the option for the ilutp.

The results of the performance of the exact and inexact RKSM for f 1 (z) = e -z , f 2 (z) = e - √ z , and f 3 (z) = z -1/2 are summarized in Tables 5.1, 5.2, and 5.3, respectively.

We show the size of the matrices n, the residual tolerance tol , the type of preconditioners (including the drop tolerance for ILUTP), the max number of RKSM steps m, the total number of GMRES iterations, the runtime for both the exact and the inexact methods, and the number of RKSM steps to converge. We chose the smallest tolerance tol (a negative integer power of 10) for each test matrix across all test functions such that the inexact RKSM can successfully converge to this tolerance for all functions of interest here. Such a problem-dependent tolerance is preferred to a uniform tolerance because an absolute tolerance for the residual (2.6) depends on the norm (and probably the condition number) of the matrix A. A uniform absolute tolerance similarly does not indicate the quality of approximation for different problems, nor does it show whether the computed approximation is close to the most accurate approximation achievable in double precision. Note that the total runtime includes the time for constructing the preconditioners, applying GMRES, the orthogonalization of the basis vectors of the RKSM, evaluating f (A m ) for the small restricted matrices, and estimating the level of relaxation for the inner linear solves. Overall, from the results in Tables 5.1-Table 5.3, the inexact methods need fewer GMRES iterations to solve the inner linear systems, so that they need less time to converge than the "exact" RKSM. However, the level of advantage of the inexact RKSM over the exact method varies for different matrices and preconditioners. In general, if the proportion of time used to construct preconditioners is small, then the relative advantage of the inexact RKSM is significant.

To demonstrate this point, we choose f 2 (z) = e - √ z as an example and show the computation time used for constructing the preconditioners and applying GMRES in Table 5.4.

For the problem af23560, the exact RKSM needs 47% of the total runtime to construct the preconditioners, and the inexact RKSM requires 35% less runtime than the exact method; by comparison, for the problem LinDir3D, the exact RKSM takes only 1% of the total runtime to construct the preconditioners, and the inexact RKSM takes 79% less runtime than the exact method. In general, excluding the cost for constructing the preconditioners, the inexact RKSM can save about 35-81% of the runtime needed for the exact method.

6. Proofs. In this section, we prove Theorems 3.2 and 3.5 in Section 3. We use a rational approximation approach called the Faber-Dzhrbashyan (FD) rational functions introduced in [25]; see also [62, Ch. XIII, Section 3] and the references therein. Our proofs of Lemma A.1 and Theorem 3.2 largely follow the ideas of those in [57], especially the introduction to FD rational functions before the proofs of Theorems 3.2 and 3.5. To make our paper self-contained, the background details are presented in Appendix A. There are two minor differences between the materials in the appendix and in [57, Section 7]. First, a more complete description of the conditions for the expansions of the FD rational functions is presented in the appendix based on the original reference [62]. Second, a sharper upper bound is constructed in Lemma A.1 by using the least number of inequalities. In the proofs of Theorems 3.2 and 3. our derivation is developed for the entries of K -1 m f (A m ) instead of A m or f (A m ), and we keep the integrals of the upper bounds to provide sharper bounds and propose efficient quadrature rules to evaluate them accurately. In addition, to the best of our knowledge, Theorem 3.5 for Markov functions and its proof in this paper are new. , for 1 ≤ j ≤ k -ℓ, and in addition we set

yielding the following expansion:

By Lemma 2.4, we have

Note that since we set α j = 0, for j > k -ℓ, from (A.1) with j ≥ k -ℓ, we have

.3) ETNA Kent State University and Johann Radon Institute (RICAM)

Combining c j in (6.1) and (6.3), for j ≥ k -ℓ, it holds that

Since we set α j = 0, for j > k -ℓ, (A.11) implies (6.5)

Note that similar to (A.11), the bound in (6.5) is valid for both analytic functions and Markov functions. Combining (6.2), (6.4), and (6.5), we get

The simplified bound in (3.3) can be derived as follows:

We emphasize that c j is independent of τ > 1. In fact, from (6.4), if we define

the residue theorem shows that

where w i (1 ≤ i ≤ r) denote all the poles of g in the set {w : |w| ≤ τ }. Since τ > 1 and |α i | < 1, these poles are 0, α 1 , . . . , α k-ℓ , which means that Res(g, w i ) is independent of the value of τ for 1 ≤ i ≤ r. Therefore, the coefficients c j are also independent of the value of τ . This independence of τ is important for us to develop reliable FFT-based composite trapezoid quadrature rule to evaluate these coefficients for analytic functions. From the definition of ϕ j+1 1 w in (6.3) and c j in (6.6), we have (6.7)

dµ(ψ(w)) ψ (w) .

Combining the upper bounds for M j (A) in (6.5) and |c j | in (6.7) with (6.2), it holds that

In this paper, we studied the residual of the RKSM for approximating the action of a function of a matrix f (A) to a vector b. We explored the decay bounds for the off-diagonal entries of a restricted matrix that arise in the RKSM approximation of f (A)b for analytic functions and Markov functions. For the inexact RKSM, upper bounds for the allowable errors for the inner linear solves are derived, and a heuristic tolerance relaxation strategy is proposed to enable that the inexact RKSM keeps track of the convergence of the exact RKSM. Numerical experiments show that the inexact RKSM can exhibit a convergence behavior similar to that of the exact method, but it entails lower computational cost thanks to the relaxed accuracy for the inner linear systems.

Appendix A. This appendix provides an introduction to the Faber-Dzhrbashyan (FD) rational functions for Theorems 3.2 and 3.5. To this end, we first review the definition of the Takenaka-Malmquist (TM) system of rational functions: where δ mn is the Kronecker delta; see, e.g., [54].

The FD rational function M j (z) is defined as the sum of the principal part and the constant in the Laurent decomposition of ϕ j (φ(z)) in the neighborhoods of the points {z k } Proof. First, if j = 1, we have from (A.1) that

For j > 1 and |w| = τ , it holds that (A.4)

Let w = τ e θ0i and α j = ρ j e θj i for some ρ j = |α j | ∈ [0, 1). We define g j (θ 0 ) := |w -α j | = τ e -θ0i -ρ j e -θj i = τ e (θj -θ0)i -ρ j = τ 2 + ρ 2 j -2τ ρ j cos (θ j -θ 0 ) ≥ τ -ρ j .

it is easy to show that h k (θ 0 ) achieves its maximum when θ k -θ 0 = π. Therefore,

Combining (A.5) and (A.6) into (A.4), we get

For j = 1, it is easy to verify that the above inequality also holds. We can write the FD rational functions M j (z) with the Faber transformation of a Takenaka-Malmquist system. For every function f continuous on the boundary of D and analytic in the interior of D, the Faber transformation is defined as