Rational Krylov for Stieltjes matrix functions: convergence and pole selection

Evaluating the action of a matrix function on a vector, that is $x=f(\mathcal M)v$, is an ubiquitous task in applications. When $\mathcal M$ is large, one usually relies on Krylov projection methods. In this paper, we provide effective choices for the poles of the rational Krylov method for approximating $x$ when $f(z)$ is either Cauchy-Stieltjes or Laplace-Stieltjes (or, which is equivalent, completely monotonic) and $\mathcal M$ is a positive definite matrix. Then, we consider the case $\mathcal M=I \otimes A - B^T \otimes I$, and $v$ obtained vectorizing a low-rank matrix, which finds application, for instance, in solving fractional diffusion equation on two-dimensional tensor grids. We see how to leverage tensorized Krylov subspaces to exploit the Kronecker structure and we introduce an error analysis for the numerical approximation of $x$. Pole selection strategies with explicit convergence bounds are given also in this case.


Introduction.
We are concerned with the evaluation of x = f (M)v, where f (z) is a Stieltjes function, which can be expressed in integral form (1.1) f (z) = ∞ 0 g(z, t)µ(t) dt, g(z, t) ∈ e −zt , 1 z + t .
The two choices for g(z, t) define Laplace-Stieltjes and Cauchy-Stieltjes functions, respectively [9,28]. The former class is a superset of the latter, and coincides with the set of completely monotonic functions, whose derivatives satisfy (−1) j f (j) 0 over R + for any j ∈ N.
We are interested in two instances of this problem; first, we consider the case M := A, where A ∈ C n×n is Hermitian positive definite with spectrum contained in [a, b], v ∈ C n is a generic vector, and a rational Krylov method [16] is used to approximate x = f (M)v. In this case, we want to estimate the Euclidean norm of the error x − x ℓ 2 , where x ℓ is the approximation returned by ℓ steps of the method. Second, we consider (1.2) M := I ⊗ A − B T ⊗ I ∈ C n 2 ×n 2 , where A, −B ∈ C n×n are Hermitian positive definite with spectra contained in [a, b], v = vec(F ) ∈ C n 2 is the vectorization of a low-rank matrix F = U F V T F ∈ C n×n , and a tensorized rational Krylov method [9] is used for computing vec(X) = f (M)vec(F ). This problem is a generalization of the solution of a Sylvester equation with a low-rank right hand side, which corresponds to evaluate the function f (z) = z −1 . Here, we are concerned with estimating the quantity X − X ℓ 2 , where X ℓ is the approximation obtained after ℓ steps.
1.1. Main contributions. This paper gives the following main contributions. We provide a convergence analysis that links the norm of the error for a given choice of poles to a parameter dependent rational approximation (with the given poles) of the kernel functions e −zt and 1 z+t . Table 1.1 Summary of the convergence rates for rational Krylov methods with the proposed poles. The convergence rate ρ [α,β] is defined by ρ [α,β] := exp(−π 2 / log(4 β α )).
These problems can be related -by relying on some Möbius transforms -to the so-called Zolotarev problem on symmetric real intervals [30]. This, in turn, allows us to give explicit expressions of poles that achieve a quasi-optimal convergence rate for the two classes of interest. These results, summarized in Table 1.1, are given for both the standard rational Krylov method and the tensorized version (for M = I ⊗ A − B T ⊗ I). Our analysis for the Cauchy case without Kronecker structure provides analogous results to the ones in [4, Section 6.1] which relies on powerful tools from logarithmic potential theory. See also [14] for similar, although less explicit, results.
It is well known that completely monotonic functions are well approximated by exponential sums [12]. Our results prove constructively that they are also well-approximated by rational functions.
A byproduct of our analysis is that, in the Kronecker structured case, the low-rank property of the right hand side numerically guarantees the low-rank approximability of the matrix X. This generalizes the well-known property of the solutions of of Sylvester equations with low-rank right hand sides [6]. In particular, we provide bounds for the exponential decay for the singular values of X, see Section 4.4.

