Rational Krylov for Stieltjes matrix functions: convergence and pole selection

TL;DR

This paper develops effective pole selection strategies for rational Krylov methods to efficiently compute matrix functions, especially for Cauchy- and Laplace-Stieltjes functions, with applications to tensor-structured problems.

Contribution

It introduces new pole selection techniques with convergence bounds for rational Krylov methods applied to Stieltjes matrix functions, including tensor-structured cases.

Findings

Effective pole selection strategies with explicit convergence bounds.

Application to tensor-structured matrices in fractional diffusion problems.

Leverages Kronecker structure for efficient computation.

Abstract

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. Relying on the same tools used to analyze the generic situation, we then focus on the case $\mathcal M=I \otimes A - B^T \otimes I$, and $v$ obtained vectorizing a low-rank matrix; this 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.

Tables

Table 1. Table 1: Summary of the convergence rates for rational Krylov methods with the proposed poles. The convergence rate ρ [ α , β ] subscript 𝜌 𝛼 𝛽 \rho_{[\alpha,\beta]} is defined by ρ [ α , β ] := exp ⁡ ( − π 2 / log ⁡ ( 4 ​ β α ) ) assign subscript 𝜌 𝛼 𝛽 superscript 𝜋 2 4 𝛽 𝛼 \rho_{[\alpha,\beta]}:=\exp(-\pi^{2}/\log(4\frac{\beta}{\alpha})) .

Function class	Argument	Error bound	Reference
Laplace-Stieltjes	$ℳ:=A$	${\parallel x-x_{\ell }\parallel }_2\sim 𝒪\left({\rho }_{[a,b]}^{\frac{\ell }{2}}\right)$	Cor. 3.21
	$ℳ:=I\otimes A-B^T\otimes I$	${\parallel X-X_{\ell }\parallel }_2\sim 𝒪\left({\rho }_{[a,b]}^{\frac{\ell }{2}}\right)$	Cor. 4.1
Cauchy-Stieltjes	$ℳ:=A$	${\parallel x-x_{\ell }\parallel }_2\sim 𝒪\left({\rho }_{[a,4b]}^{\ell }\right)$	Cor. 3.24
	$ℳ:=I\otimes A-B^T\otimes I$	${\parallel X-X_{\ell }\parallel }_2\sim 𝒪\left({\rho }_{[a,2b]}^{\ell }\right)$	Cor. 4.12

Table 2. Table 2: Comparison of the time and number of iterations required for computing A − 1 2 ​ v superscript 𝐴 1 2 𝑣 A^{-\frac{1}{2}}v with different relative tolerances using markovfunm , EDS, and extended Krylov. The argument A 𝐴 A is the 100000 × 100000 100000 100000 100000\times 100000 matrix trid ​ ( − 1 , 2 , − 1 ) trid 1 2 1 \mathrm{trid}(-1,2,-1) .

	markovfunm		EDS		extended Krylov
Relative error	Time (s)	Its	Time (s)	Its	Time (s)	Its
${10}^{-1}$	$0.37$	$18$	$0.10$	$7$	$0.32$	$20$
${10}^{-2}$	$0.76$	$29$	$0.26$	$14$	$2.17$	$64$
${10}^{-3}$	$1.17$	$38$	$0.37$	$18$	$4.79$	$106$
${10}^{-4}$	$1.64$	$47$	$0.43$	$20$	$8.16$	$144$
${10}^{-5}$	$2.01$	$53$	$0.58$	$24$	$12.32$	$180$
${10}^{-6}$	$2.56$	$61$	$0.82$	$31$	$16.72$	$212$

Equations

$$f(z)={\int }_0^{\infty }g(t,z)\mu (t)dt,g(t,z)\in \left\{e^{-tz},\frac{1}{t+z}\right\}.$$

$$\mathcal{M}:=I\otimes A-B^T\otimes I\in {\mathbb{C}}^{n^2\times n^2},$$

$$f(z)={\int }_0^{\infty }g(t,z)\mu (t)dt,$$

$$f(z)={\int }_0^{\infty }{e}^{-tz}\mu (t)dt.$$

$$f(z)$$

$${\phi }_j(z):={\int }_0^{\infty }{e}^{-tz}\frac{{\left[\max \{1-t,0\}\right]}^{j-1}}{(j-1)!}dt,j⩾1.$$

$$f(z)={\int }_0^{\infty }\frac{\mu (t)}{t+z}dt.$$

$$f(z)$$

$$f(z)={\int }_0^{\infty }\frac{{\mu }_C(t)}{t+z}dt={\int }_0^{\infty }{\int }_0^{\infty }{e}^{-s(t+z)}{\mu }_C(t)dsdt={\int }_0^{\infty }{e}^{-sz}\underset{{\mu }_L(s)}{\underbrace{{\int }_0^{\infty }{e}^{-st}{\mu }_C(t)dt}}ds,$$

$$f(A)v\approx x_{\mathcal{W}}:=Wf(A_{\mathcal{W}})v_{\mathcal{W}}.$$

$${\mathcal{K}}_{\ell }(A,v):=\mathrm{span}\{v,Av,\dots ,A^{\ell }v\}=\{p(A)v:\mathrm{deg}(p)⩽\ell \}.$$

$${\mathcal{R}\mathcal{K}}_{\ell }(A,v,\Psi )=q_{\ell }(A{)}^{-1}{\mathcal{K}}_{\ell }(A,v)=\mathrm{span}\{q_{\ell }(A{)}^{-1}v,q_{\ell }(A{)}^{-1}Av,\dots ,q_{\ell }(A{)}^{-1}{A}^{\ell }v\},$$

$$\Vert x_{\mathcal{W}}-x{\Vert }_2⩽2\cdot \Vert v{\Vert }_2\cdot \underset{r(z)\in \frac{{\mathcal{P}}_{\ell }}{\Psi }}{\min }\underset{z\in [a,b]}{\max }|f(z)-r(z)|.$$

$$f(A)v={\int }_0^{\infty }g(t,A)\mu (t)dt,g(t,A)\in \{e^{-tA},(tI+A{)}^{-1}\}$$

$$x_{\mathcal{W}}=Wf(A_{\mathcal{W}})v_{\mathcal{W}}={\int }_0^{\infty }Wg(t,A_{\mathcal{W}})v_{\mathcal{W}}\mu (t)dt.$$

$$\Vert (tI+A{)}^{-1}v-U(tI+A_{\ell }{)}^{-1}{v}_{\ell }{\Vert }_2⩽\frac{2}{t+a}\Vert v{\Vert }_2\underset{r(z)\in \frac{{\mathcal{P}}_{\ell }}{\Psi }}{\min }\frac{{\max }_{z\in [a,b]}|r(z)|}{{\min }_{z\in (-\infty ,0]}|r(z)|}$$

$$f_{\mathrm{skel}}(t,\lambda ):=\left[\begin{array}{ccc}\frac{1}{t_1+\lambda }& \dots & \frac{1}{t_{\ell }+\lambda }\end{array}\right]M^{-1}\left[\begin{array}{c}\frac{1}{t+{\lambda }_1}\\ ⋮\\ \frac{1}{t+{\lambda }_{\ell }}\end{array}\right],M_{ij}=\frac{1}{t_j+{\lambda }_i}\in {\mathbb{C}}^{\ell \times \ell },$$

$$\Vert (tI+A{)}^{-1}v-U(tI+A{)}^{-1}{v}_{\ell }{\Vert }_2⩽\frac{2}{t+a}\Vert v{\Vert }_2\frac{{\max }_{z\in [a,b]}|r(z)|}{{\min }_{z\in (-\infty ,0]}|r(z)|}.$$

$$e^{-tA}=\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }{\int }_{-\mathbf{i}T}^{\mathbf{i}T}{e}^{st}(sI+A{)}^{-1}ds.$$

$$\Vert {\mathcal{L}}^{-1}[\hat{r}(s)]{\Vert }_{L^{\infty }({\mathbb{R}}_{+})}⩽{\gamma }_{\ell ,\kappa },{\gamma }_{\ell ,\kappa }:=2.23+\frac{2}{\pi }\log \left(4\ell \cdot \sqrt{\frac{\kappa }{\pi }}\right),$$

$$\Vert e^{-tA}v-U{e}^{-t{A}_{\ell }}{v}_{\ell }{\Vert }_2⩽4{\gamma }_{\ell ,\kappa }\Vert v{\Vert }_2\underset{z\in [a,b]}{\max }|r_{\Psi }(z)|,$$

$$f_{\mathrm{skel}}(s,\lambda ):=\left[\begin{array}{ccc}\frac{1}{\lambda -{\psi }_1}& \dots & \frac{1}{\lambda -{\psi }_{\ell }}\end{array}\right]M^{-1}\left[\begin{array}{c}\frac{1}{s-{\psi }_1}\\ ⋮\\ \frac{1}{s-{\psi }_{\ell }}\end{array}\right],M_{ij}=-\frac{1}{{\psi }_i+{\psi }_j},$$

$$e^{-tA}v-U{e}^{-t{A}_{\ell }}{v}_{\ell }=\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }{\int }_{-\mathbf{i}T}^{\mathbf{i}T}{e}^{st}(sI+A{)}^{-1}v-e^{st}U(sI+A_{\ell }{)}^{-1}{v}_{\ell }ds.$$

$$e^{-tA}v-U{e}^{-t{A}_{\ell }}{v}_{\ell }$$

$$=\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }{\int }_{-\mathbf{i}T}^{\mathbf{i}T}{e}^{st}(sI+A{)}^{-1}\left[I-(sI+A)f_{\mathrm{skel}}(s,A)\right]vds$$

$$-\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }U{\int }_{-\mathbf{i}T}^{\mathbf{i}T}{e}^{st}(sI+A_{\ell }{)}^{-1}\left[I-(sI+A_{\ell })f_{\mathrm{skel}}(s,A_{\ell })\right]v_{\ell }ds$$

$$=\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }{\int }_{-\mathbf{i}T}^{\mathbf{i}T}{e}^{st}(sI+A{)}^{-1}{r}_{\Psi }(A)r_{\Psi }(-s{)}^{-1}vds$$

