REVIEW 4 minor 1 cited by
Rational Krylov for Stieltjes matrix functions: convergence and pole selection
T0 review · 0 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read Rational Krylov poles chosen from the Zolotarev problem make the action of Laplace-Stieltjes and Cauchy-Stieltjes matrix functions converge exponentially, with rates expressed purely in terms of the spectral interval and with analogous…
desk verdict Solid rational Krylov pole-selection paper with explicit bounds for Laplace- and Cauchy-Stieltjes functions; the math is credible, and the main caveats are experimental documentation and an unproven sharp constant that does not affect the rates. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The engine of the argument is the reduction of the projection error to a rational approximation problem. For the resolvent kernel $(tI+A)^{-1}$, the skeleton approximation yields the explicit error formula $1-(t+\lambda)f_{\mathrm{skel}}(t,\lambda)=r(\lambda)/r(-t)$, and for the exponential kernel the Laplace transform converts the same skeleton error into a bound involving the inverse Laplace transform of $(1/s)p(s)/p(-s)$. The optimal rational functions solving the third Zolotarev problem, that is, the poles of $p(z)/p(-z)$ that minimize $\max_{[a,b]}|r(z)|/\min_{[a,b]}|r(-z)|$, provide the pole sets; Möbius transformations map the relevant intervals $[-\infty,0]\cup[a,b]$ onto the symmetric Zolotarev configuration. In the Kronecker case the machinery is the tensorized subspace $U\otimes V$, where the projected matrix retains the same Kronecker-sum form, and the Cauchy-Stieltjes analysis passes through uniformly accurate solution of parameter-dependent Sylvester equations $(tI+A)X_t-X_tB=F$ using rational Galerkin residual bounds.
What would settle it
Pick the 1D Laplacian $A$ of size $10^5$ with known spectrum in $[a,b]\approx[1/n,4]$, set $f(z)=e^{-z}$, and run rational Krylov with the paper's Zolotarev poles for $\ell=5,10,15,20$. The log-error should fall with slope at least $-\pi^2/(2\log(4b/a))$; a markedly flatter slope, or failure of the error to decay once past the constant regime, would falsify the bound.
Extended reading notes
Core claim
For Hermitian positive definite $A$ with spectrum in $[a,b]$, the rational Krylov iterate with Zolotarev poles satisfies $\|f(A)v-x_\ell\|_2\le C\rho_{[a,b]}^{\ell/2}$ when $f$ is Laplace-Stieltjes and $\|f(A)v-x_\ell\|_2\le C\rho_{[a,4b]}^{\ell}$ when $f$ is Cauchy-Stieltjes. For $M=I\otimes A-B^T\otimes I$ with $A,-B$ Hermitian positive definite and $v=\mathrm{vec}(F)$ for low-rank $F$, tensorized rational Krylov spaces give $\|X-X_\ell\|_2\le C\rho_{[a,b]}^{\ell/2}$ for Laplace-Stieltjes functions and $\|X-X_\ell\|_2\le C\rho_{[a,2b]}^{\ell}$ for Cauchy-Stieltjes functions. The constants depend only on $f(0+)$ or $f(2a)$, the spectral ratio $\kappa=b/a$, and a mild logarithmic factor $\gamma_{\ell,\kappa}$, but not on the density defining $f$. These are explicit, parameter-free convergence rates, and the paper also derives matching singular-value decay bounds for the low-rank solution $X$.
Load-bearing premise
The convergence rates assume $A$ and $-B$ are symmetric positive definite with spectra exactly contained in a known interval $[a,b]$, and the poles are built from those endpoints; if the interval endpoints are inaccurate, the predicted rates no longer apply.
Editorial extensions
If this is right
- Ill-conditioned arguments are where the gain appears: the rate depends on the spectral ratio only through $\log(4b/a)$, so a widening spectrum costs little extra.
- Nested sequences of poles with the same asymptotic rate can be generated by equidistributed sequences, so the method can be run adaptively without fixing $\ell$ in advance.
- For Laplace-Stieltjes functions and rank-$k$ right-hand sides, the solution $X$ of $\mathrm{vec}(X)=f(I\otimes A-B^T\otimes I)\mathrm{vec}(F)$ has singular values bounded by $O(\rho_{[a,b]}^{\ell/2})$ after the $k\ell$-th, justifying low-rank compression.
- For Cauchy-Stieltjes functions of Kronecker-sum matrices, the analogous singular-value decay $O(\rho_{[a,2b]}^{\ell})$ holds without the condition-number factor $(1+\kappa)$.
- The projected function evaluation can be done by diagonalizing $A_U$ and $B_V$, costing $O(k^3)$ plus $k^2$ scalar function evaluations rather than $O(k^6)$.
Reading between the lines
- A natural stress test is the non-normal case the authors leave open: if the field of values lies in the right half-plane, the same Zolotarev intervals should still yield some convergence, but the sharp rates and the diagonalization trick for the projected problem would need separate treatment.
- The bounds can be read as a stopping-criterion blueprint: since the error after $\ell$ steps is dominated by an explicit power of $\rho$, one can choose $\ell$ in advance from the spectrum and target tolerance, and the nested pole sequences let the same choice be revised online.
- The singular-value decay theorems imply that solutions of parameter-dependent Sylvester equations inherit low-rank approximability uniformly in the shift parameter, which is exactly what makes low-rank time-stepping for fractional diffusion feasible on tensor grids.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The manuscript develops pole selection strategies for rational Krylov subspace methods that approximate x = f(M)v when f is a Laplace-Stieltjes or Cauchy-Stieltjes function and M is positive definite. For a Hermitian positive definite A with spectrum contained in [a,b], the authors prove explicit convergence bounds: Laplace-Stieltjes functions converge like O(rho^{ell/2}_{[a,b]}) using Zolotarev poles (Corollary 3.16) and Cauchy-Stieltjes functions converge like O(rho^ell_{[a,4b]}) using a Mobius-transformed Zolotarev set (Corollary 3.18). They also construct nested pole sequences via equidistributed sequences (Section 3.5). For Kronecker-structured matrices M = I otimes A - B^T otimes I with low-rank right-hand sides, they analyze tensorized rational Krylov subspaces and prove analogous rates: O(rho^{ell/2}_{[a,b]}) for Laplace-Stieltjes functions (Corollary 4.1) and O(rho^ell_{[a,2b]}) for Cauchy-Stieltjes functions (Corollary 4.8). Finally, they derive singular-value decay bounds for the solution matrix X, extending known low-rank approximability results for Sylvester equations (Theorems 4.10 and 4.11). Numerical experiments on finite-difference matrices, including fractional diffusion and exponential integrator examples, support the predicted rates.
Significance. If the results are correct, this is a valuable contribution to the rational Krylov literature. The paper provides explicit, non-fitted convergence rates and concrete pole sets for two important classes of matrix functions, and it extends the analysis to Kronecker-structured problems with low-rank right-hand sides, yielding both practical tensorized rational Krylov algorithms and new low-rank approximability bounds. The proofs are built on standard machinery -- Zolotarev theory, skeleton approximation, inverse Laplace transform estimates, and ADI error analysis -- and the constants are explicit rather than fitted to numerical data. The authors are appropriately candid about the limits of the analysis: Remark 3.7 acknowledges that the constant gamma_{ell,kappa} in Lemma 3.5 is likely not sharp, and Section 5 states that the non-normal case is left open. These caveats do not undermine the main rate claims, because gamma_{ell,kappa} only affects constants and not the exponential factor rho^{ell/2}.
minor comments (4)
- [Section 3.2, Lemma 3.5 and Remark 3.7] The constant gamma_{ell,kappa} is presented as an upper bound, but Remark 3.7 reports numerical evidence suggesting that the true constant is 1 for all admissible poles. Since Corollaries 3.16, 4.1, and Theorem 4.10 all inherit this constant, the authors should state more explicitly that (i) the exponential rate in Table 1.1 is unaffected by gamma, and (ii) the logarithmic factor in gamma is subexponential and therefore does not change the O(rho^{ell/2}) classification when ell is treated as the growing parameter.
- [Section 4.3, proof of Corollary 4.8] The phrase 'By setting I_A = I, I_B = -I in the statement of Theorem 4.6' is confusing: I_A and I_B are intervals, not identity matrices. Please rewrite this sentence using the normalized intervals, e.g. I_A = [\tilde{a},1] and I_B = [-1,-\tilde{a}], so that the Mobius transformation argument is unambiguous.
- [Section 3.6 and Section 5] The numerical experiments on tridiag(-1,2,-1) use the exact spectral interval, so they do not probe how the method behaves when the available bounds on [a,b] are over- or underestimated. Since all theorems are stated in terms of the interval used to construct the poles, adding a sentence in Section 5 on practical estimation of [a,b] (e.g., by Lanczos or Ritz values) would help the reader apply the bounds in realistic settings.
- [Throughout] There are a few minor typos: the acknowledgment reads 'The author wish' instead of 'The authors wish', and the introduction contains a duplicated 'of of' in the phrase 'the solution of of Sylvester equations' (Section 1.1, item 5). These should be corrected in a final pass.
Circularity Check
No significant circularity: the convergence bounds are derived from external rational-approximation results and the Stieltjes integral representation, with no fitted constants or load-bearing self-citations.
full rationale
The paper's central claims are upper bounds on rational Krylov errors for Laplace-Stieltjes and Cauchy-Stieltjes matrix functions. The derivation chain is: (i) the Stieltjes integral representation (1.1) expresses f(A)v as an integral of matrix exponentials or resolvents; (ii) Theorem 3.4 bounds resolvent projection errors via the skeleton approximation from Oseledets [26] and Druskin et al. [14]; (iii) Theorem 3.6 bounds exponential projection errors via the Laplace transform and the inverse-Laplace estimate of Lemma 3.5, whose proof is included in Appendix B; (iv) the resulting integrals are bounded using the Zolotarev rational approximation estimate from Beckermann/Beckermann-Townsend [3,5] and Mobius invariance. The pole sets are chosen to minimize or nearly minimize explicit rational approximation problems, and the error bounds are parameter-free upper bounds: no constant is fitted to the reported numerical errors, and the rates are not used to define the poles. The numerical experiments compare measured errors against the theoretical bounds, but the bounds are not calibrated to the data. The self-citations [23-25] appear only in motivational paragraphs about low-rank solvers and are not used in the convergence theorems. Remark 3.7 openly notes that the constant gamma_{ell,kappa} in Lemma 3.5 may not be sharp and that experiments suggest gamma = 1; this is a potential looseness in the bound, not a circularity, because the bound remains a proven upper bound independent of the numerics. The skeptical concern about the saddle-point constant is therefore a sharpness question rather than a reduction of the prediction to its inputs. Overall, the derivation is self-contained given the cited external mathematical results, and no prediction is equivalent to an input by construction.
Assumptions & free parameters
free parameters (2)
- Equidistributed sequence seed zeta =
1/sqrt(2)
- Shift parameter eta for Laplace-Stieltjes functions with f(0+)=infinity =
user-chosen in (0,a)
assumptions (6)
- standard math Zolotarev's solution of the third Zolotarev problem for [-b,-a] union [a,b] with the estimate rho_{[a,b]} = exp(-pi^2/log(4b/a))
- standard math Bernstein's theorem identifying Laplace-Stieltjes functions with completely monotonic functions
- standard math The quasi-optimality bound (3.1) for Galerkin projection onto rational Krylov subspaces (Guettel [18])
- standard math The skeleton approximation error formula 1-(t+lambda)f_skel(t,lambda)=r(lambda)/r(-t) (Oseledets [26])
- standard math The inverse Laplace representation (3.3) and the ADI error representation for Sylvester equations (Beckermann [3], Beckermann-Townsend [5])
- domain assumption A (and -B) are Hermitian positive definite with spectra contained in [a,b]
Cite this review
Pith. "Pith review of Rational Krylov for Stieltjes matrix functions: convergence and pole selection." pith.science (2026). https://pith.science/paper/S6HOYCVZ
@misc{pith2026190802032,
author = {Pith},
title = {Pith review of: Rational Krylov for Stieltjes matrix functions: convergence and pole selection},
year = {2026},
howpublished = {\url{https://pith.science/paper/S6HOYCVZ}},
note = {Machine review of arXiv:1908.02032}
}
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.
Figures
Forward citations
Cited by 1 Pith paper
-
A Shared Observation Shields Collective Fluctuations while Preserving Local Independence
Conditioning an exchangeable population on its own noisy collective record produces a rank-one, nonpositive Schur shield correction to collective fluctuations while keeping pair correlations at O(1/z), giving a calcul...
Reference graph
Works this paper leans on
-
[1]
M. Abramowitz and I. A. Stegun , Handbook of mathematical functions: with formulas, graphs , and mathematical tables , vol. 55, Courier Corporation, 1965
work page 1965
-
[2]
R. H. Bartels and G. W. Stew art , Algorithm 432: Solution of the Matrix Equation AX + XB = C, Comm. ACM, 15 (1972), pp. 820–826
work page 1972
-
[3]
B. Beckermann , An error analysis for rational Galerkin projection applied to the Sylvester equation , SIAM J. Numer. Anal., 49 (2011), pp. 2430–2450, https://doi.org/10.1137/110824590, https://doi.org/10. 1137/110824590
-
[4]
B. Beckermann and L. Reichel , Error estimates and evaluation of matrix functions via the F aber trans- form, SIAM J. Numer. Anal., 47 (2009), pp. 3849–3883, https://doi.org/10.1137/080741744, https:// doi.org/10.1137/080741744
-
[5]
B. Beckermann and A. Townsend , Bounds on the singular values of matrices with displacement structure, SIAM Rev., 61 (2019), pp. 319–344, https://doi.org/10.1137/19M1244433, https://doi.org/10.1137/ 19M1244433. Revised reprint of ”On the singular values of matrices with displacement structure” [ MR3717820]
-
[6]
P. Benner, R.-C. Li, and N. Truhar , On the ADI method for Sylvester equations , J. Comput. Appl. Math., 233 (2009), pp. 1035–1045, https://doi.org/10.1016/j.cam.2009.08.108, https://doi.org/10.1016/j.cam. 2009.08.108
-
[7]
M. Benzi and C. Klymko , Total communicability as a centrality measure , Journal of Complex Networks, 1 (2013), pp. 124–149
work page 2013
-
[8]
M. Benzi and V. Simoncini , Approximation of functions of large matrices with Kronecke r structure , Nu- mer. Math., 135 (2017), pp. 1–26, https://doi.org/10.1007/s00211-016-0799-9, https://doi.org/10.1007/ s00211-016-0799-9
Show all 33 references
-
[9]
Berg , Stieltjes-Pick-Bernstein-Schoenberg and their connecti on to complete monotonicity , Positive Def- inite Functions: From Schoenberg to Space-Time Challenges , (2008), pp
C. Berg , Stieltjes-Pick-Bernstein-Schoenberg and their connecti on to complete monotonicity , Positive Def- inite Functions: From Schoenberg to Space-Time Challenges , (2008), pp. 15–45
2008
-
[10]
Bernstein , Sur les fonctions absolument monotones , Acta Math., 52 (1929), pp
S. Bernstein , Sur les fonctions absolument monotones , Acta Math., 52 (1929), pp. 1–66, https://doi.org/ 10.1007/BF02547400, https://doi.org/10.1007/BF02547400
1929 doi
-
[11]
Braess , Nonlinear approximation theory , vol
D. Braess , Nonlinear approximation theory , vol. 7, Springer Science & Business Media, 2012
2012
-
[12]
Breiten, V
T. Breiten, V. Simoncini, and M. Stoll , Low-rank solvers for fractional differential equations , Electron. Trans. Numer. Anal., 45 (2016), pp. 107–132
2016
-
[13]
Druskin and L
V. Druskin and L. Knizhnerman , Extended Krylov subspaces: approximation of the matrix squ are root and related functions , SIAM Journal on Matrix Analysis and Applications, 19 (1998 ), pp. 755–771
1998
-
[14]
Druskin, L
V. Druskin, L. Knizhnerman, and M. Zaslavsky , Solution of large scale evolutionary problems using rational Krylov subspaces with optimized shifts , SIAM J. Sci. Comput., 31 (2009), pp. 3760–3780, https://doi.org/10.1137/080742403, https://doi.org/10.1137/080742403
2009 doi
-
[15]
Druskin, C
V. Druskin, C. Lieberman, and M. Zaslavsky , On adaptive choice of shifts in rational Krylov subspace reduction of evolutionary problems , SIAM Journal on Scientific Computing, 32 (2010), pp. 2485–2 496
2010
-
[16]
Druskin and V
V. Druskin and V. Simoncini , Adaptive rational Krylov subspaces for large-scale dynami cal systems , Sys- tems Control Lett., 60 (2011), pp. 546–560, https://doi.org/10.1016/j.sysconle.2011.04.013, https://doi. org/10.1016/j.sysconle.2011.04.013
2011 doi
-
[17]
Frommer, S
A. Frommer, S. G ¨uttel, and M. Schweitzer , Efficient and stable Arnoldi restarts for matrix functions based on quadrature , SIAM J. Matrix Anal. Appl., 35 (2014), pp. 661–683, https://doi.org/10.1137/ 13093491X, https://doi.org/10.1137/13093491X
2014 doi
-
[18]
G ¨uttel, Rational Krylov approximation of matrix functions: numeri cal methods and optimal pole se- lection, GAMM-Mitt., 36 (2013), pp
S. G ¨uttel, Rational Krylov approximation of matrix functions: numeri cal methods and optimal pole se- lection, GAMM-Mitt., 36 (2013), pp. 8–31, https://doi.org/10.1002/gamm.201310002, https://doi.org/ 10.1002/gamm.201310002
2013 doi
-
[19]
G ¨uttel and L
S. G ¨uttel and L. Knizhnerman , A black-box rational Arnoldi variant for Cauchy-Stieltjes matrix func- tions, BIT, 53 (2013), pp. 595–616, https://doi.org/10.1007/s10543-013-0420-x, https://doi.org/10. 1007/s10543-013-0420-x
2013 doi
-
[20]
Hochbruck and A
M. Hochbruck and A. Ostermann , Exponential integrators, Acta Numer., 19 (2010), pp. 209–286, https:// doi.org/10.1017/S0962492910000048, https://doi.org/10.1017/S0962492910000048
2010 doi
-
[21]
R. A. Horn and F. Kittaneh , Two applications of a bound on the Hadamard product with a Cau chy matrix , Electron. J. Linear Algebra, 3 (1998), pp. 4–12, https://doi.org/10.13001/1081-3810.1010, https://doi. org/10.13001/1081-3810.1010. Dedicated to Hans Schneider on the occas...
1998
-
[22]
G. A. Kalugin, D. J. Jeffrey, R. M. Corless, and P. B. Borwein , Stieltjes and other integral represen- tations for functions of Lambert W , Integral Transforms and Special Functions, 23 (2012), pp. 581–593
2012
-
[23]
Kressner, S
D. Kressner, S. Massei, and L. Robol , Low-rank updates and a divide-and-conquer method for linea r ma- trix equations, SIAM J. Sci. Comput., 41 (2019), pp. A848–A876, https://doi.org/10.1137/17M1161038, https://doi.org/10.1137/17M1161038
2019 doi
-
[24]
Massei, M
S. Massei, M. Mazza, and L. Robol , Fast solvers for two-dimensional fractional diffusion equa tions using rank structured matrices, SIAM J. Sci. Comput., 41 (2019), pp. A2627–A2656, https://doi.org/10.1137/ 30 STEF ANO MASSEI, LEONARDO ROBOL 18M1180803, https://doi.org/10.1137...
2019 doi
-
[25]
Massei, D
S. Massei, D. Palitta, and L. Robol , Solving rank-structured Sylvester and Lyapunov equations , SIAM J. Matrix Anal. Appl., 39 (2018), pp. 1564–1590, https://doi.org/10.1137/17M1157155, https://doi.org/ 10.1137/17M1157155
2018 doi
-
[26]
I. V. Oseledets , Lower bounds for separable approximations of the Hilbert ke rnel, Mat. Sb., 198 (2007), pp. 137–144, https://doi.org/10.1070/SM2007v198n03ABEH003842, https://doi.org/10.1070/ SM2007v198n03ABEH003842
2007 doi
-
[27]
M. Schweitzer , Restarting and error estimation in polynomial and extended Krylov subspace methods for the approximation of matrix functions , PhD thesis, Universit¨ atsbibliothek W uppertal, 2016
2016
-
[28]
Simoncini , Computational methods for linear matrix equations , SIAM Rev., 58 (2016), pp
V. Simoncini , Computational methods for linear matrix equations , SIAM Rev., 58 (2016), pp. 377–441, https://doi.org/10.1137/130912839, https://doi.org/10.1137/130912839
2016 doi
-
[29]
Susnjara, N
A. Susnjara, N. Perraudin, D. Kressner, and P. V andergheynst , Accelerated filtering on graphs using Lanczos method , arXiv preprint arXiv:1509.04537, (2015)
2015 arXiv
-
[30]
Townsend and S
A. Townsend and S. Olver , The automatic solution of partial differential equations us ing a global spectral method, J. Comput. Phys., 299 (2015), pp. 106–123, https://doi.org/10.1016/j.jcp.2015.06.031, https:// doi.org/10.1016/j.jcp.2015.06.031
2015 doi
-
[31]
D. V. Widder , The Laplace Transform , Princeton Mathematical Series, v. 6, Princeton Universit y Press, Princeton, N. J., 1941
1941
-
[32]
Q. Yang, I. Turner, F. Liu, and M. Ili ´c, Novel numerical methods for solving the time-space fractio nal diffusion equation in two dimensions , SIAM J. Sci. Comput., 33 (2011), pp. 1159–1180, https://doi. org/10.1137/100800634, https://doi.org/10.1137/100800634
2011 doi
-
[33]
Zolotarev , Application of elliptic functions to questions of function s deviating least and most from zero, Zap
E. Zolotarev , Application of elliptic functions to questions of function s deviating least and most from zero, Zap. Imp. Akad. Nauk. St. Petersburg, 30 (1877), pp. 1–59
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.