Motivating problems.
Computing the action of a matrix function on a vector is a classical task in numerical analysis, and find applications in several fields, such as complex networks [8], signal processing [25], numerical solution of ODEs [17], and many others.
Matrices with the Kronecker sum structure as in (1.2) often arise from the discretization of differential operators on tensorized 2D grids. Applying the inverse of such matrices to a vector is equivalent to solving a matrix equation. When the right hand side is a smooth function or has small support, the vector v is the vectorization of a numerically low-rank matrix. The latter property has been exploited to develop several efficient solution methods, see [24] and the references therein. Variants of these approaches have been proposed under weaker assumptions, such as when the smoothness is only available far from the diagonal x = y, as it happens with kernel functions [19,21].
In recent years, there has been an increasing interest in models involving fractional derivatives. For 2D problems on rectangular grids, discretizations by finite differences or finite elements lead to linear systems that can be recast into the solution of matrix equations with particularly structured coefficients [13,20]. However, a promising formulation which simplifies the design of boundary conditions relies on first discretizing the 2D Laplacian on the chosen domain, and then consider the action of the matrix function z −α (with the Laplacian as argument) to the right hand side. This is known in the literature as matrix transform method [29]. In this framework, one has 0 < α < 1, and therefore z −α is a Cauchy-Stieltjes function, a property that has been previously exploited for designing fast and stable restarted polynomial Krylov methods for its evaluation [23]. The algorithm proposed in this paper allows to exploit the Kronecker structure of the 2D Laplacian on rectangular domains in the evaluation of the matrix function.
Another motivation for our analysis stems from the study of exponential integrators, where it is required to evaluate the ϕ j (z) functions [17], which are part of the Laplace-Stieltjes class. This has been subject to a deep study concerning (restarted) polynomial and rational Krylov methods [15,23]. However, to the best of our knowledge the Kronecker structure, and the associated low-rank preservation, has not been exploited in these approaches, despite being often present in discretization of differential operators [26].
2. Laplace-Stieltjes and Cauchy-Stieltjes functions. We recall the definition and the properties of Laplace-Stieltjes and Cauchy-Stieltjes functions that are relevant for our analysis. Functions expressed as Stieltjes integrals admit a representation of the form: where µ(t)dt is a (non-negative) measure on [0, ∞], and g(z, t) is integrable with respect to that measure. The choice of g(z, t) determines the particular class of Stieltjes functions under consideration (Laplace-Stieltjes or Cauchy-Stieltjes), and µ(t) is called the density of f (z). µ(t) can be a proper function, or a distribution, e.g. a Dirac delta. We refer the reader to [28] for further details.
Examples of Laplace-Stieltjes functions include: The last example is an instance of a particularly relevant class of Laplace-Stieltjes functions, with applications to exponential integrators. These are often denoted with ϕ j (z), and can be defined as follows: A famous theorem of Bernstein states the equality between the set of Laplace-Stieltjes functions and those of completely monotonic functions in (0, ∞) [11], that is (−1) j f (j) (z) is positive over [0, ∞] for any j ∈ N.
From the algorithmic point of view, the explicit knowledge of the Laplace density µ(t) will not play any role. Therefore, for the applications of the algorithms and projection methods described here, it is only relevant to know that a function is in this class.
A few examples of Cauchy-Stieltjes functions are: The rational functions with poles on the negative real semi-axis do not belong to this class if one requires µ(t) to be a function, but they can be obtained by setting µ(t) = h j=1 α j δ(t−β j ), where δ(·) is the Dirac delta with unit mass at 0. For instance, z −1 is obtained setting µ(t) := δ(t).
Since Cauchy-Stieltjes functions are also completely monotonic in (0, ∞) [10], the set of Cauchy-Stieltjes functions is contained in the one of Laplace-Stieltjes functions. Indeed, assuming that f (z) is a Cauchy-Stieltjes function with density µ C (t), one can construct a Laplace-Stieltjes representation as follows: where µ L (s) defines the Laplace-Stieltjes density. In particular, note that if µ C (t) is positive, so is µ L (s). For a more detailed analysis of the relation between Cauchy-and Laplace-Stieltjes functions we refer the reader to [28,Section 8.4].
As in the Laplace case, the explicit knowledge of µ(t) is not crucial for the analysis and is not used in the algorithm.
3. Rational Krylov for evaluating Stieltjes functions. Projection schemes for the evaluation of the quantity f (A)v work as follows: an orthonormal basis W for a (small) subspace W ⊆ C n is computed, together with the projections A W := W * AW and v W := W * v. Then, the action of f (A) on v is approximated by: Intuitively, the choice of the subspace W is crucial for the quality of the approximation. Usually, one is interested in providing a sequence of subspaces W 1 ⊂ W 2 ⊂ W 3 ⊂ . . . and study the convergence of x Wj to f (A)v as j increases. A common choice for the space W j are Krylov subspaces.
3.1. Krylov subspaces. Several functions can be accurately approximated by polynomials. The idea behind the standard Krylov method is to generate a subspace that contains all the quantities of the form p(A)v for every p(z) polynomial of bounded degree.
Definition 3.1. Let A be an n × n matrix, and v ∈ C n×s be a (block) vector. The Krylov subspace of order ℓ generated by A and v is defined as Projection on Krylov subspaces is closely related to polynomial approximation. Indeed, if f (z) is well approximated by p(z), then p(A)v is a good approximation of f (A)v, in the sense Rational Krylov subspace are their rational analogue, and can be defined as follows.
Definition 3.2. Let A be a n×n matrix, v ∈ C n×s be a (block) vector and Ψ = (ψ 1 , . . . , ψ ℓ ), with ψ j ∈ C := C ∪ {∞}. The rational Krylov subspace with poles Ψ generated by A and v is defined as where q ℓ (z) = ℓ j=1 (z − ψ j ) and if ψ j = ∞, then we omit the factor (z − ψ j ) from q ℓ (z). Note that, a Krylov subspace is a particular rational Krylov subspace where all poles are chosen equal to ∞: RK ℓ (A, v, (∞, . . . , ∞)) = K ℓ+1 (A, v). A common strategy of poles selection consists in alternating 0 and ∞. The resulting vector space is known in the literature as the extended Krylov subspace.
It is well-known that Krylov subspaces contain the action of related rational matrix functions of A on the (block) vector v, if the poles of the rational functions are a subset of the poles used to construct the approximation space. Lemma 3.3 (Exactness property). Let A be a n × n matrix, v ∈ C n×s be a (block) vector and Ψ = {ψ 1 , . . . , ψ ℓ } ⊂ C. If U P , U R are orthonormal bases of K ℓ (A, v) and RK ℓ (A, v, Ψ), respectively, then: The optimal choice of poles for generating the rational Krylov subspaces is problem dependent and linked to the rational approximation of the function f (z) on [a, b]. We investigate how to perform this task when f is either a Laplace or Cauchy-Stieltjes function.
3.2. Simultaneous approximation of resolvents and matrix exponentials. The integral expression (1.1) reads as when evaluated at a matrix argument. Since the projection is a linear operation we have This suggests to look for a space approximating uniformly well, in the parameter t, matrix exponentials and resolvents, respectively. The following result -concerning the approximation of the resolvent -will be key in our construction, and is obtained adapting some ideas by Oseledets in [22], also used by Druskin, Knizhnerman and Zaslavsky in [14].
Theorem 3.4. Let A be Hermitian positive definite with spectrum contained in [a, b] and U be an orthonormal basis of U R = RK ℓ (A, v, Ψ). Then, ∀t ∈ [0, ∞), we have the following inequality: Proof. Following the construction in [22], we consider the function g(λ, t) defined by where (λ i , t j ) is a ℓ × ℓ grid of interpolation nodes. The function g(λ, t) is usually called Skeleton approximation and it is practical for approximating 1 λ+t ; indeed its relative error takes the explicit form: z−λj z+tj . If this ratio of rational functions is small, g(λ, t) is a good approximation of 1 λ+t and -consequently -g(A, t) is a good approximation of (tI + A) −1 . Note that, for every fixed t, g(λ, t) is a rational function in λ with poles −t 1 , . . . , −t ℓ . Therefore, using the poles ψ j = −t j , j = 1, . . . , ℓ for the projection we may write, thanks to the exactness property of Lemma 3.3: Taking the minimum over the possible choices of the parameters λ j we get (3.2).
Concerning the rational approximation of the (parameter dependent) exponential, the idea is to rely on its Laplace transform that involves the resolvent: In this formulation, it is possible to exploit the Skeleton approximation of 1 s+λ in order to find a good choice of poles, independently on the parameter t. For proving the main result we need the following technical lemma whose proof is given in the appendix.
Proof. We consider the Skeleton approximation of 1 λ+s by restricting the choice of poles in both variables to Ψ Then, by using (3.3) for A and A ℓ we get Adding and removing the term g(A, s)v = U g(A ℓ , s)U * v (where the equality holds thanks to Lemma 3.3) we obtain the error expression Since A and A ℓ are normal, the above integrals can be controlled by the maximum of the corresponding scalar functions on the spectrum of A (and A ℓ ), which yields the bound We note that r Ψ (λ) can be pull out of the integral, since it does not depend on s, and thus the above can be rewritten as where p(s) is as in Lemma 3.5 and δ(t) indicates the Dirac delta function. Since the 1-norm of Therefore, we need to estimate the infinity norm of L −1 1 s p(s) p(−s) (t). Such inverse Laplace transform can be uniformly bounded in t using Lemma 3.5 with a constant that only depends on ℓ and b/a: This completes the proof.
Remark 3.7. The constant provided by Lemma 3.5 is likely not optimal. Indeed, experimentally it seems to hold that γ ℓ,κ = 1 for any choice of poles in the negative real axis -not necessarily contained in [−b, −a] -and this has been verified on many examples. If this is proved, then the statement of Theorem 3.4 can be made sharper by removing the factor γ ℓ,κ .
3.3. Bounds for the rational approximation problems. Theorem 3.4 and Theorem 3.6 show the connection between the error norm and certain rational approximation problems. In this section we discuss the optimal values of such problems in the cases of interests.
Definition 3.8. Let Ψ ⊂ C be a finite set, and I 1 , I 2 ⊆ C. Then, we define 1 The θ ℓ functions enjoy some invariance and inclusion properties, which we report here, and will be extensively used in the rest of the paper. Lemma 3.9. Let I 1 , I 2 be subsets of the complex plane, and Ψ ⊂ C. Then, the map θ ℓ satisfies the following properties: is monotonic with respect to the inclusion on the parameters I 1 and I 2 : Proof. Property (i) follows by (iii) because applying a shift is a particular Möbius transformation. Note that, generically, when we compose a rational function r(z) = p(z) h(z) ∈ P ℓ Ψ with M −1 (z) we obtain another rational function of (at most) the same degree and with poles M (Ψ). Hence, we obtain Property (ii) follows easily from the fact that the maximum taken on a larger set is larger, and the minimum taken on a larger set is smaller. Now, we consider the related optimization problem, obtained by allowing Ψ to vary: The latter was posed and studied by Zolotarev in 1877 [30], and it is commonly known as the third Zolotarev problem. We refer to [3] for a modern reference where the theory is used to recover bounds on the convergence of rational Krylov methods and ADI iterations for solving Sylvester equations.
In the case which admits the following explicit estimate.
In addition, the optimal rational function r Explicit expression for the elements of Ψ We use Theorem 3.10 and the Möbius invariance property as building blocks for bounding ] -for some a ∈ (0, 1)with a Möbius transformation; then make use of Theorem 3.10 and Lemma 3.9-(iii) to provide a convenient choice of Ψ for the original problem.
Remark 3.12. We note that the estimate ρ [ a, 1] ρ [a,4b] is asymptotically tight, that is the limit of ρ [ a,1] /ρ [a,4b] → 1 as b/a → ∞. For instance, if b/a = 10 the relative error between the two quantities is about 2 · 10 −2 , and for b/a = 1000 about 5 · 10 −5 . Since the interest for this approach is in dealing with matrices not well-conditioned, we consider the error negligible in practice.
In light of Theorem 3.10 and Lemma 3.11, we consider the choice in Theorem 3.4. This yields the following estimate.
When considering Laplace-Stieltjes functions, we may choose as poles Ψ Proof. The proof relies on the fact that the optimal Zolotarev function evaluated on the interval [a, b] can be bounded by 2ρ Since its zeros and poles are symmetric with respect to the imaginary axis and real, we can apply Theorem 3.6 and conclude.