$$-\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }U{\int }_{-\mathbf{i}T}^{\mathbf{i}T}{e}^{st}(sI+A_{\ell }{)}^{-1}{r}_{\Psi }(A_{\ell })r_{\Psi }(-s{)}^{-1}{v}_{\ell }ds.$$

$$\Vert e^{-tA}v-U{e}^{-t{A}_{\ell }}{v}_{\ell }{\Vert }_2⩽2\underset{\lambda \in [a,b]}{\max }|h(t,\lambda )|,h(t,\lambda ):=\frac{1}{2\pi \mathbf{i}}\underset{T\to \infty }{\lim }{\int }_{-iT}^{iT}{e}^{st}\frac{1}{s+\lambda }\frac{r_{\Psi }(\lambda )}{r_{\Psi }(-s)}ds.$$

$$h(t,\lambda )=r_{\Psi }(\lambda )\cdot {\mathcal{L}}^{-1}\left[\frac{1}{\lambda +s}\frac{p(s)}{p(-s)}\right](t)$$

$$=r_{\Psi }(\lambda )\cdot (\delta (t)-\lambda e^{-\lambda t})\star {\mathcal{L}}^{-1}\left[\frac{1}{s}\frac{p(s)}{p(-s)}\right](t),$$

$$|h(\lambda ,t)|⩽2{\gamma }_{\ell ,\kappa }|r_{\Psi }(\lambda )|.$$

Full text

Relying on the same tools used to analyze the generic situation, we then focus on the case $\mathcal{M}=I\otimes A-B^{T}\otimes I$, and $v$ obtained vectorizing a low-rank matrix; this 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.

1 Introduction

We are concerned with the evaluation of $x=f(\mathcal{M})v$, where $f(z)$ is a Stieltjes function, which can be expressed in integral form

[TABLE]

The two choices for $g(t,z)$ define Laplace-Stieltjes and Cauchy-Stieltjes functions, respectively [8, 31]. 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)}\geqslant 0$ over $\mathbb{R}_{+}$ for any $j\in\mathbb{N}$.

We are interested in two instances of this problem; first, we consider the case $\mathcal{M}:=A$, where $A\in\mathbb{C}^{n\times n}$ is Hermitian positive definite with spectrum contained in $[a,b]$, $v\in\mathbb{C}^{n\times s}$ is a generic (block) vector, and a rational Krylov method [18] is used to approximate $x=f(\mathcal{M})v$. In this case, we want to estimate the Euclidean norm of the error $\lVert x-x_{\ell}\rVert_{2}$, where $x_{\ell}$ is the approximation returned by $\ell$ steps of the method. Second, we consider

[TABLE]