Convergence bounds for Stieltjes functions.
Let us consider f (z) Stieltjes function of the general form (1.1). Then the error of the rational Krylov method for approximating f (A)v is given by where g(A, t) is either a parameter dependent exponential or resolvent functions. Therefore Corollary 3.13 and Corollary 3.14 provide all the ingredients to study the error of the rational Krylov projection, when the suggested pole selection strategy is adopted. ) and . where γ ℓ,κ is defined as in Theorem 3.6.
Proof. Since f (z) is a Laplace-Stieltjes function, we can express the error as follows: where we applied (3.1) and Corollary 3.14 to obtain the second inequality.
Remark 3.16. Corollary 3.15 requires the function f (z) to be finite over [0, ∞), which might not be the case in general (consider for instance x −α , which is both Cauchy and Laplace-Stieltjes). Nevertheless, the result can be applied to f (z + α), which is always completely monotonic for a positive α, by taking 0 < α < a. A value of α closer to a gives a slower decay rate, but a smaller constant. Similarly, if f (z) happens to be completely monotonic on an interval larger than [0, ∞), bounds with a faster asymptotic convergence rate (but a larger constant) can be obtained considering α < 0.
Example 3.18. We test the proposed scheme for the evaluation of the Laplace-Stieltjes function f (z) = ϕ 1 (z) that arises in the use of exponential integrators. We consider the 1D diffusion problem over [0, 1] with zero Dirichlet boundary conditions discretized using central finite differences in space with 1000 points, and integrated by means of the exponential Euler method. We choose the diffusion coefficient as K = 10 −2 , and a time step of ∆t = 0.1. This requires to evaluate ϕ 1 ( 1 h 2 ∆tA)v, where A is the tridiagonal matrix A = tridiag(−1, 2, −1), and v is a vector. To test the performance of the algorithm, we use a random vector, and we consider the poles from Corollary 3.14, the extended Krylov method, and polynomial Krylov. The results are reported in Figure 3.1. It is visible that, since A is ill-conditioned (the condition number is about 5 · 10 5 ), polynomial Krylov performs poorly. On the other hand, extended Krylov is a valid competitor. As expected the choice of poles from Corollary 3.14 yields the fastest convergence rate.
Example 3.19. To test the performances of the poles for Cauchy-Stieltjes functions, we consider the evaluation of f (z) = z − 1 4 to the 1D Laplacian used also in Example 3.18. The tested matrix is 1000 × 1000, and the results are reported in Figure 3.2. The plot shows that, also in this case, polynomial Krylov has a very slow convergence rate (due to the high condition number of the matrix), whereas the poles for Cauchy-Stieltjes functions offer the best convergence rate. We have also tested the ones for Laplace-Stieltjes (since Cauchy-Stieltjes functions are also Laplace-Stieltjes), which converge at about half the speed, as predicted. Extended Krylov is again a possible alternative, but in this case there is a larger gain in convergence speed by choosing the poles Ψ C,ℓ from Corollary 3.13. 4. Evaluating Stieltjes functions of matrices with Kronecker structure. We consider the task of computing f (M)v where M = I ⊗ A − B T ⊗ I. This problem often stems from the discretizations of 2D differential equations, such as the matrix transfer method used for fractional diffusion equations [29].
We assume that v = vec(F ), where F = U F V T F where U F and V F are tall and skinny matrices.

Iterations (ℓ)
x − x ℓ 2 Bound from Corollary 3.14 Extended Krylov Polynomial Krylov Poles from Corollary 3.14 numerically inherited by X [6]. For more general functions than z −1 , a projection scheme that preserves the Kronecker structure has been proposed in [9] using polynomial Krylov methods. We briefly review it in Section 4.1. The method proposed in [9] uses tensorized polynomial Krylov subspaces, so it is not well-suited when A and B are ill-conditioned, as it often happens discretizing differential operators. Therefore, we propose to replace the latter with a tensor product of rational Krylov subspaces and we provide a strategy for the pole selection. This enables a faster convergence and an effective scheme for the approximation of the action of such matrix function in a low-rank format.
The case of Laplace-Stieltjes functions, described in Section 4.2, follows easily by the analysis performed for the pole selection with a generic matrix A. The error analysis for Cauchy-Stieltjes functions, presented in Section 4.3, requires more care and builds on the theory for the solution of Sylvester equations.
4.1. Projection methods that preserve Kronecker structure. If A, B are n×n matrices, applying the projection scheme described in Section 3 requires to build an orthonormal basis W for a (low-dimensional) subspace W ⊆ C n 2 , together with the projections of W * MW = H and v W = W * v. Then the action of f (M) on v is approximated by: The trick at the core of the projection scheme proposed in [9] consists in choosing a tensorized subspace of the form W := U ⊗ V, spanned by an orthonormal basis of the form W = U ⊗ V , where U and V are orthonormal bases of U ⊆ C n and V ⊆ C n , respectively. With this choice, the projection of M onto U ⊗ V retains the same structure, that is Since in our case v = vec(F ) and F = U F V T F , this enables to exploit the low-rank structure as well. Indeed, the projection of F onto U ⊗ V can be written as v W = vec(F W ) = vec((U * U F )(V T F V )). The high-level structure of the procedure is sketched in Algorithm 4.1.
U, V ← orthonormal bases for the selected subspaces.