where $A,-B\in\mathbb{C}^{n\times n}$ are Hermitian positive definite with spectra contained in $[a,b]$, $v=\mathrm{vec}(F)\in\mathbb{C}^{n^{2}}$ is the vectorization of a low-rank matrix $F=U_{F}V_{F}^{T}\in\mathbb{C}^{n\times n}$, and a tensorized rational Krylov method [8] is used for computing $\mathrm{vec}(X)=f(\mathcal{M})\mathrm{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 $\lVert X-X_{\ell}\rVert_{2}$, where $X_{\ell}$ is the approximation obtained after $\ell$ steps.

1.1 Main contributions

This paper discusses the connection between rational Krylov evaluation of Stieltjes matrix functions and the parameter dependent rational approximation (with the given poles) of the kernel functions $e^{-tz}$ and $\frac{1}{t+z}$.

The contributions of this work are the following:

Corollary 3.21 provides a choice of poles for the rational Krylov approximation of $f(\mathcal{M})v$, where $f(z)$ is Laplace-Stieltjes, with an explicit error bound depending on the spectrum of $A$. 2. 2.

Similarly, for Cauchy-Stieltjes functions, we show (in Corollary 3.24) how leveraging an approach proposed in [14] allows to recover a result previously given in [4] using different theoretical tools. 3. 3.

In Section 3.5, we obtain new nested sequences of poles by applying the approach of equidistributed sequences to the results in Corollary 3.21–3.24. 4. 4.

In the particular case where $\mathcal{M}:=I\otimes A-B^{T}\otimes I$ we extend the analysis recently proposed in [8] to rational Krylov subspaces. Also in this setting, we provide explicit choices for the poles and explicit convergence bounds. For Laplace-Stieltjes functions a direct consequence of the analysis mentioned above leads to Corollary 4.1; in the Cauchy case, we describe a choice of poles that enables the simultaneous solution of a set of parameter dependent Sylvester equations. This results in a practical choice of poles and an explicit error bound given in Corollary 4.12. 5. 5.

Finally, we give results predicting the decay in the singular values of $X$ where $\mathrm{vec}(X)=f(\mathcal{M})\mathrm{vec}(F)$, $F$ is a low-rank matrix, and $f(z)$ is either Laplace-Stieltjes (Theorem 4.15) or Cauchy-Stieltjes (Theorem 4.17). This generalizes the well-known low-rank approximability of the solutions of of Sylvester equations with low-rank right hand sides [5]. The result for Laplace-Stieltjes follows by the error bound for the rational Krylov method and an Eckart-Young argument. The one for Cauchy-Stieltjes requires to combine the integral representation with the ADI approximant for the solution of matrix equations.

The error bounds obtained are summarized in Table 1.

We recall that completely monotonic functions are well approximated by exponential sums [11]. Another consequence of our results in the Laplace-Stieltjes case is to constructively show that they are also well-approximated by rational functions.

1.2 Motivating problems

Computing the action of a matrix function on a vector is a classical task in numerical analysis, and finds applications in several fields, such as complex networks [7], signal processing [29], numerical solution of ODEs [20], and many others.

Matrices with the Kronecker sum structure as in (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 [28] 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 [23, 25].

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 as the solution of matrix equations with particularly structured coefficients [12, 24]. However, a promising formulation which simplifies the design of boundary conditions relies on first discretizing the 2D Laplacian on the chosen domain, and then considers the action of the matrix function $z^{-\alpha}$ (with the Laplacian as argument) on the right hand side. This is known in the literature as the matrix transform method [32]. In this framework, one has $0<\alpha<1$, and therefore $z^{-\alpha}$ is a Cauchy-Stieltjes function, a property that has been previously exploited for designing fast and stable restarted polynomial Krylov methods for its evaluation [27]. 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 $\varphi_{j}(z)$ functions [20], which are part of the Laplace-Stieltjes class. This has been the subject of deep studies concerning (restarted) polynomial and rational Krylov methods [17, 27]. 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 [30].

The paper is organized as follows. In Section 2 we recall the definitions and main properties of Stieltjes functions. Then, in Section 3 we recall the rational Krylov method and then we analyze the simultaneous approximation of parameter dependent exponentials and resolvents; this leads to the choice of poles and convergence bounds for Stieltjes functions given in Section 3.4. In Section 4 we provide an analysis of the convergence of the method proposed in [8] when rational Krylov subspaces are employed. In particular, in Section 4.4 we provide decay bounds for the singular values of $X$ such that $\mathrm{vec}(X)=f(\mathcal{M})\mathrm{vec}(F)$. We give some concluding remarks and outlook in Section 5.

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:

[TABLE]

where $\mu(t)dt$ is a (non-negative) measure on $[0,\infty]$, and $g(t,z)$ is integrable with respect to that measure. The choice of $g(t,z)$ determines the particular class of Stieltjes functions under consideration (Laplace-Stieltjes or Cauchy-Stieltjes), and $\mu(t)$ is called the density of $f(z)$. $\mu(t)$ can be a proper function, or a distribution, e.g. a Dirac delta. In particular, we can restrict the domain of integration to a subset of $(0,\infty)$ imposing that $\mu(t)=0$ elsewhere. We refer the reader to [31] for further details.

2.1 Laplace-Stieltjes functions

Laplace-Stieltjes functions are obtained by setting $g(t,z)=e^{-tz}$ in (3).

Definition 2.1**.**

Let $f(z)$ be a function defined on $(0,+\infty)$. Then, $f(z)$ is a Laplace-Stieltjes function if there is a positive measure $\mu(t)dt$ on $\mathbb{R}_{+}$ such that

[TABLE]

Examples of Laplace–Stieltjes functions include:

[TABLE]

The last example is an instance of a particularly relevant class of Laplace-Stieltjes functions, with applications to exponential integrators. These are often denoted by $\varphi_{j}(z)$, and can be defined as follows:

[TABLE]

A famous theorem of Bernstein states the equality between the set of Laplace–Stieltjes functions and those of completely monotonic functions in $(0,\infty)$ [10], that is $(-1)^{j}f^{(j)}(z)$ is positive over $(0,\infty)$ for any $j\in\mathbb{N}$.

From the algorithmic point of view, the explicit knowledge of the Laplace density $\mu(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.

2.2 Cauchy-Stieltjes functions

Cauchy-Stieltjes functions form a subclass of Laplace-Stieltjes functions, and are obtained from (3) by setting $g(t,z)=(t+z)^{-1}$.

Definition 2.2**.**

Let $f(z)$ be a function defined on $\mathbb{C}\setminus\mathbb{R}_{-}$. Then, $f(z)$ is a Cauchy-Stieltjes function if there is a positive measure $\mu(t)dt$ on $\mathbb{R}_{+}$ such that

[TABLE]

A few examples of Cauchy-Stieltjes functions are:

[TABLE]

The rational functions with poles on the negative real semi-axis do not belong to this class if one requires $\mu(t)$ to be a function, but they can be obtained by setting $\mu(t)=\sum_{j=1}^{h}\alpha_{j}\delta(t-\beta_{j})$, where $\delta(\cdot)$ is the Dirac delta with unit mass at [math]. For instance, $z^{-1}$ is obtained setting $\mu(t):=\delta(t)$.

Since Cauchy-Stieltjes functions are also completely monotonic in $(0,\infty)$ [9], 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 $\mu_{C}(t)$, one can construct a Laplace-Stieltjes representation as follows:

[TABLE]

where $\mu_{L}(s)$ defines the Laplace-Stieltjes density. In particular, note that if $\mu_{C}(t)$ is positive, so is $\mu_{L}(s)$. For a more detailed analysis of the relation between Cauchy- and Laplace-Stieltjes functions we refer the reader to [31, Section 8.4].

As in the Laplace case, the explicit knowledge of $\mu(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 $\mathcal{W}\subseteq\mathbb{C}^{n}$ is computed, together with the projections $A_{\mathcal{W}}:=W^{*}AW$ and $v_{\mathcal{W}}:=W^{*}v$. Then, the action of $f(A)$ on $v$ is approximated by:

[TABLE]

Intuitively, the choice of the subspace $\mathcal{W}$ is crucial for the quality of the approximation. Usually, one is interested in providing a sequence of subspaces $\mathcal{W}_{1}\subset\mathcal{W}_{2}\subset\mathcal{W}_{3}\subset\dots$ and study the convergence of $x_{\mathcal{W}_{j}}$ to $f(A)v$ as $j$ increases. A common choice for the space $\mathcal{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\times n$ matrix, and $v\in\mathbb{C}^{n\times s}$ be a (block) vector. The Krylov subspace of order $\ell$ generated by $A$ and $v$ is defined as

[TABLE]

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 that $\lVert f(A)v-p(A)v\rVert_{2}\leqslant\max_{z\in[a,b]}|f(z)-p(z)|\cdot\lVert v\rVert_{2}$.

Rational Krylov subspaces are their rational analogue, and can be defined as follows.

Definition 3.2**.**

Let $A$ be a $n\times n$ matrix, $v\in\mathbb{C}^{n\times s}$ be a (block) vector and $\Psi=(\psi_{1},\ldots,\psi_{\ell})$, with $\psi_{j}\in\overline{\mathbb{C}}:=\mathbb{C}\cup\{\infty\}$. The rational Krylov subspace with poles $\Psi$ generated by $A$ and $v$ is defined as

[TABLE]

*where $q_{\ell}(z)=\prod_{j=1}^{\ell}(z-\psi_{j})$ and if $\psi_{j}=\infty$, then we omit the factor $(z-\psi_{j})$ from $q_{\ell}(z)$. *

Note that, a Krylov subspace is a particular rational Krylov subspace where all poles are chosen equal to $\infty$: $\mathcal{RK}_{\ell}(A,v,(\infty,\ldots,\infty))=\mathcal{K}_{\ell}(A,v)$. A common strategy of pole selection consists in alternating [math] and $\infty$. The resulting vector space is known in the literature as the extended Krylov subspace [13].

We denote by $\mathcal{P}_{\ell}$ the set of polynomials of degree at most $\ell$, and by $\mathcal{R}_{\ell,\ell}$ the set of rational functions $g(z)/l(z)$ with $g(z),l(z)\in\mathcal{P}_{\ell}$. Given $\Psi=\{\psi_{1},\ldots,\psi_{\ell}\}\subset\overline{\mathbb{C}}$, we indicate with $\frac{\mathcal{P}_{\ell}}{\Psi}$ the set of rational functions of the form $g(z)/l(z)$, with $g(z)\in\mathcal{P}_{\ell}$ and $l(z):=\prod_{\psi_{j}\in\Psi\setminus\{\infty\}}(z-\psi_{j})$.

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\times n$ matrix, $v\in\mathbb{C}^{n\times s}$ be a (block) vector and $\Psi=\{\psi_{1},\ldots,\psi_{\ell}\}\subset\overline{\mathbb{C}}$. If $U_{\mathcal{P}},U_{\mathcal{R}}$ are orthonormal bases of $\mathcal{K}_{\ell}(A,v)$ and $\mathcal{RK}_{\ell}(A,v,\Psi)$, respectively, then:

$f(z)\in\mathcal{P}_{\ell}\quad\Longrightarrow\quad f(A)v=U_{\mathcal{P}}f(A_{\ell})(U_{\mathcal{P}}^{*}v)\in\mathcal{K}_{\ell}(A,v),\quad A_{\ell}=U^{*}_{\mathcal{P}}AU_{\mathcal{P}}$, 2. 2.

$f(z)\in\frac{\mathcal{P}_{\ell}}{\Psi}\quad\Longrightarrow\quad f(A)v=U_{\mathcal{R}}f(A_{\ell})(U_{\mathcal{R}}^{*}v)\in\mathcal{RK}_{\ell}(A,v,\Psi)$, $\quad A_{\ell}=U^{*}_{\mathcal{R}}AU_{\mathcal{R}}$.

Lemma 3.3 enables to prove the quasi-optimality of the Galerkin projection described in Section 3. Indeed, if $\mathcal{W}:=\mathcal{RK}(A,v,\Psi)$, then [18]

[TABLE]

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-Stieltjes or Cauchy-Stieltjes function.

3.2 Simultaneous approximation of resolvents and matrix exponentials

The integral expression (1) reads as

[TABLE]

when evaluated at a matrix argument. Since the projection is a linear operation we have

[TABLE]

This suggests to look for a space approximating uniformly well, in the parameter $t$, matrix exponentials and resolvents, respectively. A result concerning the approximation error in the $L^{2}$ norm for $t\in\mathbf{i}\mathbb{R}$ is given in [14, Lemma 4.1]. The proof is obtained exploiting some results on the skeleton approximation of $\frac{1}{t+\lambda}$ [26]. We provide a pointwise error bound, which can be obtained by following the same steps of the proof of [14, Lemma 4.1]. We include the proof for completeness.

Theorem 3.4**.**

Let $A$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $U$ be an orthonormal basis of $\mathcal{U_{R}}=\mathcal{RK}_{\ell}(A,v,\Psi)$. Then, $\forall t\in[0,\infty)$, we have the following inequality:

[TABLE]

*where $A_{\ell}=U^{*}AU$ and $v_{\ell}=U^{*}v$. *

Proof 3.5**.**

Following the construction in [26], we consider the function $f_{\mathrm{skel}}(\lambda,t)$ defined by

[TABLE]

where $M_{ij}$ are the entries of $M$ and $(t_{j},\lambda_{i})$ is an $\ell\times\ell$ grid of interpolation nodes. The function $f_{\mathrm{skel}}(t,\lambda)$ is usually called Skeleton approximation and it is practical for approximating $\frac{1}{t+\lambda}$; indeed its relative error takes the explicit form: $1-(t+\lambda)f_{\mathrm{skel}}(t,\lambda)=\frac{r(\lambda)}{r(-t)}$ with $r(z)=\prod_{j=1}^{\ell}\frac{z-\lambda_{j}}{t_{j}+z}$. If this ratio of rational functions is small, then $f_{\mathrm{skel}}(t,\lambda)$ is a good approximation of $\frac{1}{t+\lambda}$ and — consequently — $f_{\mathrm{skel}}(t,A)$ is a good approximation of $(tI+A)^{-1}$. Note that, for every fixed $t$, $f_{\mathrm{skel}}(t,\lambda)$ is a rational function in $\lambda$ with poles $-t_{1},\ldots,-t_{\ell}$. Therefore, using the poles $\psi_{j}=-t_{j}$, $j=1,\ldots,\ell$ for the projection we may write, thanks to (6):

[TABLE]

*Taking the minimum over the possible choices of the parameters $\lambda_{j}$ we get (7). *

Concerning the rational approximation of the (parameter dependent) exponential, the idea is to rely on its Laplace transform that involves the resolvent:

[TABLE]

In this formulation, it is possible to exploit the Skeleton approximation of $\frac{1}{s+\lambda}$ 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 B.

Lemma 3.6**.**

Let $\mathcal{L}^{-1}[\widehat{r}(s)]$ be the inverse Laplace transform of $\widehat{r}(s)=\frac{1}{s}\frac{p(s)}{p(-s)}$, where $p(s)$ is a polynomial of degree $\ell$ with positive real zeros contained in $[a,b]$. Then,

[TABLE]

*where $\kappa=\frac{b}{a}$. *

Theorem 3.7**.**

Let $A$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $U$ be an orthonormal basis of $\mathcal{U_{R}}=\mathcal{RK}_{\ell}(A,v,\Psi)$, where $\Psi=\{\psi_{1},\dots,\psi_{\ell}\}\subseteq[-b,-a]$. Then, $\forall t\in[0,\infty)$, we have the following inequality:

[TABLE]

*where $A_{\ell}=U^{*}AU$, $v_{\ell}=U^{*}v$, $\kappa:=\frac{b}{a}$, $r_{\Psi}(z)\in\mathcal{R}_{\ell,\ell}$ is the rational function defined by $r_{\Psi}(z):=\prod_{j=1}^{\ell}\frac{z+\psi_{j}}{z-\psi_{j}}$ and $\gamma_{\ell,\kappa}$ is the constant defined in Lemma 3.6. *

Proof 3.8**.**

We consider the Skeleton approximation of $\frac{1}{s+\lambda}$ by restricting the choice of poles in both variables to $\Psi$

[TABLE]

where $M_{ij}$ denote the entries of $M$. Then, by using (8) for $A$ and $A_{\ell}$ we get

[TABLE]

Adding and removing the term $e^{st}f_{\mathrm{skel}}(s,A)v=e^{st}Uf_{\mathrm{skel}}(s,A_{\ell})U^{*}v$ inside the integral (the equality holds thanks to Lemma 3.3) we obtain the error expression

[TABLE]

Since $A$ and $A_{\ell}$ are normal, the above integrals can be controlled by the maximum of the corresponding scalar functions on the spectrum of $A$ (and $A_{\ell}$), which yields the bound

[TABLE]

We note that $r_{\Psi}(\lambda)$ can be pulled out of the integral, since it does not depend on $s$, and thus the above can be rewritten as

[TABLE]

where $p(s)$ is as in Lemma 3.6 and $\delta(t)$ indicates the Dirac delta function. Since the $1$-norm of $\delta(t)-\lambda e^{-t\lambda}$ is equal to $2$, using Young’s inequality we can bound $\lVert h(t,\lambda)\rVert_{\infty}\leqslant 2\lVert\mathcal{L}^{-1}\left[\frac{1}{s}\frac{p(s)}{p(-s)}\right]\rVert_{\infty}$. Therefore, we need to estimate the infinity norm of $\mathcal{L}^{-1}\left[\frac{1}{s}\frac{p(s)}{p(-s)}\right](t)$. Such inverse Laplace transform can be uniformly bounded in $t$ by using Lemma 3.6 with a constant that only depends on $\ell$ and $b/a$:

[TABLE]

*This completes the proof. *

Remark 3.9**.**

*The constant provided by Lemma 3.6 is likely not optimal. Indeed, experimentally it seems to hold that $\gamma_{\ell,\kappa}=1$ for any choice of poles in the negative real axis — not necessarily contained in $[-b,-a]$ — and this has been verified in many examples. If this is proved, then the statement of Theorem 3.4 can be made sharper by removing the factor $\gamma_{\ell,\kappa}$. *

3.3 Bounds for the rational approximation problems

Theorem 3.4 and Theorem 3.7 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.10**.**

Let $\Psi\subset\overline{\mathbb{C}}$ be a finite set, and $I_{1},I_{2}$ closed subsets of $\overline{\mathbb{C}}$. Then, we define111We allow the slight abuse of notation of writing $|r(\infty)|$ as the limit of $|r(z)|$ as $|z|\to\infty$, in the case either $I_{1}$ or $I_{2}$ contains the point at infinity.

[TABLE]

The $\theta_{\ell}$ functions enjoy some invariance and inclusion properties, which we report here, and will be extensively used in the rest of the paper.

Lemma 3.11**.**

Let $I_{1},I_{2}$ be subsets of the complex plane, and $\Psi\subset\overline{\mathbb{C}}$. Then, the map $\theta_{\ell}$ satisfies the following properties:

1. (i)

(shift invariance) For any $t\in\mathbb{C}$, it holds $\theta_{\ell}(I_{1}+t,I_{2}+t,\Psi+t)=\theta(I_{1},I_{2},\Psi)$. 2. (ii)

(monotonicity) $\theta_{\ell}(I_{1},I_{2},\Psi)$ is monotonic with respect to the inclusion on the parameters $I_{1}$ and $I_{2}$:

[TABLE] 3. (iii)

(Möbius invariance) If $M(z)$ is a Möbius transform, that is a rational function $M(z)=(\alpha z+\beta)/(\gamma z+\delta)$ with $\alpha\delta\neq\beta\gamma$, then

[TABLE]

Proof 3.12**.**

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)=\frac{p(z)}{h(z)}\in\frac{\mathcal{P}_{\ell}}{\Psi}$ with $M^{-1}(z)$ we obtain another rational function of (at most) the same degree and with poles $M(\Psi)$. Hence, we obtain

[TABLE]

*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 $\Psi$ to vary:

[TABLE]

The latter was posed and studied by Zolotarev in 1877 [33], 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 $I_{1}=-I_{2}=[a,b]$ (10) simplifies to

[TABLE]

which admits the following explicit estimate.

Theorem 3.13** (Zolotarev).**

Let $I=[a,b]$, with $0<a<b$. Then

[TABLE]

In addition, the optimal rational function $r^{[a,b]}_{\ell}(z)$ that realizes the minimum has the form

[TABLE]

*We denote by $\Psi_{\ell}^{[a,b]}:=\{\psi_{1,\ell}^{[a,b]},\dots,\psi_{\ell,\ell}^{[a,b]}\}$ the set of poles of $r^{[a,b]}_{\ell}(z)$. *

Explicit expression for the elements of $\Psi_{\ell}^{[a,b]}$ are available in terms of elliptic functions, see [14, Theorem 4.2].

Remark 3.14**.**

*The original version of Zolotarev’s result involves $\mathrm{exp}(-\frac{\pi^{2}}{\mu(\kappa^{-1})})$ in place of $\rho_{[a,b]}$, where $\mu(\cdot)$ is the Grötzsch ring function. For simplicity, in this paper we prefer the slightly suboptimal form involving the logarithm. We remark that for large $\kappa$ (which is usually the case when considering rational Krylov methods) the difference is negligible [1, Section 17.3]. *

We use Theorem 3.13 and the Möbius invariance property as building blocks for bounding (7). The idea is to map the set $[-\infty,0]\cup[a,b]$ into $[-1,-\widehat{a}]\cup[\widehat{a},1]$ — for some $\widehat{a}\in(0,1)$ — with a Möbius transformation; then make use of Theorem 3.13 and Lemma 3.11-(iii) to provide a convenient choice of $\Psi$ for the original problem.

Lemma 3.15**.**

The Möbius transformation

[TABLE]

maps $[-\infty,0]\cup[a,b]$ into $[-1,-\widehat{a}]\cup[\widehat{a},1]$, with $\widehat{a}:=\frac{\Delta+a-b}{\Delta-a+b}=\frac{b-\Delta}{\Delta+b}$. The inverse map $T_{C}(z)^{-1}$ is given by:

[TABLE]

*Moreover, for any $0<a<b$ it holds $\widehat{a}^{-1}\leqslant\frac{4b}{a}$, and therefore $\rho_{[\widehat{a},1]}\leqslant\rho_{[a,4b]}$. *

Proof 3.16**.**

By direct substitution, we have $T_{C}(-\infty)=-1$, and $T_{C}(b)=1$; moreover, again by direct computation one verifies that $T_{C}(0)+T_{C}(a)=0$, which implies that $T_{C}([-\infty,0])=[-1,-\widehat{a}]$ and $T_{C}([a,b])=[\widehat{a},1]$. Then, we have

[TABLE]

*Using the relation $\sqrt{1-t}\leqslant 1-\frac{t}{2}$ for any $0\leqslant t\leqslant 1$, we obtain that $\widehat{a}^{-1}\leqslant\frac{2b-\frac{a}{2}}{\frac{a}{2}}\leqslant 4\frac{b}{a},$ which concludes the proof. *

Remark 3.17**.**

*We note that the estimate $\rho_{[\widehat{a},1]}\leqslant\rho_{[a,4b]}$ is asymptotically tight, that is the limit of $\rho_{[\widehat{a},1]}/\rho_{[a,4b]}\to 1$ as $b/a\to\infty$. For instance, if $b/a=10$ then the relative error between the two quantities is about $2\cdot 10^{-2}$, and for $b/a=1000$ about $5\cdot 10^{-5}$. Since the interest for this approach is in dealing with matrices that are not well-conditioned, we consider the error negligible in practice. *

In light of Theorem 3.13 and Lemma 3.15, we consider the choice

[TABLE]

in Theorem 3.4. This yields the following estimate.

Corollary 3.18**.**

Let $A$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $U$ be an orthonormal basis of $\mathcal{U_{R}}=\mathcal{RK}_{\ell}(A,v,\Psi_{C,\ell}^{[a,b]})$. Then, $\forall t\in[0,\infty)$

[TABLE]

*where $A_{\ell}=U^{*}AU$ and $v_{\ell}=U^{*}v$. *

When considering Laplace-Stieltjes functions, we may choose as poles $\Psi_{\ell}^{[a,b]}$ which are the optimal Zolotarev poles on the interval $[a,b]$. This enables to prove the following result, which builds on Theorem 3.7.

Corollary 3.19**.**

Let $A$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $U$ be an orthonormal basis of $\ \mathcal{U_{R}}=\mathcal{RK}_{\ell}(A,v,\Psi_{\ell}^{[a,b]})$. Then, $\forall t\in[0,\infty)$

[TABLE]

*where $A_{\ell}=U^{*}AU$ and $v_{\ell}=U^{*}v$. *

Proof 3.20**.**

*The proof relies on the fact that the optimal Zolotarev function evaluated on the interval $[a,b]$ can be bounded by $2\rho_{[a,b]}^{\frac{\ell}{2}}$ [5, Theorem 3.3]. Since its zeros and poles are symmetric with respect to the imaginary axis and real, we can apply Theorem 3.7 to obtain (13). *

3.4 Convergence bounds for Stieltjes functions

Let us consider $f(z)$ a Stieltjes function of the general form (1). Then the error of the rational Krylov method for approximating $f(A)v$ is given by

[TABLE]

where $g(t,A)$ is either a parameter dependent exponential or resolvent function. Therefore Corollary 3.18 and Corollary 3.19 provide all the ingredients to study the error of the rational Krylov projection, when the suggested pole selection strategy is adopted.

Corollary 3.21**.**

Let $f(z)$ be a Laplace-Stieltjes function, $A$ be Hermitian positive definite with spectrum contained in $[a,b]$, $U$ be an orthonormal basis of $\mathcal{U_{R}}=\mathcal{RK}_{\ell}(A,v,\Psi_{\ell}^{[a,b]})$ and $x_{\ell}=Uf(A_{\ell})v_{\ell}$ with $A_{\ell}=U^{*}AU$ and $v_{\ell}=U^{*}v$. Then

[TABLE]

*where $\gamma_{\ell,\kappa}$ is defined as in Theorem 3.7, and $f(0^{+}):=\lim_{z\to 0^{+}}f(z)$. *

Proof 3.22**.**

Since $f(z)$ is a Laplace-Stieltjes function, we can express the error as follows:

[TABLE]

*where we applied (6) and Corollary 3.19 to obtain the second inequality. *

Remark 3.23**.**

*In order to be meaningful, Corollary 3.21 requires the function $f(z)$ to be finite over $[0,\infty)$, which might not be the case in general (consider for instance $x^{-\alpha}$, which is both Cauchy and Laplace-Stieltjes). Nevertheless, the result can be applied to $f(z+\eta)$, which is always completely monotonic for a positive $\eta$, by taking $0<\eta<a$. A value of $\eta$ closer to $a$ gives a slower decay rate, but a smaller constant $f(0^{+})$. Similarly, if $f(z)$ happens to be completely monotonic on an interval larger than $[0,\infty)$, then bounds with a faster asymptotic convergence rate (but a larger constant) can be obtained considering $\eta<0$. *

Corollary 3.18 allows to state the corresponding bound for Cauchy-Stieltjes functions. The proof is analogous to the one of Corollary 3.21.

Corollary 3.24**.**

Let $f(z)$ be a Cauchy-Stieltjes function, $A$ be Hermitian positive definite with spectrum contained in $[a,b]$, $U$ be an orthonormal basis of $\mathcal{U_{R}}=\mathcal{RK}_{\ell}(A,v,\Psi_{C,\ell}^{[a,b]})$ with $\Psi_{C,\ell}^{[a,b]}$ as in (11) and $x_{\ell}=Uf(A_{\ell})v_{\ell}$ with $A_{\ell}=U^{*}AU$ and $v_{\ell}=U^{*}v$. Then

[TABLE]

3.5 Nested sequences of poles

From the computational perspective, it is more convenient to have a nested sequence of subspaces $\mathcal{W}_{1}\subseteq\ldots\mathcal{W}_{j}\subseteq\mathcal{W}_{j+1}\subseteq\ldots$, so that $\ell$ can be chosen adaptively. For example, in [19] the authors propose a greedy algorithm for the selection of the poles taylored to the evaluation of Cauchy-Stieltjes matrix functions. See [15, 16] for greedy pole selection strategies to be applied in different — although closely related — contexts.

The choices of poles proposed in the previous sections require to a priori determine the degree $\ell$ of the approximant $x_{\ell}$. Given a target accuracy, one can use the convergence bounds in Corollary 3.21–3.24 to determine $\ell$. Unfortunately, this is likely to overestimate the minimum value of $\ell$ that provides the sought accuracy.

An option, that allows to overcome this limitation, is to rely on the method of equidistributed sequences (EDS), as described in [14, Section 4]. The latter exploits the limit — as $\ell\to\infty$ — of the measures generated by the set of points $\Psi_{\ell}^{[a,b]},\Psi_{C,\ell}^{[a,b]}$ to return infinite sequences of poles that are guaranteed to provide the same asymptotic rate of convergence. More specifically, the EDS $\{\widetilde{\sigma}_{j}\}_{j\in\mathbb{N}}$ associated with $\Psi_{\ell}^{[a,1]}$ is obtained with the following steps:

1. (i)

Select $\zeta\in\mathbb{R}^{+}\setminus\mathbb{Q}$ and generate the sequence $\{s_{j}\}_{j\in\mathbb{N}}:=\{0,\ \zeta-\lfloor\zeta\rfloor,\ 2\zeta-\lfloor 2\zeta\rfloor,\ 3\zeta-\lfloor 3\zeta\rfloor,\ \dots\}$, where $\lfloor\cdot\rfloor$ indicates the greatest integer less than or equal to the argument; this sequence has as asymptotic distribution (in the sense of EDS) the Lebesgue measure on $[0,1]$. 2. (ii)

Compute the sequence $\{t_{j}\}_{j\in\mathbb{N}}$ such that $g(t_{j})=s_{j}$ where

[TABLE] 3. (iii)

Define $\widetilde{\sigma}_{j}:=\sqrt{t_{j}}$.

More generally, the EDS associated with $\Psi_{\ell}^{[a,b]},\Psi_{C,\ell}^{[a,b]}$ are obtained by applying either a scaling or the Möbius transformation (11) to the EDS for $\Psi_{\ell}^{[a,1]}$.

In our implementation, only the finite portion $\{\widetilde{\sigma}_{j}\}_{j=0,\dots,\ell-1}$ is — incrementally — generated for computing $x_{\ell}$. As starting irrational number we select $\zeta=\frac{1}{\sqrt{2}}$ and each quantity $t_{j}$ is approximated by applying the Newton’s method to the equation $g(t_{j})-s_{j}=0$. The initialization of the Newton iteration is done by approximating $\hat{t}\mapsto g(e^{\hat{t}})-s_{j}$ with a linear function on the domain of interest, and then using the exponential of its only root as starting point. This is done beforehand selecting $t=a^{2}$ and $t=a$ as interpolation points; in our experience, with such starting point Newton’s method converges in a few steps.

3.6 Some numerical tests

3.6.1 Laplace-Stieltjes functions

Let us consider the 1D diffusion problem over $[0,1]$ with zero Dirichlet boundary conditions

[TABLE]

discretized using central finite differences in space with $50000$ points, and integrated by means of the exponential Euler method with time step $\Delta t=0.1$. This requires to evaluate the action of the Laplace-Stieltjes matrix function $\varphi_{1}(\frac{\epsilon}{h^{2}}\Delta tA)v$, where $A$ is the tridiagonal matrix $A=\mathrm{tridiag}(-1,2,-1)$. We test the convergence rates of various choices of poles by measuring the absolute error when using a random vector $v$. Figure 1-(left) reports the results associated with: the poles from Corollary 3.19, the corresponding EDS computed as described in Section 3.5 and the extended Krylov method. It is visible that the three approximations have the same convergence rate, although the choice of poles from Corollary 3.19 and the EDS performs slightly better. We mention that, since $A$ is ill-conditioned, polynomial Krylov performs poorly on this example.

We keep the same settings and we test the convergence rates for the Laplace-Stieltjes function $z^{-\frac{3}{2}}W(z)$ where $W(z)$ is the Lambert W function [22]. The plot in Figure 1-(right) shows that after about $10$ iterations the convergence rate of the extended Krylov method deteriorates, while the poles from Corollary 3.19 and the EDS provide the best convergence rate.

3.6.2 Inverse square root

Let us test the pole selection strategies for Cauchy-Stieltjes functions, by considering the evaluation of $f(z)=z^{-\frac{1}{2}}$ on the $n\times n$ matrix $\mathrm{tridiag}(-1,2,-1)$, for $n=10^{4},5\cdot 10^{4},10^{5}$. The list of methods that we consider includes: the poles $\Psi_{C,\ell}^{[a,b]}$ from Corollary 3.18, the associated EDS, the extended Krylov method and the adaptive strategy proposed in [19] for Cauchy-Stieltjes functions. The latter is implemented in the markovfunmv package available at http://guettel.com/markovfunmv/ which we used for producing the results reported in Figure 2. The poles from Corollary 3.18 and the extended Krylov method yield the best and the worst convergence history, respectively, for all values of $n$. The EDS and markovfunm perform similarly for $n=10^{4}$, but as $n$ increases the decay rate of markovfunm deteriorates significantly.

We consider a second numerical experiment which keeps the same settings apart from the size of the matrix argument which is fixed to $n=10^{5}$. Then, we measure the number of iterations and the computational time needed by the methods using nested sequences of poles, i.e. EK, EDS, markovfunm, to reach different target values for the relative error $\frac{\lVert x-x_{\ell}\rVert_{2}}{\lVert x\rVert_{2}}$. The EK method has the cheapest iteration cost because it exploits the pre-computation of the Cholesky factor of the matrix $A$ for the computation of the orthogonal basis. However, as testified by the results in Table 2, the high number of iterations makes EK competitive only for the high relative error $10^{-1}$. The iteration costs of EDS and markovfunm is essentially the same since they only differ in the computation of the poles, which requires a negligible portion of the computational time. Therefore, the comparison between EDS and markovfunm goes along with the number of iterations which makes the former more efficient222To make a fair comparison between the methods, for this test we relied on the rational Arnoldi implementation found in markovfunm for the implementation of Algorithm 1 using EDS poles.. We remark that in the situation where precomputing the Cholesky gives a larger computational benefit, and memory is not an issue, EK may be competitive again.

We conclude the numerical experiments on the inverse square root by considering matrix arguments with different distributions of the eigenvalues. More precisely, we set $A$ as the diagonal matrix of dimension $n=5\cdot 10^{4}$ with the following spectrum configurations:

• $(i)$

Equispaced values in the interval $[\frac{1}{n},1]$,

• $(ii)$

Eigenvalues of $\mathrm{trid}(-1,2+10^{-3},-1)$ (shifted Laplacian),

• $(iii)$

$20$ Chebyshev points in $[10^{-3},10^{-1}]$ and $n-20$ Chebyshev points in $[10,10^{3}]$.

The convergence histories of the different projection spaces are reported in Figure 3. For all the eigenvalues configurations, EDS and markovfunm provide comparable performances. The poles from Corollary 3.18 performs as EDS and markovfunm on $(ii)$ and slightly better on $(i)$ and $(iii)$. Once again, EK is the one providing the slowest convergence rate on all examples.

3.6.3 Other Cauchy-Stieltjes functions

Finally, we test the convergence rate of the different pole selection strategies for the Cauchy-Stieltjes functions $\frac{1-e^{-\sqrt{z}}}{z},z^{-0.2},z^{-0.8}$ and the matrix argument $A=\mathrm{trid}(-1,2,-1)$.

The results reported in Figure 4 show that in all cases the poles from Corollary 3.18 and the extended Krylov method provide the best and the worst convergence rates, respectively. The EDS converges faster than markovfunm apart from the case of $z^{-0.2}$ where the two strategies perform similarly.

4 Evaluating Stieltjes functions of matrices with Kronecker structure

We consider the task of computing $f(\mathcal{M})v$ where $\mathcal{M}=I\otimes A-B^{T}\otimes I$. This problem often stems from the discretizations of 2D differential equations, such as the matrix transfer method used for fractional diffusion equations [32].

We assume that $v=\mathrm{vec}(F)$, where $F=U_{F}V_{F}^{T}$ where $U_{F}$ and $V_{F}$ are tall and skinny matrices. For instance, when $f(z)=z^{-1}$, this is equivalent to solving the matrix equation $AX-XB=F$. It is well-known that, if the spectra of $A$ and $B$ are separated, then the low-rank property is numerically inherited by $X$ [5]. For more general functions than $z^{-1}$, a projection scheme that preserves the Kronecker structure has been proposed in [8] using polynomial Krylov methods. We briefly review it in Section 4.1. The method proposed in [8] 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\times n$ matrices, applying the projection scheme described in Section 3 requires to build an orthonormal basis $W$ for a (low-dimensional) subspace $\mathcal{W}\subseteq\mathbb{C}^{n^{2}}$, together with the projections of $W^{*}\mathcal{M}W=H$ and $v_{\mathcal{W}}=W^{*}v$. Then the action of $f(\mathcal{M})$ on $v$ is approximated by:

[TABLE]

The trick at the core of the projection scheme proposed in [8] consists in choosing a tensorized subspace of the form $\mathcal{W}:=\mathcal{U}\otimes\mathcal{V}$, spanned by an orthonormal basis of the form $W=U\otimes V$, where $U$ and $V$ are orthonormal bases of $\mathcal{U}\subseteq\mathbb{C}^{n}$ and $\mathcal{V}\subseteq\mathbb{C}^{n}$, respectively. With this choice, the projection of $\mathcal{M}$ onto $\mathcal{U}\otimes\mathcal{V}$ retains the same structure, that is

[TABLE]

where $A_{\mathcal{U}}=U^{*}AU$ and $B_{\mathcal{V}}=V^{*}BV$.

Since in our case $v=\mathrm{vec}(F)$ and $F=U_{F}V_{F}^{T}$, this enables to exploit the low-rank structure as well. Indeed, the projection of $F$ onto $\mathcal{U}\otimes\mathcal{V}$ can be written as $v_{\mathcal{W}}=\mathrm{vec}(F_{\mathcal{W}})=\mathrm{vec}((U^{*}U_{F})(V_{F}^{T}V))$. The high-level structure of the procedure is sketched in Algorithm 1.

At the core of Algorithm 1 is the evaluation of the matrix function on the projected matrix $I\otimes A_{\mathcal{U}}-B_{\mathcal{V}}^{T}\otimes I$. Even when $U,V$ have a low dimension $k\ll n$, this matrix is $k^{2}\times k^{2}$, so it is undesirable to build it explicitly and then evaluate $f(\cdot)$ 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\otimes A_{\mathcal{U}}-B_{\mathcal{V}}^{T}\otimes I$. For a more general function, we can still design a $\mathcal{O}(k^{3})$ procedure for the evaluation of $f(\cdot)$ in our case. Indeed, since $A_{\mathcal{U}}$ and $B_{\mathcal{V}}$ are Hermitian, we may diagonalize them using a unitary transformation as follows:

[TABLE]

Then, the evaluation of the matrix function $f(z)$ with argument $I\otimes A_{\mathcal{U}}-B_{\mathcal{V}}^{T}\otimes I$ can be recast to a scalar problem by setting

[TABLE]

where $\mathcal{D}:=I\otimes D_{A}-D_{B}\otimes I$. If we denote by $X=\mathrm{vec}^{-1}(f(\mathcal{M})\mathrm{vec}(F))$ and with $D$ the matrix defined by $D_{ij}=(D_{A})_{ii}-(D_{B})_{jj}$, then

[TABLE]

where $\circ$ denotes the Hadamard product and $f^{\circ}(\cdot)$ the function $f(\cdot)$ applied component-wise to the entries of $D$: $[f^{\circ}(D)]_{ij}=f(D_{ij})$.

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 $\mathcal{O}(c_{f}\cdot k^{2}+k^{3})$, where $c_{f}$ denotes the cost of a single function evaluation. The procedure is described in Algorithm 2.

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}\otimes I$ and $I\otimes A$ commute, we have

[TABLE]

Consider projecting the matrix $\mathcal{M}$ onto a tensorized subspace spanned by the Kronecker products of unitary matrices $U\otimes V$. This, combined with Algorithm 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.

Corollary 4.1**.**

Let $f(z)$ be a Laplace-Stieltjes function, $A,-B$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $X_{\ell}$ be the approximation of $X=\mathrm{vec}^{-1}(f(\mathcal{M})\mathrm{vec}(F))$ obtained using Algorithm 1 with $U\otimes V$ orthonormal basis of $\mathcal{U_{R}}\otimes\mathcal{V_{R}}=\mathcal{RK}_{\ell}(A,U_{F},\Psi_{\ell}^{[a,b]})\otimes\mathcal{RK}_{\ell}(B^{T},V_{F},\Psi_{\ell}^{[a,b]})$. Then,

[TABLE]

Proof 4.2**.**

If $f(z)$ is a Laplace-Stieltjes function, we may express the error matrix $X-X_{\ell}$ as follows:

[TABLE]

where $A_{\ell}=U^{*}AU$ and $B_{\ell}=V^{*}BV$. Adding and subtracting the quantity $Ue^{-tA_{\ell}}U^{*}Fe^{tB}$ yields the following inequalities:

[TABLE]

*where in the last step we used Corollary 3.19 for both addends. *

Example 4.3**.**

To test the proposed projection spaces we consider the same matrix $A$ of Example 3.6.1, and we evaluate the function $\varphi_{1}$ to $\mathcal{M}=I\otimes A+A\otimes I$, applied to a vector $v=\mathrm{vec}(F)$, where $F$ is a random rank $1$ matrix, generated by taking the outer product of two unit vectors with normally distributed entries. The results are reported in Figure 5.

4.3 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 $\mathcal{M}=I\otimes A-B^{T}\otimes I$ which is Kronecker structured, this acquires a particular meaning.

In fact, relying on the integral representation (3) of $f(z)$ we obtain:

[TABLE]

where $X_{t}:=\mathrm{vec}^{-1}((tI+\mathcal{M})^{-1}\mathrm{vec}(F))$ solves the matrix equation

[TABLE]

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 $\mathcal{K}_{m}(tI+A,U_{F})\otimes\mathcal{K}_{m}(B^{T},V_{F})$ for solving (16) at a certain point $t$. Thanks to the shift invariance of polynomial Krylov subspaces, all these subspaces coincide with $\mathcal{U}_{P}\otimes\mathcal{V}_{P}=\mathcal{K}_{m}(A,U_{F})\otimes\mathcal{K}_{m}(B^{T},V_{F})$. This observation is at the core of the strategy proposed in [8], which makes use of $\mathcal{U}_{P}\otimes\mathcal{V}_{P}$ in Algorithm 1. This allows to use the convergence theory for linear matrix equations to provide error bounds in the Cauchy-Stieltjes case, see [8, 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 $\mathcal{U}_{P}\otimes\mathcal{V}_{P}$. However, they are not shift invariant, and this makes the analysis not straightforward. Throughout this section, we denote by $U\otimes V$ the orthonormal basis of the tensorized rational Krylov subspace

[TABLE]

where $\Psi:=\{\psi_{1},\dots,\psi_{\ell}\}$ and $\Xi:=\{\xi_{1},\dots,\xi_{\ell}\}$ are the prescribed poles. We define the following polynomials of degree (at most) $\ell$:

[TABLE]

and we denote by $A_{\ell}=U^{*}AU$, $B_{\ell}=V^{*}BV$ the projected $(\ell k\times\ell 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 1 applied to $f(z)=z^{-1}$ coincides with the Galerkin projection method for Sylvester equations [28], 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. The proof is reported in the Appendix A.

Theorem 4.4**.**

Let $A,-B$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $X_{\ell}$ be the approximate solution returned by Algorithm 1 using $f(z)=z^{-1}$ and the orthonormal basis $U\otimes V$ of $\mathcal{U}_{R}\otimes\mathcal{V}_{R}=\mathcal{RK}_{\ell}(A,U_{F},\Psi)\otimes\mathcal{RK}_{\ell}(B^{T},V_{F},\Xi)$, then

[TABLE]

Remark 4.5**.**

Using the mixed norm inequality $\lVert AB\rVert_{F}\leqslant\lVert A\rVert_{F}\lVert B\rVert_{2}$, one can state the bound in the Frobenius norm as well:

[TABLE]

*which is tighter than the one in [3]. *

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 $\lVert X-X_{\ell}\rVert_{2}$ causes the appearances of an additional condition number.

Corollary 4.6**.**

If $X_{\ell}$ is the approximate solution of the linear matrix equation $AX-XB=F$ returned by Algorithm 1 as in Theorem 4.4, then

[TABLE]

Proof 4.7**.**

*We note that $X_{\ell}-X$ solves the Sylvester equation $A(X_{\ell}-X)-(X_{\ell}-X)B=R,$ where $R:=AX_{\ell}-X_{\ell}B-F$ verifies $\lVert R\rVert_{2}\leqslant\left(1+\kappa\right)\max\{\theta_{\ell}(I_{A},I_{B},\Psi),\ \theta_{\ell}(I_{B},I_{A},\Xi)\}\lVert F\rVert_{2}$, thanks to Theorem 4.4. In view of [21, Theorem 2.1] $\lVert X_{\ell}-X\rVert_{2}$ is bounded by $\frac{1}{2a}\lVert R\rVert_{2}$. *

4.3.2 Error analysis for Cauchy-Stieltjes functions

In view of Equation 16, the evaluation of Cauchy-Stieltjes function is closely related to solving (in a uniformly accurate way) parameter-dependent Sylvester equations. This connection is clarified by the following result.

Theorem 4.8**.**

Let $f(z)$ be a Cauchy-Stieltjes function, $A,-B$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $X_{\ell}$ be the approximate evaluation of $f(z)$ returned by Algorithm 1 using the orthonormal basis $U\otimes V$ of the subspace $\mathcal{U}_{R}\otimes\mathcal{V}_{R}=\mathcal{RK}_{\ell}(A,U_{F},\Psi)\otimes\mathcal{RK}_{\ell}(B^{T},V_{F},\Xi)$. Then,

[TABLE]

*where $\kappa=\frac{b}{a}$ and $\theta_{\ell}(\cdot,\cdot,\cdot)$ is as in Definition 3.10. *

Proof 4.9**.**

Applying the definition of $f(\mathcal{M})$ we have $f(\mathcal{M})\mathrm{vec}{(F)}=\int_{0}^{\infty}(tI+\mathcal{M})^{-1}\mathrm{vec}{(F)}\mu(t)\ dt.$ We note that, for any $t\geqslant 0$, the vector $\mathrm{vec}(X_{t}):=(tI+\mathcal{M})^{-1}\mathrm{vec}(F)$ is such that $X_{t}$ solves the Sylvester equation $(tI+A)X_{t}-X_{t}B=F$. Then, we can write $X$ as $X=\int_{0}^{\infty}X_{t}\mu(t)\ dt.$

Let us consider the approximation $UY_{t}V^{*}$ to $X_{t}$ obtained by solving the projected Sylvester equation $(tI+U^{*}AU)Y_{t}-Y_{t}(V^{*}BV)=U^{*}FV$, and $Y=\int_{0}^{\infty}Y_{t}\mu(t)\ dt$. We remark that $\mathcal{RK}_{\ell}(A,U_{F},\Psi)=\mathcal{RK}_{\ell}(tI+A,U_{F},\Psi+t)$.

Then, relying on Corollary 4.6, we can bound the error $R_{t}:=\lVert X_{t}-UY_{t}V^{*}\rVert_{2}$ with

[TABLE]

where $C(t):=\frac{2(t+a+b)}{(t+2a)^{2}}$. Making use of Lemma 3.11-(i) we get:

[TABLE]

An estimate for the error on $X$ is obtained by integrating $R_{t}$:

[TABLE]

*where we used that the function $\frac{2(t+a+b)}{t+2a}$ is maximum over $[0,\infty]$ at $t=0$. *

Inspired by Theorem 4.4, we look at the construction of rational functions that make the quantities $\theta_{\ell}(I_{A},I_{B}-t,\Psi)$ and $\theta_{\ell}(I_{B},I_{A}+t,\Xi)$ small. If we choose $\Xi=-\Psi$ then (19) simplifies to

[TABLE]

because $\theta_{\ell}(I_{A},I_{B}-t,\Psi)=\theta_{\ell}(I_{A},-I_{A}-t,\Psi)=\theta_{\ell}(-I_{A},I_{A}+t,-\Psi)=\theta_{\ell}(I_{B},I_{A}+t,-\Psi)$, in view of Lemma 3.11-(iii).

Similarly to the analysis done for Cauchy-Stieltjes function for a generic matrix $A$, we may consider a Möbius transform that maps the Zolotarev problem involving the point at infinity in a more familiar form. More precisely, we aim at mapping the set $[-\infty,-a]\cup[a,b]$ into $[-1,-\widetilde{a}]\cup[\widetilde{a},1]$ — for some $\widetilde{a}\in(0,1)$. Then, we make use of Theorem 3.13 and Lemma 3.11-(iii) to provide a choice of $\Psi$ that makes the quantity $\theta_{\ell}(I_{A},I_{B}-t,\Psi)$ small, independently of $t$.

Lemma 4.10**.**

The Möbius transformation

[TABLE]

maps $[-\infty,-a]\cup[a,b]$ into $[-1,-\widetilde{a}]\cup[\widetilde{a},1]$, with $\widetilde{a}:=\frac{\Delta+a-b}{\Delta-a+b}$. The inverse map $T(z)^{-1}$ is:

[TABLE]

*In addition, we have $\widetilde{a}^{-1}\leqslant 2b/a$, and therefore $\rho_{[\widetilde{a},1]}\leqslant\rho_{[a,2b]}$. *

Proof 4.11**.**

*The proof can be easily obtained following the same steps of Lemma 3.15. As in that case, the overestimate introduced by the inequality $\rho_{[\widetilde{a},1]}\leqslant\rho_{[a,2b]}$ is negligible in practice (see Remark 3.17). *

In light of the previous result, we consider Theorem 4.8 with the choice of poles

[TABLE]

where $\Psi_{\ell}^{[\widetilde{a},1]}$ indicates the set of optimal poles and zeros — provided by Theorem 3.13 — for the domain $[-1,\widetilde{a}]\cup[\widetilde{a},1]$. This yields the following.

Corollary 4.12**.**

Let $f(z)$ be a Cauchy-Stieltjes function with density $\mu(t)$, $A,-B$ be Hermitian positive definite with spectrum contained in $[a,b]$ and $X_{\ell}$ the approximate evaluation of $f(z)$ returned by Algorithm 1 using the orthonormal basis $U\otimes V$ of the subspace $\mathcal{RK}_{\ell}(A,U_{F},\Psi_{C_{2},\ell}^{[a,b]}))\otimes\mathcal{RK}_{\ell}(B^{T},V_{F},-\Psi_{C_{2},\ell}^{[a,b]}))$, where $\Psi_{C_{2},\ell}^{[a,b]}$ is as in (21). Then,

[TABLE]

Proof 4.13**.**

By setting $I_{A}=I$, $I_{B}=-I$ in the statement of Theorem 4.8 we get (20), so that we just need the bound

[TABLE]

*where the first inequality follows from Lemma 3.11-(ii) and the last equality from Lemma 3.11-(iii) applied with the map $T(z)$. The claim follows combining this inequality $\rho_{[\widetilde{a},1]}\leqslant\rho_{[a,2b]}$ from Lemma 4.10. *

Example 4.14**.**

We consider the same matrix $A$ of Example 3.6.2, and we evaluate the inverse square root of $\mathcal{M}=I\otimes A+A\otimes I$, applied to a vector $v=\mathrm{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, in Figure 6, the bound from Corollary 4.12 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+\kappa)$ 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.

4.4.1 Laplace-Stieltjes functions

In the Laplace-Stieltjes case, we may employ Corollary 4.1 directly to provide an estimate for the decay in the singular values.

Theorem 4.15**.**

Let $f(z)$ be a Laplace-Stieltjes function and $\mathcal{M}=I\otimes A-B^{T}\otimes I$ where $A,-B$ are Hermitian positive definite with spectra contained in $[a,b]$. Then, if $\mathrm{vec}(X)=f(\mathcal{M})\mathrm{vec}(F)$, with $F=U_{F}V_{F}^{T}$ of rank $k$, we have

[TABLE]

Proof 4.16**.**

*We note that the approximation $X_{\ell}$ obtained using the rational Krylov method with the poles given by Corollary 4.1 has rank (at most) $\ell k$, and $\lVert X-X_{\ell}\rVert_{2}\leqslant 16\gamma_{\ell,\kappa}f(0^{+})\rho_{[a,b]}^{\frac{\ell}{2}}$. The claim follows by applying the Eckart-Young theorem. *

4.4.2 Cauchy-Stieltjes functions

In the case of Cauchy-Stieltjes function, the error estimate in Corollary 4.12 would provides a result completely analogue to Theorem 4.15. However, the bound obtained this way involves the multiplicative factor $1+\kappa$; this can be avoided relying on an alternative strategy.

The idea is to consider the close connection between the rational problem (10) and the approximate solution returned by the factored Alternating Direction Implicit method (fADI) [3, 5, 6]. More specifically, for $t\geqslant 0$ let us denote with $X_{t}$, the solution of the shifted Sylvester equation

[TABLE]

In view of (16), $X_{t}$ is such that $X=\int_{0}^{\infty}X_{t}\mu(t)dt$. Running fADI for $\ell$ iterations, with shift parameters $T^{-1}(\Psi_{\ell}^{[\widetilde{a},1]})=\{\alpha_{1},\dots,\alpha_{\ell}\}$ and $T^{-1}(-\Psi_{\ell}^{[\widetilde{a},1]})=\{\beta_{1},\dots,\beta_{\ell}\}$, provides an approximate solution $X^{ADI}_{\ell}(t)$ of (22) such that its column and row span belong to the spaces

[TABLE]

Note that the space $V_{\ell}^{\text{ADI}}$ does not depend on $t$ because the right coefficient of (22) does not depend on $t$. If we denote by $U_{\ell}^{\text{ADI}}(t)$ and $V_{\ell}^{\text{ADI}}$ orthonormal bases for these spaces, we have $X^{\text{ADI}}_{\ell}(t)=U_{\ell}^{\text{ADI}}(t)Y_{\ell}^{\text{ADI}}(t)(V_{\ell}^{\text{ADI}})^{*}$, and using the ADI error representation [3, 5] we obtain $\lVert X_{t}-X^{\text{ADI}}_{\ell}(t)\rVert_{2}\leqslant\lVert X_{t}\rVert_{2}\rho_{[a,2b]}^{\ell}.$

In particular, $X^{ADI}_{\ell}(t)$ is a uniformly good approximation of $X_{t}$ having rank (at most) $\ell k$ and its low-rank factorization has the same right factor $\forall t\geqslant 0$.

Theorem 4.17**.**

Let $f(z)$ be a Cauchy-Stieltjes function and $X=\mathrm{vec}^{-1}(f(\mathcal{M})\mathrm{vec}(F))$, with $\mathcal{M}:=I\otimes A-B^{T}\otimes I$, where $A,-B$ are Hermitian positive definite with spectra contained in $[a,b]$. Then the singular values $\sigma_{j}(X)$ of the matrix $X$ verifies:

[TABLE]

Proof 4.18**.**

Let us define $\widehat{X}_{\ell}:=\int_{0}^{\infty}X^{\text{ADI}}_{\ell}(t)\mu(t)dt=\int_{0}^{\infty}U^{\text{ADI}}_{\ell}(t)\mu(t)dt\cdot Y^{\text{ADI}}_{\ell}(V^{\text{ADI}}_{\ell})^{T}$. Since $V^{\text{ADI}}_{\ell}$ does not depend on $t$ we can take it out from the integral, and therefore $\widehat{X}_{\ell}$ has rank bounded by $\ell k$. Then, applying the Eckart-Young theorem we have the inequality

[TABLE]

5 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. The theory of equidistributed sequences has been used to obtained a nested sequence of poles with the same asymptotic convergence rate. Then, the approach presented in [8] 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 $\mathrm{vec}^{-1}(f(I\otimes A-B^{T}\otimes I)\mathrm{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. Also in the Kronecker structured case, it has been shown that relying on equidistributed sequences is an effective practical choice.

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.

Acknowledgment

The author wish to thank Paul Van Dooren and André Ran for fruitful discussions about Lemma 3.6.

Appendix A Proof of Theorem 4.4

According to [3, Theorem 2.1], the residue $R:=AX_{\ell}-X_{\ell}B-F$ can be written333In 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. as $\rho=\rho_{12}+\rho_{21}$, with

[TABLE]

where $r_{A}^{G}(z):=\det(zI-A_{\ell})/p(z)$, and $r_{B}^{G}(z)=\det(zI-B_{\ell})/q(z)$, with $p(z),q(z)$ defined as in (18). In addition, it is shown that $\rho_{12}=UU^{*}\rho(I-VV^{*})$ and $\rho_{21}=(I-UU^{*})\rho VV^{*}$.

Moreover, the proof of [3, Theorem 2.1] shows that, for any choice of $(\ell,\ell)$-rational function $r_{B}(x)$ with poles $z_{1},\ldots,z_{\ell}$, we can further decompose $\rho_{12}$ as $\rho_{12}=U(J_{1}-J_{2})$, where

[TABLE]

with $\mathcal{S}_{A,B}(X):=AX-XB$ and $\Gamma_{A}$ a path encircling once the interval $I_{A}$ but not $I_{B}$.

With a direct integration we get

[TABLE]

which yields $\lVert J_{1}\rVert_{2}\leqslant\lVert F\rVert_{2}\cdot\lVert r_{B}(A_{\ell})^{-1}\rVert_{2}\lVert r_{B}(B_{\ell})\rVert_{2}$. Let $\widetilde{B}:=VB_{\ell}V^{*}-c(I-VV^{*})$. Then,444The matrix $\widetilde{B}$ is not used in the original proof of [3], which contains a minor typo. There, the operator $\mathcal{S}_{A_{\ell},\widetilde{B}}$ is replaced by $\mathcal{S}_{A_{\ell},B_{\ell}}$ which does not have compatible dimensions.

[TABLE]

where we used that $V^{*}\widetilde{B}=B_{\ell}V^{*}$ and that the integral on the path $\Gamma_{A}$ of $(zI-B_{\ell})^{-1}/r_{B}(z)$ vanishes. Notice that $\lVert\mathcal{S}_{A,B}(X)\rVert_{2}\leqslant(\lVert A\rVert_{2}+\lVert B\rVert_{2})\lVert X\rVert_{2}$ and $\lVert\mathcal{S}_{A,B}^{-1}(X)\rVert_{2}\leqslant\lVert X\rVert_{2}/\min_{i,j}|\lambda_{i}(A)-\lambda_{j}(B)|$ [21, Theorem 2.1]. We get $\lVert J_{2}\rVert_{2}\leqslant\kappa\lVert r_{B}(A_{\ell})^{-1}\rVert_{2}\lVert r_{B}(B_{\ell})\rVert_{2}\lVert F\rVert_{2}$ and consequently

[TABLE]

Taking the minimum over all $(\ell,\ell)$-rational functions with poles $\Xi$ provides $\lVert\rho_{12}\rVert_{2}\leqslant(1+\kappa)\theta_{\ell}(I_{B},I_{A},\Xi)\lVert F\rVert_{2}.$ Analogously one obtains the similar estimate for $\rho_{21}$ swapping the role of $A$ and $B$. Since $\rho_{12}$ and $\rho_{21}$ have orthonormal rows and columns, we have $\lVert\rho_{12}+\rho_{21}\rVert_{2}=\max\{\lVert\rho_{12}\rVert_{2},\lVert\rho_{21}\rVert_{2}\}$, which concludes the proof.

Appendix B Bounding an inverse Laplace transform

The proof of Theorem 3.4 requires to bound the infinity norm of an inverse Laplace transform of a particular rational function, given in Lemma 3.6. The purpose of this appendix is to provide the details of its proof, that uses elementary arguments even though it is quite long.

Let us consider the following functions, usually called sine integral functions, that will be useful in the following proofs:

[TABLE]

It is known that $\mathrm{si}(x)+\mathrm{Si}(x)=\frac{\pi}{2}$, and that $0\leqslant\mathrm{Si}(x)\leqslant 1.852$ (see [1, Section 6.16]), and therefore $|\mathrm{si}(x)|\leqslant\frac{\pi}{2}$. We will need the following result, which involved integral of the sinc function by some particular measure.

Lemma B.1**.**

Let $g(t)$ be a decreasing and positive $\mathcal{C}^{1}$ function over an interval $[0,\gamma]$. Then, the following inequality holds:

[TABLE]

Proof B.2**.**

Integrating by parts yields I=\mathrm{Si}(s)g(s)\Big{|}_{0}^{\gamma}-\int_{0}^{\gamma}\mathrm{Si}(s)g^{\prime}(s)\ ds. The first term is equal to $\mathrm{Si}(\gamma)g(\gamma)$, which can be bounded by $1.852\cdot g(\gamma)$. The second part can be bounded in modulus with

[TABLE]

where we have used that $g^{\prime}(s)$ is negative, so $|g^{\prime}(s)|=-g^{\prime}(s)$. Combining the two inequalities we have

[TABLE]

Given a set of positive real points $\gamma_{j}$ enclosed in a interval $[a,b]$ with $a>0$, we define the rational function

[TABLE]

Note that $\widehat{r}(s)$ has poles enclosed in the negative half-plane which ensures that $\lim_{t\to\infty}\mathcal{L}^{-1}[\widehat{r}(s)]=1$. In particular $\mathcal{L}^{-1}[\widehat{r}(s)]$ is bounded on $\mathbb{R}_{+}$. We are now ready to prove Lemma 3.6.

Proof B.3** (Proof of Lemma 3.6).**

We write the inverse Laplace transform as follows:

[TABLE]

The integration path needs to be chosen to keep all the poles on its left, including zero. Therefore, we choose the path $\gamma_{\epsilon}$ that goes from $-\mathbf{i}T$ to $-\mathbf{i}\epsilon$, follows a semicircular path around [math] on the right, and then connects $\mathbf{i}\epsilon$ to $\mathbf{i}T$. It is sketched in the following figure:

-\mathbf{i}T$$\mathbf{i}T$$-\mathbf{i}\epsilon$$\mathbf{i}\epsilon

Splitting the integral in the three segments we obtain the formula:

[TABLE]

Concerning the first term, it is immediate to see that the integrand uniformly converges to the $1/s$ for $\epsilon\to 0$, and therefore the first terms goes to $\frac{1}{2}$ for $\epsilon$ small. We now focus on the second part.

We can rephrase the ratio of polynomials defining $r(s)$ as follows:

[TABLE]

Then, we note that the above ratio restricted to the points of the form $\mathbf{i}s$ yields a complex number of modulus one, that must have the form $\frac{p(\mathbf{i}s)}{p(-\mathbf{i}s)}=e^{\mathbf{i}\theta(s)}$, where

[TABLE]

In particular, $\lim_{s\to\infty}\theta(s)=-\ell\pi$ and for $s>0$ it holds

[TABLE]

This allows to rephrase the integrals of (24) in the more convenient form

[TABLE]

Since we are summing the integral between $[\epsilon,\infty]$ and $[-\infty,\epsilon]$ we can drop the odd part of the integrand, and rewrite their sum as follows:

[TABLE]

The above integral is well-defined even if we let $\epsilon\to 0$, we can can take the limit in (24) which yields exactly the value $\frac{1}{2}$ for the first term, and we have reduced the problem to estimate $f(t)=\frac{1}{2}+\frac{(-1)^{\ell}}{\pi}\int_{0}^{\infty}\frac{\sin(st+\theta(s))}{s}\ ds.$ To bound the integral, we split the integration domain in three parts:

[TABLE]

where we choose $\nu=a\tan(\frac{\pi}{4\ell})$ and $\xi=4\ell b$. For later use, we note that $\frac{a\pi}{4\ell}\leqslant\nu\leqslant\frac{a}{\ell}$. Concerning $I_{1}$, we can further split the integral as $I_{1}=\frac{1}{\pi}\int_{0}^{\nu}\frac{\sin(st)\cos(\theta(s))}{s}ds+\frac{1}{\pi}\int_{0}^{\nu}\frac{\cos(st)\sin(\theta(s))}{s}ds.$ Note that $|\theta(s)|\leqslant 2s\sum_{j=1}^{\ell}\gamma_{j}^{-1}$, which can be obtained making use of the inequality $|\mathrm{atan}(t)|\leqslant|t|$. We can bound the second integral term as follows:

[TABLE]

where we have used $\nu\leqslant\frac{a}{\ell}$ and $\sum_{j}^{\ell}\gamma_{j}^{-1}\leqslant\frac{\ell}{a}$. The first part can be bounded making use of Lemma B.1, by introducing the change of variable $y=st$, which yields

[TABLE]

Note that on $[0,t\nu]$ the function $\cos(\theta(y/t))$ is indeed decreasing, thanks to our choice of $\nu$, and therefore the above can be bounded in modulus by $\frac{1}{\pi}\left|\int_{0}^{t\nu}\frac{\sin(y)\cos(\theta(y/t))}{y}dy\right|\leqslant\frac{1.852}{\pi},$ where we have used that $\cos(\theta(0))=1$, and applied Lemma B.1.

Concerning $I_{2}$ we have

[TABLE]

Concerning $I_{3}$, we perform the same splitting for a since of a sum that we had for $I_{1}$, yielding

[TABLE]

By using (26) we have that $\forall s\in[\xi,\infty)$:

[TABLE]

Using the Lagrange expression for the residual of the Taylor expansion we get $\cos(\varphi(s))\underbrace{=}_{\psi\in[0,\varphi(s)]}1-\sin(\psi(s))\varphi(s).$ This enables bounding $I_{4}$ as follows:

[TABLE]

Analogously, for bounding $I_{5}$, we remark that by using (26) we have that $\forall s\in[\xi,\infty)$:

[TABLE]

Hence,

[TABLE]

Combining all these inequalities, we have

[TABLE]