2:
A U ← U * AU 3: At the core of Algorithm 4.1 is the evaluation of the matrix function on the projected matrix I ⊗ A U − B T V ⊗ I. Even when U, V have a low dimension k ≪ n, this matrix is k 2 × k 2 , so it is undesirable to build it explicitly and then evaluate f (·) on it.
When f (z) = z −1 , it is well-known that such evaluation can be performed in k 3 flops by the Bartels-Stewart algorithm [2], in contrast to the k 6 complexity that would be required by a generic dense solver for the system defined by I ⊗ A U − B T V ⊗ I. For a more general function, we can still design a O(k 3 ) procedure for the evaluation of f (·) in our case. Indeed, since A U and B V are Hermitian, we may diagonalize them using a unitary transformation as follows: Then, the evaluation of the matrix function f (z) with argument I ⊗ A U − B T V ⊗ I can be recast to a scalar problem by setting where • denotes the Hadamard product and f • (·) the function f (·) applied component-wise to the entries of D: Assuming that the matrices Q A , Q B and the corresponding diagonal matrices D A , D B , are available, this step requires k 2 scalar function evaluation, plus 4 matrix-matrix multiplications, for a total computational cost bounded by O(c f · k 2 + k 3 ), where c f denotes the cost of a single function evaluation. The procedure is described in Algorithm 4.2.
for i, j = 1, . . . , n do 6: end for 8: return vec(Q A XQ * B ) 9: end procedure 4.2. Convergence bounds for Laplace-Stieltjes functions of matrices with Kronecker structure. The study of approximation methods for Laplace-Stieltjes functions is made easier by the following property of the matrix exponential: whenever M, N commute, then e M+N = e M e N . Since the matrices B T ⊗ I and I ⊗ A commute, we have Consider projecting the matrix M onto a tensorized subspace spanned by the Kronecker products of unitary matrices U ⊗ V . This, combined with Algorithm 4.1, yields an approximation whose accuracy is closely connected with the one of approximating e −tA by projecting using U , and e tB using V . As discussed in Section 3, there exists a choice of poles that approximates uniformly well the matrix exponential, and this can be leveraged here as well.
Proof. If f (z) is a Laplace-Stieltjes function, we may express the error matrix X − X ℓ as follows: where A ℓ = U * AU and B ℓ = V * BV . Adding and subtracting the quantity U e −tA ℓ U * F e tB yields the following inequalities: where in the last step we used Corollary 3.14 for both addends.
Example 4.2. To test the algorithm we consider the same matrix A of Example 3.18, and we evaluate the matrix function to M = I ⊗ A + A ⊗ I, applied to a vector v = vec(F ), where F is a random rank 1 matrix, generated by taking the outer product of two unit vectors with normally distributed entries.

Convergence bounds for
Cauchy-Stieltjes functions of matrices with Kronecker structure. As already pointed out in Section 3, evaluating a Cauchy-Stieltjes function requires a space which approximates uniformly well the shifted inverses of the matrix argument under consideration. When considering a matrix M = I ⊗ A − B T ⊗ I which is Kronecker structured, this acquires a particular meaning.
In fact, relying on the integral representation (2.1) of f (z) we obtain: where X t := vec −1 ((tI + M) −1 vec(F )) solves the matrix equation Therefore, to determine a good projection space for the function evaluation, we should aim at determining a projection space where these parameter dependent Sylvester equations can be solved uniformly accurate. We note that, unlike in the Laplace-Stieltjes case, the evaluation of the resolvent does not split into the evaluation of the shifted inverses of the factors, and this does not allow to apply Theorem 3.4 for the factors A and B. A possible strategy to determine an approximation space is using polynomial Krylov subspaces K m (A + tI, U F ) ⊗ K m (B T , V F ) for solving (4.1) at a certain point t. Thanks to the shift invariance of polynomial Krylov subspaces, all these subspaces coincide with U P ⊗ V P = K m (A, U F ) ⊗ K m (B T , V F ). This observation is at the core of the strategy proposed in [9], which makes use of U P ⊗ V P in Algorithm 4.1. This allows to use the convergence theory for linear matrix equations to provide error bounds in the Cauchy-Stieltjes case, see [9, Section 6.2].
Since rational Krylov subspaces are usually more effective in the solution of Sylvester equations, it is natural to consider their use in place U P ⊗ V P . However, they are not shift invariant, and this makes the analysis not straightforward. Throughout this section, we denote by U ⊗ V the orthonormal basis of the tensorized rational Krylov subspace where Ψ := {ψ 1 , . . . , ψ ℓ } and Ξ := {ξ 1 , . . . , ξ ℓ } are the prescribed poles. We define the following polynomials of degree (at most) ℓ: and we denote by A ℓ = U * AU , B ℓ = V * BV the projected (ℓk × ℓk)-matrices, where k is the number of columns of U F and V F . In Section 4.3.1, we recall and slightly extend some results about rational Krylov methods for Sylvester equations i.e., the case f (z) = z −1 . This will be the building block for the convergence analysis of the approximation of Cauchy-Stieltjes functions in Section 4.3.2.

4.3.1.
Convergence results for Sylvester equations. Algorithm 4.1 applied to f (z) = z −1 coincides with the Galerkin projection method for Sylvester equations [24], whose error analysis can be found in [3]; the results in that paper relate the Frobenius norm of the residual to a rational approximation problem. We state a slightly modified version of Theorem 2.1 in [3], that enables to bound the residual in the Euclidean norm.
In the proof of [3, Theorem 2.1] it is shown that, for any choice of (ℓ, ℓ)-rational function r B (x) with poles z 1 , . . . , z ℓ , we can further decompose ρ 12 as ρ 12 = U (J 1 − J 2 ), where with S A,B (X) := AX − XB and Γ A a path encircling once the interval I A but not I B .
With a direct integration we get In the original statement of [3, Theorem 2.1] the residual is decomposed in three parts; the missing term is equal to zero whenever the projection subspace contains the right hand side, which is indeed our case.
where we used that V * B = B ℓ V * and that the integral on the path Γ A of (zI − and consequently Taking the minimum over all (ℓ, ℓ)-rational functions with poles Ξ provides Analogously one obtains the similar estimate for ρ 21 swapping the role of A and B. Since ρ 12 and ρ 21 have orthonormal rows and columns, we have ρ 12 + ρ 21 2 = max{ ρ 12 2 , ρ 21 2 }, which concludes the proof.
Remark 4.4. Using the mixed norm inequality AB F A F B 2 , one can state the bound in the Frobenius norm as well: which is tighter than the one in [3], which involves the Frobenius norm of the evaluated matrix functions.
For our analysis, it is more natural to bound the approximation error of the exact solution X, instead of the norm of the residual. Since the residual is closely related with the backward error of the underlying linear system, bounding the forward error X − X ℓ 2 causes the appearances of an additional condition number.
Proof. We note that X ℓ − X solves the Sylvester equation where κ = b a and θ ℓ (·, ·, ·) is as in Definition 3.8. Proof. Applying the definition of f (M) we have We note that, for any t 0, the vector vec(X t ) := (M + tI) −1 vec(F ) is such that X t solves the Sylvester equation (A + tI)X t − X t B = F . Then, we can write X as Let us consider the approximation U Y t V * to X t obtained by solving the projected Sylvester Then, relying on Corollary 4.5, we can bound the error R t : where C(t) := 2(a+b+t) (2a+t) 2 . Making use of Lemma 3.9-(i) we get: An estimate for the error on X is obtained by integrating R t : where we used that the function 2(a+b+t) 2a+t is maximum over [0, ∞] at t = 0.
Inspired by Theorem 4.3, we look at the construction of rational functions that make the quantities θ ℓ (I A , I B − t, Ψ) and θ ℓ (I B , I A + t, Ξ) small. If we choose Ξ = −Ψ then (4.4) simplifies to (4.5) X The inverse map T (z) −1 is: In addition, we have a −1 2b/a, and therefore ρ [ a,1] ρ [a,2b] .
Proof. The proof can be easily obtained following the same steps of Lemma 3.11. As in that case, the overestimate introduced by the inequality ρ [ a, 1] ρ [a,2b] is negligible in practice (see Remark 3.12).
In light of the previous result, we consider Theorem 4.6 with the choice of poles C2,ℓ is as in (4.6). Then, Proof. By setting I A = I, I B = −I in the statement of Theorem 4.6 we get (4.5), so that we just need the bound 1] , where the first inequality follows from Lemma 3.9-(ii) and the last equality from Lemma 3.9-(iii) applied with the map T (z). The claim follows combining this inequality ρ [ a, 1] ρ [a,2b] from Lemma 4.7. Example 4.9. We consider the same matrix A of Example 3.19, and we evaluate the matrix function on M = I⊗A+A⊗I, applied to a vector v = vec(F ), where F is a random rank 1 matrix, generated by taking the outer product of two unit vectors with normally distributed entries. We note that the bound from Corollary 4.8 accurately predicts the asymptotic convergence rate, even though it is off by a constant; we believe that this is due to the artificial introduction of (1 + κ) in the Galerkin projection bound, which is usually very pessimistic in practice [3].

4.4.
Low-rank approximability of X. The Kronecker-structured rational Krylov method that we have discussed provides a practical way to compute the evaluation of the matrix function under consideration. However, it can be used also theoretically to predict the decay in the singular values of the computed matrix X, and therefore to describe its approximability properties in a low-rank format.
Proof. We note that the approximation X ℓ obtained using the rational Krylov method with the poles given by Corollary 4.1 has rank (at most) ℓk, and X − X ℓ 2 16γ ℓ,κ f (0)ρ ℓ 2 [a,b] . The claim follows by applying the Eckart-Young theorem.

Cauchy-Stieltjes functions.
In the case of Cauchy-Stieltjes function, the error estimate in Corollary 4.8 would provides a result completely analogue to Theorem 4.10. However, the bound obtained this way involves the multiplicative factor 1 + κ; this can be avoided relying on an alternative strategy.
The idea is to consider the close connection between the rational problem (3.5) and the approximate solution returned by the factored Alternating Direction Implicit method (fADI) [3,7,27]. More specifically, for t 0 let us denote with X t , the solution of the shifted Sylvester equation . . , β ℓ }, provides an approximate solution X ADI ℓ (t) of (4.7) such that its column and row span belong to the spaces Note that the space V ADI ℓ does not depend on t because the right coefficient of (4.7) does not depend on t. If we denote by U ADI ℓ (t) and V ADI ℓ orthonormal bases for these spaces, we have X ADI ) * , and using the ADI error representation [3,27] we obtain In particular, X ADI ℓ (t) is a uniformly good approximation of X t having rank (at most) ℓk and its low-rank factorization has the same right factor ∀t 0.
Proof. Let us define X ℓ : does not depend on t we can take it out from the integral, and therefore X ℓ has rank bounded by ℓk. Then, applying the Eckart-Young theorem we have the inequality

Conclusions, possible extensions and open problems.
We have presented a pole selection strategy for the rational Krylov methods when approximating the action of Laplace-Stieltjes and Cauchy-Stieltjes matrix functions on a vector. The poles have been shown to provide a fast convergence rate and explicit error bounds have been established. Then, the approach presented in [9] that addresses the case of a matrix argument with a Kronecker sum structure has been extended to use rational Krylov subspaces. We have proposed a pole selection strategy that ensures a good exponential rate of convergence of the error norm. From the theoretical perspective we established decay bounds for the singular values of vec −1 (f (I ⊗A−B T ⊗I)vec(F )) when F is low-rank. This generalizes the well known low-rank approximability property of the solutions of Sylvester equations with low-rank right hand side.
There are some research lines that naturally stem from this work. For instance, we have assumed for simplicity to be working with Hermitian positive definite matrices. This assumption might be relaxed, by considering non-normal matrices with field of values included in the positive half plane. Designing an optimal pole selection for such problems would require the solution of Zolotarev problems on more general domains, and deserves further study. In addition, since the projected problem is also non-normal, the fast diagonalization approach for the evaluation proposed in Section 4.1 might not be applicable or stable, and therefore an alternative approach would need to be investigated. Concerning the first term, it is immediate to see that the integrand uniformly converges to the 1/s for ǫ → 0, and therefore the first terms goes to 1 2 for ǫ small. We now focus on the second part.
We can rephrase the ratio of polynomials defining r(s) as follows: Then, we note that the above ratio restricted to the points of the form is yields a complex number of modulus one, that must have the form p The above integral is well-defined even if we let ǫ → 0, we can can take the limit in (A.2) which yields exactly the value 1 2 for the first term, and we have reduced the problem to estimate f (t) = 1 2 + (−1) ℓ π ∞ 0 sin(st + θ(s)) s ds.
To bound the integral, we split the integration domain in three parts: 1 π where we choose ν = a tan( π 4ℓ ) and ξ = 4ℓb. For later use, we note that aπ 4ℓ ν a ℓ . Concerning I 1 , we can further split the integral as Note that |θ(s)| 2s ℓ j=1 γ −1 j , which can be obtained making use of the inequality |atan(t)| |t|. We can bound the second integral term as follows: sin(y) cos(θ(y/t)) y dy.
Note that on [0, tν] the function cos(θ(y/t)) is indeed decreasing, thanks to our choice of ν, and therefore the above can be bounded in modulus by 1 π tν 0 sin(y) cos(θ(y/t)) y dy 1.852 π , where we have used that cos(θ(0)) = 1, and applied Lemma A.