Nearly Optimal Approximation of Matrix Functions by the Lanczos Method
Authors: Noah Amsel
Tyler Chen
Anne Greenbaum
Cameron N Musco
Christopher Musco
Venue: NeurIPS 2024 spotlight (2024)
Keywords: linear algebra
lanczos
matrix function
krylov
iterative
rational
approximation
matvec
Primary area: other
TL;DR: We provide theoretical justification for why the Lanczos-FA method outperforms other methods for approximating the action of a matrix function on a vector, f(A)b.

Abstract: Approximating the action of a matrix function $f(\vec{A})$ on a vector $\vec{b}$ is an increasingly important primitive in machine learning, data science, and statistics, with applications such as sampling high dimensional Gaussians, Gaussian process regression and Bayesian inference, principle component analysis, and approximating Hessian spectral densities.
Over the past decade, a number of algorithms enjoying strong theoretical guarantees have been proposed for this task.
Many of the most successful belong to a family of algorithms called Krylov subspace methods.
Remarkably, a classic Krylov subspace method, called the Lanczos method for matrix functions (Lanczos-FA), frequently outperforms newer methods in practice. Our main result is a theoretical justification for this finding: we show that, for a natural class of rational functions, Lanczos-FA matches the error of the best possible Krylov subspace method up to a multiplicative approximation factor. 
The approximation factor depends on the degree of $f(x)$'s denominator and the condition number of $\vec{A}$, but not on the number of iterations $k$. Our result provides a strong justification for the excellent performance of Lanczos-FA, especially on functions that are well approximated by rationals, such as the matrix square root.

Nearly Optimal Approximation of Matrix Functions by the Lanczos Method
Noah Amsel"
Tyler Chen"
Anne Greenbaum Cameron Musco Christopher Musco"
Abstract
Approximating the action of matrix function f(A) vector b is an increas- ingly important primitive machine learning, data science, and statistics, with applications such sampling high dimensional Gaussians, Gaussian process re- gression and Bayesian inference, principle component analysis, and approximat- ing Hessian spectral densities. Over the past decade, number of algorithms joying strong theoretical guarantees have been proposed for this task. Many of the most successful belong family of algorithms called Krylov subspace methods. Remarkably, a classic Krylov subspace method, called the Lanczos method for matrix functions (Lanczos-FA), frequently outperforms newer methods in prac- tice. Our main result is theoretical justification for this finding: we show that, for natural class of rational functions, Lanezos-Fa matches the error of the best possible Krylov subspace method up to multiplicative approximation factor. The approximation factor depends on the degree of f(x)'s denominator and the con- dition number of A, but not on the number of iterations k. Our result provides strong justification for the excellent performance of Lanezos- especially on functions that are well approximated by rationals, such the matrix square root.
Introduction
Given symmetric matrix A E Raxd with cigendecomposition. A Σd Auu.u. the matrix function f(A) corresponding scalar function f - is defined as
f(A) Σ f(i)u.u
Matrix functions arise throughout machine learning, data science, and statistics. For instance, the matrix square root is used in sampling Gaussians, Bayesian modeling, and Gaussian processes [4, 5, 59], general fractional matrix powers are used in Markov chain modeling and Rényi entropy estima- tion [68, 42, 43, 19]. the matrix logarithm is used for determinantal point processes, kernel learning, and approximating deg-determinants for Gaussian process regression and Bayesian inference [17, 39, 31], the matrix sign is used in principal components regression and spectral density estimation [29, 18, 44, 58, 32, 69, 13, 7, 14], and the matrix exponential is used network science [3, 67, 48). In many of these applications, on a vector E 2 i.c., f(A)b. This task can be performed much more efficiently than computing an eigendecomposition of and forming f(A) using (1).
do not need to compute the matrix ((A itself, but rather its action
Perhaps the first general purpose method for approximating (AA) is the Lanezos method for matrix function approximation (Lanczos-FA) [22, 30), which is the focus of the present paper. Over the past decade, a number of special purpose algorithms, designed for single function or class functions,
New York University (noah. amsel@nyu. edu, tyler chen@nyu. edu, cmusco@nyu.edu.
University of Washington (greenbauêuw.c edu) t University Massachusetts Amherst (cmusco0cs umass.edu)
38th Conference on Neural Information Processing Systems (NcurIPS 2024).

---

have been developed [23, 16, 26, 44, 11]. These newer algorithms often satisfy strong theoretical guarantees better than the best-known bounds for Lanczos-FA. The present paper is motivated by the following remarkable observation:
Despite being arguably the simplest and most general algorithm for computing f(A)b, Lanczos-FA frequently outperforms special purpose algorithms, some- times by orders of magnitude, on common test problems.
For instance, in Figures 4, 5, and we compare Lanczos-FA with specialized algorithms [44, 59, 11] that satisfy the best-known theoretical guarantees for computing f(A)b for the particular functions they were designed for. In these experiments Lanczos-FA drastically outperforms these methods, despite the fact that it was not designed for any particular function. Because of its outstanding performance, Lanczos-FA is perhaps the most commonly used algorithm for computing f(A)b in practice. The main goal of this paper is to improve our theoretical understanding of why Lanczos-FA performs so well, in order to help close the theory-practice gap.
1.1 Krylov subspace methods
Lanczos-FA falls into class of algorithms called Krylov Subspace Methods (KSMs). KSMs among the most powerful and widely used algorithms for a broad range of computational tasks in- cluding solving linear systems, computing eigenvalues/vectors, and low-rank approximation [62, 38, 47, 65]. Like other KSMs for computing I (A)b, Lanczos-FA iteratively constructs an approxima- tion from the Krylov subspace
Kg (A, b) == (b, A b).
(2)
The Lanczos algorithm [46] produces orthonormal basis Q [qi..... q& ] for the Krylov sub- space (A, b) and symmetric tridiagonal matrix T satisfying T = Q AQ [62]. The Lanczos-FA algorithm uses this Q and to approximate f(A)b.
Definition 1. The Lanczos-FA iterate for a problem instance (f. A. k) is defined
lan, (f: == Q
where Q and T are above.
In our analysis we assume exact arithmetic. We do not discuss the implementation of Lanczos or Lanczos-FA since there are many resources on this topic; see for instance [52, 8, 9]. Fortunately, the Lanczos algorithm is implemented properly, its finite-precision behavior closely follows its exact-arithmetic behavior nearby problem [33]. See Section 5 for further discussion.
1.2 Optimality guarantees for Krylov subspace methods
For wide variety of problem instances, Lanczos-FA observed to converge almost as quickly as the best approximation to f(A)b that could be returned by any KSM run for the same number of iterations. In particular, all KSMs output approximations that lie the span of Kg (A, b), that is, approximations of the form p for polynomial of degree less than k. Given a problem instance (f, k), the best possible approximation returned by a Krylov method is*
opt, (/: A, b) == argmin IIí(A)b - x||2-
(A.b)
By definition, the error of this Krylov optimal iterate can be characterized as follows:
(A) opt, (f: A. b)ll2
min IIf(A) - p(A)bbl2.
(3)
deg(p)<k
For general matrix functions, no efficient algorithm for computing opt, (f: is known, but as shown in Figure 1, the solution returned by Lanczos-FA often nearly matches the error this best approximation. It is thus natural to ask if, at least for some class of problems, Lanczos-FA satisfies following strong notion of approximate optimality:
We choose to measure error in the Euclidean norm, and thus define optimality with respect to this norm. While other norms have been considered in prior work on matrix matrix-function function approximation [1]]. the Euclidean simple starting point.

---

A b
exp(-A/10)
log(A)b
10'
10
10-16
Rolatiee Error
10 33
10:21
36
10-4h
10:3
10
10-91
Fact
10"6
10-122
10-
Lanczos-FA
Instance Optimal
Optimality Ratio
1.5
1.1
1.0
40
20
Number of itterations (k)
Figure Lanczos-FA error (A) lan (f: A. b)ll2 each iteration for several functionssspectra. "Instance Optimal" is the right hand side of Definition 2 with C' 1, which is lower bound for all KSMs, including Lanczos-FA, Lanczos-FA performs nearly instance optimally on wide range of problems, far better than Fact 3 predicts. This is easily seen in the bottom plots, which show the ratio the error of the Lanczos-I'A iterate and Krylov optimal iterate, opt, (f:A.b).
Definition 2 (Near Instance Optimality). For problem instance (f, A. b, k), we say that a Krylov method is nearly instance optimal with parameters C and c if
IIf(A) alg, (f;, C
min deg(p)<ckk
(If(A)b (F(A)b-p(A)bllz
Above, alg, (f: A, b) denotes the output of an algorithm (e.g., Lanczos-FAl obtained from the Krylov subspace K2 (A, b), ie., the output after k iterutions.
c s allow some slack comparing to the Krylov optimal iterate. The right hand side depends on the entire problem instance, (f:A, b), which is why we call the guarantee "instance optimal".
Definition 2,
1.3 Existing near-optimality analyses of Lanczos-FA
The best bound for Lanczos-FA applying to broad class of functions is the following common bound for Lanczos-FA: see for instance [61, 54, 52, 10].
Fact 3. For all problem instances (f.A. k), Lanczos-FA FA satisfies
- -fg (f:A. s 2|[b||2 min deg( (mm
max
1f(x) ((())-((x)))
Fact 3 is in some sense optimality guarantee; it compares the convergence of Lanczos-FA to the best possible uniform polynomial approximation fx).. However, it does not take into account properties of A such as isolated or clustered eigenvalues, and seen in Figure 1, it typically only gives loose upper bound on the performance of Lanczos-FA.
To date, near-instance-optinalit guarantees for Lanczos-FA akin to those of Definition are known only few special cases. The most well-known when f(x) = 1/x and is positive definite, in which case Lanczos-FA is mathematically equivalent to the celebrated Conjugate Gradient al- gorithm, and therefore exactly optimal in the A-norm [47]. The instance-optimality guarantee for Lanczos-FA FA in this setting is used to prove well-known super-exponential convergence in certain settings [6, 47, 8]. contrast, bound like Fact only provides exponential convergence and is widely understood by the numerical linear algebra community to not accurately describe the actual behavior of the algorithm in most cases [35, 47, 8). The only other near-optimality guarantees for Lanczos-FA of which we are aware concern A b for nonsymmetric A [12] and the matrix exponen- tial exp(-/A)b [20]. both cases, Lanczos-FA satisfies guarantees that are reminiscent of (although weaker than) Definition 2. Further discussion is given in Appendix B.1.

---

We also remark that there are number of works which aim to relate the convergence of of Lanczos- FA for functions with certain integral representations to the convergence of conjugate gradient on linear systems [41, 24, 27, 28, 25, 10]. While these analyses provides spectrum dependent conver- gence guarantees, they are weaker than Definition 2 because KSM approximation to a given function relates to the convergence of conjugate gradient. These bounds were mostly developed for use as a posteriori stopping criteria rather than as a theoretical explication for the behavior of Lanczos-FA.
not clear how the best possible
1.4 Our contributions
In Section we prove near instance optimality for broad class of rational functions (Theorem 4). To the best of our knowledge, this is the first true instance-opimmalit guarantee for Lanczos-FA -FA to be proven for any function besides f(x) 1/x. In Section 2.2, we discuss how results of this kind imply related guarantees for functions that are uniformly well-approximated by rationals, which includes many functions of interest in machine learning. In Section 3, we additionally show that Lanczos- FA satisfies a weaker version of near optimality for two crucial non-rational matrix functions the square root and inverse square root (Theorems and 7). Appendix of optimality to Definition 2 (near instance optimality) and to that of Fact 3. In Section 4, we present experimental evidence showing that, for many natural problem instances, our bounds are significantly sharper than the standard bound of Fact 3. These experiments also demonstrate that despite its generality, Lanczos-Fa FA often converges significantly faster than other methods practice. We conclude with discussion of next steps and open problems in Section 5.
compares the this version
2 Near optimality for rational functions
In this section, we study the Lanczos-FA approximation to r(A)b, where r(x) is rational function with raal-valued poles that lie in R \I, where I := [Amin mmx)] Specifically,
r(x) := n(xx//(xx),
where n(x) is any degree p polynomial and m(x) = (x - 2x)(x z2) ... (x zq),
Ij GERI
Since 2j +(A 2 2,1) is either positive definite or negative definite. For convenience, define
+ 1) (A 1
2j < d
Ay
Zj Amax
Our main result on rational functions the following pear-instance-optimaatyy bound, which holds under mild assumption the number of iterations k:
Theorem 4. Let r(x) n(x)/m(x)) be a degree (p, q) -rational function as in (4) and define A as in (5). Then, if max q 1}, the Lanczos-FA iterate satisfies the bound
min [|[(A)b lan (r: A. b) $ x(A1) x(A2). ... ...(Ag))) deg(p)<k-q+)
[|r(A)b- p(A)bll2-
We prove this theorem Appendix A. Above, x(A;) is the condition number of A,, the ratio of the largest to smallest magnitude eigenvalues of A,- Theorem 4 shows that Lanczos-FA used to approximate r r(A)b satisfies near instance optimality, as in Definition 2, with
C = q <(Ai) x(A2) x(Ag).
c= c= (1- 1)/k.
In particular, when A is positive definite and each of the 2 negative (as is the case for rational function approximations of many functions including the square root [37]), then C S qk(A)*.
As discussed in Section 5, we expect that Theorem 4 can be tightened by significantly reducing the prefactor. Nevertheless, as illustrated in Figure 2, even in its current form, the bound improves on the standard bound of Fact 3 in many natural cases: i.e., it more tightly characterizes the observed convergence of Lanezos-FA Finally, as discussed further tional functions being an interesting function class in their own right, Theorem 4 has implications for understanding the convergence of Lanczos-FA for functions like the square root which are much casier to approximate with rational functions than polynomials.
Section 2.2, we note that, beyond ra-

---

A
(A)b exp(-A/10jp (deg=S)
r(A)b log(A)b (deg=10)
10a
10-10
10-a
Relative Eroor
10
10-4
10-14
10-3
Fact
10
10
Theorem
10-6
10-11
10-57
Lanezos -A
Instance Optimal
Number of iterations (R)
Figure Despite its large prefactor, the bound of Theorem 4 qualitatively captures the convergence behavior of Lanczos-FA for rational functions. It can be tighter than the standard bound of Fact 3, even for a moderate number of iterations k. We use rational approximations to exp(-x/10) and log(x) for comparison with Figure I; see Section 4 for more details.
2.1 Proof Sketch
Our proof of Theorem 4 leverages the near optimality of Lanczos-Fi FA in computing A b to obtain bound for general rational functions. For illustration, consider the simplest possible rational function for which we are unaware of any previous near-optimality bounds: A - when A is positive definite. The Lanczos-FA upproximation for this function is QT = QT QTQT=QQT Using the triangle inequality and submultiplicatiiity, we can bound the error of the approximation
||A b lan, (x
b)ll2
||A b
QT Q A bll2 + QQT - Q
|| b QT QQ AA bbll QQT-'QT1l: ||A-'b- QT -
AA -AT''AAAA-1 AA
Using the normal equations and the fact that Q is a basis of the Krylov subspace, it is possible to show that Q QQQTAQQ* AQ) -A the A-norm projector onto Krylov subspace. Hence
||A - - QT QT - bblA min ||A - x|Ia x=XR(A,b)
min ||A - p(A)blAA deg(p)<k
Therefore, using ||x||2 (1/VAmin) |x]|A and ||x||A Vdmax ||]2, the first term the right hand side of (6) can be bounded as
|A b QT bbll s Vk(A)
min ||A p(A)bll2. dk
Similarly, the second factor of the second term on the right hand side of (6) can be bounded as
| Q & Vx(A)
min deg(p)<
||A-'D (A-'D--(A)bll>> (AA)bll2.
Then, since xp(x) is polynomial one degree larger than p(x),
min deg(
||A (0()bll S min IIAA -2b Ap(A)bll2
deg(p)<k-1
S Amax min deg(p)<k-1
-- (AA)bl22
Plugging the above bounds into (6) and using the fact that the eigenvalues T are contained in [dmins Amax] so that ||QT-IQT -TT T I| $ 1/dmin-
A 2 lan (x 2, b) 2x(A)3/2 deg(p)<k-1
min 2 p(A)bll2.
5 'The A noro is defined for positive definite A by ||x||A J|A 1/2x/2
5

---

Bounding x(A) by x(A)² gives the bound in Theorem 4.
Our proof of Theorem 4 generalizes the above approach. We write the Lanczos-FA error for ap- proximating r(A)b in terms of the error of the optimal approximations to a set of simpler rational functions, and then reduce polynomial approximation of r(x) to polynomial approximation of each of these particular functions.
2.2 Implications for non-rational functions
Theorem 4 can be used to derive guarantees for other (non-rational) functions. In particular, consider any function f(xx) that is uniformly well approximated by low-degree rational function r(x) on I [mmin- Amax] the interval containing all of the eigenvalues of A; i... suppose ||r ffllz == maxxer [xx) is small. A natural way to approximate f(A)b, used in [4, 5, 591. is to construct r(x) and output r(A)b, using some iterative linear solver to quickly apply the denominators of the partial fraction decomposition of r r(x). Alternatively, if we have an instance-opoe instance-optimality guarantee for Lanczos-FA on rational functions like Theorem 4, we could use Lanczos-FA The following analysis shows that simply using Lanczos-FA on fx) itself cannot be much worse.
compute r( (A)bb.
Lemma 5. Assume the we have the following instance-optimality guarantee for rational functions:
[|r() lan (r: A, b)lle S C, min deg(p)<cck
[|r(A)b - p(A)bil2-
(7)
Here, Cr and cr depend on the choice of approximant r(x). Then the error of Lanczos-FA on f(x) is bounded as follows:
II((A)b lan (f: A. b)ll2
S min (C,, +2) Ilʃ - -rll + C, min degk
I1/(A)D- - P(A)bll2).
(8)
We prove this lemma, which follows directly from the triangle inequality, in Appendix A.5. Com- pare the form of this bound to that of Definitions and 13. It is close well-approximated by a guarantee, except for the first term, which requires f(x) to be uniformly rational function r(x) on [mmins A This is still much stronger than Fact 3, which requires f(x) to be uniformly well-approximated by a polynomial to guarantee that Lanczos-FA FA provides a good approximation. There are many functions with lower degree rational approximations than polyno- mial approximations, even when we require the rational function r(x) to have poles only R \ (as in our Theorem 4). Such rational approximations are obtainable by the Remez algorithm [64, Chapter 24), and for many important functions they are also known explicitly. For example, a uni- form polynomial approximation the square root strictly positive interval [Amin: mmux requires degree 22 /mmax/Amin) [64, Chapter 8). On the other hand, uniform rational approximation can be obtained with degree only O(log(Amax/Amin
near-instancc-optimalll
[37, 64, Chapter 25J. Likewise, uniform polynomial approximation to exp(-x) (x) on the interval [0, B] requires degree Q () [1], but uniform rational ap- proximations can be constructed with no dependence on [63]. For such functions, we expect (8) to be stronger than Fact 3.
Notice also that, while in Lemma 5, we assume fxx is well-approximated by rational function, we are not required to actually construct the approximation. Indeed, since it holds for any r(x), instead of fixing rational approximation of certain degree, (8) automatically balances fll which decreases as the degree grows, with C,, which may increase as the degree grows (see Figure 4).
3 Near Spectrum Optimality for A 11/b
In the previous section, we proved that Lanczos-FA nearly instance optimal for rational functions in the sense of Definition 2. In this section, we prove that Lanczos-FÄ -FA satisfies weaker form of near optimality for two important non-rational functions: square root and inverse square root. We term this weaker form of guarantee "near spectrum optimality". In Appendix C, we formally define this notion and compare it to Definition 2. We first state our bound for the inverse square root.
Theorem 6. Let A be the spectrum of A. Then for k 2, the Lanczos-Ta FA iterate satisfies the bound
3
- -(xs))
A 1/2b lan, (x (x-12; 112
x( (A)|[bll2
min (e)(p)<kk22
max XEA
Vzk
VX

---

That is, Lanczos-FA applied to the inverse square root satisfies Definition 12 ("near spectrum opti- mality") with
3 VIk
1
C= <(A),
c=
We prove this theorem in Appendix D. The proof relies on comparing the error of the kth Lanczos iterate for 1/2 that of the Lanezos iterate for x- First, applying a bound from [10], we use the Cauchy integral formula upper bound the error of Lanczos-FA (Lemma 16). Second, as Equation (17) shows, Lanczos-FA is nearly instance optimal for the func- tion.x that is, outputs ((A) where p (nearly) the degree & polynomial that best approximates x= Third, the best degree & polynomial approximation to x must have lower error than the best degree 2/2 approximation to x 1/2 This is because any degree k/2 approximation to x 1/2 can be squared to yield good degree k polynomial approximation to x Combining these three steps upper bounds the error of the kth Lanczos-FA iterate for x- by the error of the best degree k/2 polynomial approximation of x 1/2
x=1/2 by its error on x-1
Nearly the same argument can be used to prove spectrum optimality of Lanczos-FA for the function A A for any л N with pendix C, we can convert Theorem 6 into a near- www.iinceooooimaaallieeee guarantee at the price of strong dependence of C on b. We next state our optimality result for the matrix square root.
(2"
x(A) (Ak/V Furthermore, using Lemma 15 of Ap-
Theorem 7. Let A be the spectrum of A. Then for 2, the Lanczos-PP -FA iterate satisfies the bound
3x(A)²
A 1/b lan (x A, S
|Bb]2
min deg(p)<k/24)
max (xEAU(0)
(VF-p(x))
R3/2
This bound resembles Definition 12 with
3x(A)²
C
c=
/3/2
However, it is slightly weaker in that the maximization is taken over A U (0) instead of only A.
The proof of Theorem 7 is nearly the same as that of Theorem 6, and it likewise appears in Ap-
pendix D. Ideally, if p is polynomial approximation x would like to claim that (p(x)xx)
yields a good polynomial approximation to x= However, since this function is not necessarily a
polynomial, must instead use use pxxx)-(())) which introduces the need to include {0} in the maximization on the right-hand side.
4 Experiments
We now present several numerical experiments to assess the quality of our instance-optimality bounds, Theorem and Lemma 5. Our results show that, despite the large prefactor C. our bounds already supersede the standard uniform approximation bound (Fact 3) in many cases. We also com- pare Lanezos-FA against several recently proposed algorithms for computing matrix functions with strong theoretical guarantees. We find that, in practice, Lanczos-FA FA performs better than all of them. We implement Lanczos-FA in high precision arithmetic using the flamp library, which is based on mpmath [51], in order to mitigate any potential impacts of finite precision arithmetic and observe the convergence behavior of the algorithm beyond the standard machine precision. 6
In Figure 1, we compare the performance of Lanczos to the instance-optimal KSM (which we can compute by direct methods) and against Fact for various matrix functions and spectra. We use three test matrices A E R 100x 100 which all have condition number 100. The first has a uniformly-spant spectrum, the second has a geometrically-spaeeee uniformly-space in a narrow interval except for ten very small eigenvalues. We compute the bound from Fact using the Remez algorithm and compute the instance-optimal approximation using least squares regression onto the Krylov basis Q. In Figure 1, in almost all cases we tried, Lanczos-FA FA
spectrum, and the third has eigenvalues that are all
Cooee for our experiments is available at httpge///gtttaoccco////aaeeee//p cree/neurips2024 Barpio24___ar_________ cottrrrttimaity.....
7

---

Dax Opitzoouryy Xaso (C)
10P
107
107
10P 101
10P
Figure 3: The maximum observed ratio between the error of Lanczos-Fii FA and the optimal error over choices of b when approximating A 4 for matrices with varying condition number Each point corresponds to a pair (k, q). Points with the same color have the same value K. On the left, the dotted line plots v@k for the maximum considered (10%). On the right, the dotted line plots vęk for the maximum q considered (2%). Overall, the optimality ratio appears to scale at least as Q2( VqK).
performs nearly as well the instance-opting optimal approximation. For instance, the error is never more than small multiple of the optimal error the experiments we did.
To better understand Theorem 4, we compare the bound to the true convergence curve of Lanczos- FA for various rational functions in Figure 2. We also plot Fact for reference. We use the same matrices and b vectors as in Figure I; results are similar if b is instead chosen as uniform linear combination of A's 's eigenvectors. We choose rational functions to match the functions used for Figure 1. We construct degree 5 rational approximation to exp(-x/Inn 10) following [63]. We construct degree 10 approximation to log(x) using the BRASIL algorithm [40] and verify that has real poles outside the interval of the eigenvalues. Despite Theorem 4's exponential dependence on the degree of the rational function being applicd, Figure 2 shows that matches the shape of the convergence curve well and is tighter than Fact 3 when the number of iterations is large. That said, in all cases plotted, Lanczos-Fi FA always returns an approximation much closer to optimal than predicted by Theorem 4, suggesting that the leading coefficient in our bound is pessimistic.
4.1 Dependence the rational function degree
Theorem 4 upper bounds the optimality ratio by C= (= <(A1) that this bound is far from tight. However, the following experiment provides evidence that it is not possible to entirely eliminate the dependence on the rational function's denominator degree q. In particular, for parameters (x, q), consider approximating A =b where A has spectrum d1 and .2..... evenly spaced between 0. 0.999995. x and k. In this case, k(A,) = = k(Ag) x(A) K. We generate grid of problems by picking different combinations of (k, Q) and tuning in limited way to maximize the maximum ratio between the error Lanczos-Fi FA and the optimal Krylov error
x(AA))... We conjecture
all iterations. In particular, its first entry, using grid search to maximize the optimality ratio. The results, plotted in Figure 3, suggest that the optimality ratio grows at least as 2 2(Vq x(A)). We have yet to find harder problem instances than this.
took b be all ones vector, except
varied
4.2 Non-rational functions
As noted in Section 2.2, an instance-optimality guarantee for rational functions also implies that Lanczos-FA performs well on functions that are well-approximated by rationals. As an example, we consider the function f(x) = x -0.4 for which a rational approximation in any degree can be found using the BRASIL algorithm [40]. Figure 4 shows how applying Lanczos-FA to these tional approximations compares to applying Lanczos-FA directly to the f(x) itself to approximate A -0.4b. When the number of iterations small, both methods perform nearly optimally, as the accuracy is limited more by the small size of the Krylov subspace than by the difference between f(x and the rational approximant (that is, second term in (8) dominates the first term). As number of iterations grows, the error due approximating f(x) in the Krylov subspace continues to decrease while the error of uniformly approximating fxx by some fixed rational function remains fixed (that is, the first term in (8) dominates the first term); however, increasing the degree of the

---

10-1
Relative Hroor
10-1
10-21
deg=5
deg=10
10-31
15g-15
deg-20
10
Lanezos-FA
10
15
20
25
30
35
40
Number of iterations (k)
Figure 4: Applying Lanczos-FA the function A - and rational approximations of various de- grees found using the BRASIL algorithm [40]. In this experiment, the spectrum A contains two clusters: 10 eigenvalues uniformly spaced near and 90 eigenvalues uniformly spaced near 100. As predicted by the bound in Section 2.2, convergence of Lanczos-Fa for this function appears to closely track that of high degree rational approximant.
Eigengup Uniform
Eigengap Skewed
No Eigengap Skewed
10'
10P
10
Rolative Eroor
10
10"
10-7
10
10
10
10 12
100
200
100
200
100
200
Lanezos FA
"slanczos" deg=16
"rational" deg=4
"slanczos" deg-32
"rational" deg-8
Number of matrix vector products (or equivalent in vector vector products)
Figure 5: A comparison Lanczos-Fa with two methods from [44] ("rational" and "slanczos") for computing the matrix sign function, which work by using a stochastic iterative method to ap- proximate rational approximations to the step function of various degrees. The "rational" method the main one studied in [44]. while "slanczos" is included because it is the best performing in their experiments. Each panel corresponds to one of the test problems from [44]. Iterations these meth- ods are counted in number of inner products with rows of A rather than number of matrix-vector products with A as a whole. To compare these with Lanczos-FA, we consider such inner products to be equivalent to matrix-vector product.
rational approximant decreases this source of error. This shows that understanding the convergence of Lanczos-FA for the entire family of rational approximations goes long way toward explaining its convergence for non-rational functions. In the limit, Lanczos-FA applied to fxx itself appears to automatically inherit the instance optimality of suitably high-degree rational approximation.
This result has additional implication. A number of papers use explicit rational approximations to compute non-rational matrix functions [4, 5. 36, 59]. These approximations are often applied by applying conjugate gradient (or related method) to each of the terms in the sum [36, 59]. In the case conjugate gradient used, the resulting algorithm mathematically identical to Lanczos- FA used compute the the rational approximation. However, Figure 4 suggests that simply using Lanczos-FA on the original function both simpler and converges faster (though memory usage and other factors may need be considered).
Another line of work uses specialized iterative methods that exploit problem structure to apply the rational approximations [29, 2, 53, 44J. In Figure 5, compare Lanczos-FA to two such methods from [44] for computing SAg,(A), for A of the form A = B AI. While they achieve better theoretical bounds than are known for Lanczos-FA. Lanczos-FA far outperforms them on the test problems used in [44].

---

sign(A)b
Ah
(5 a A2
100
Rdlatve Eitor
10"2
10
10"1
10
10-
Lanczup RA
10-1
10-1
Optinal
100
20000
20000
20
Figure 6: Convergence of Lanczos-Fa for rational functions with poles in A's eigenvalue range or that are imaginary. The optimality ratio be very large for some iterations. Similar behavior is seen for functions like sign(A) A) that have discontinuity in the interval of 's eigenvalues. However, the "overall" convergence of Lanczos-FF -FA still appears to closely track the instance-opting solution.
Additional Experiments. In Section 3 we introduce bounds for the matrix square root and inverse square root, and in Appendix E.1 provide numerical tests study the sharpness of the bounds, verifying that they can improve on Fact In Appendix E.2, we demonstrate the convergence be- havior of Lanczos-FA
rational functions with poles inside the range of eigenvalues (Figure 6). This illustrates why a bound like Theorem is not possible, but suggests a weaker bound, such as the bound in [12] for f(x) 1/x and indefinite A, may be possible. Appendix E.3 shows that, un- like Lanczos-FA, a related algorithm called Lanczos-OR [1 [11] (which is exactly optimal for rational functions, though not in the Euclidean norm) can perform poorly on high degree rational functions
when the error is measured in the Euclidean norm.
Outlook
This paper provides instancy-optimality guarantees for Lanczos-FA applied to range of rational functions. We conclude with open questions that we believe are worthy of further study.
Extension to other function classes. Empirically, Lanczos-FA seems to be nearly instance opti- mal for wide variety of functions beyond those considered in this paper, such rational functions with conjugate pairs of complex poles whose real parts lie outside [Amin- 6, the error of Lanezos-FA between being very large and nearly optimal. We discuss this more in Appendix E.2. would be valuable provide bounds explaining these behaviors.
As seen in Fig- functions with real poles in [Amin- Amax] is intriguing, oscillating
would be also natural to try to extend our bounds to Sticltjes/Markov functions, which can be viewed as a certain type of infinite degree rational function approximations with poles in (-oo, 0], and includes important functions like the inverse square root and shifted logartihm. Our bound Theorem 4 cannot be directly applied to this class due to the dependence the rational function denominator degree q.
Construction of hard instances refined upper bounds. Theorem has an exponential depen- dence on the degree of the rational function's 's denominator 4, which limits the practicality of our bounds. It unclear if and when this dependence can be improved. The experiment in Section 4.1 provides strong evidence that some dependence have depend on vq. instead of the current bound of <(A) guaranteed by Theorem 4. is an open question whether Theorem can be tightened, or whether matching hard instances exist.
is necessary, but the hardest examples
Finite precision arithmetic. Our analysis concerns the behavior of the Lanczos algorithm when run in exact arithmetic. In practice, the implementation of the Lanczos algorithm very important; for instance, practical implementations often output a Q which is far from orthogonal [50, 8). While this instability can be mitigated with more expensive implementations, theoretical work shows that, surprisingly, Lanczos and Lanczos-FA -AA can work well despite it [57, 56, 55, 33, 21, 20]. For example, finite precision arithmetic
[21, 52, 9] show that Fact 3 still holds up to close approximation for any bounded matrix function. It would be valuable to study whether stronger near-optimality guarantees like those proven Theorem 4 are also robust to finite precision. For f(x) 1/x, this problem has been studied in [33].
10

---

Acknowledgments: This research was supported in part by NSF Awards 2045590 (Chen, Ch. Musco), 2046235 (Ca. Musco). 2427363 (Chen, Ca. Musco, Ch. Musco), and 2234660 (Amsel).
References
[1] Amol Aggarwal and Josh Alman. "Optimal-Degree Polynomial Approximations for Expo- nentials and Gaussian Kernel Density Estimation". In: 37th Computational Complexity Con- ference (CCC 2022). 2022 (cit. on p. 6).
[2] Zeyuan Allen-Zhu and Yuanzhi Li. "Faster Principal Component Regression and Stable Ma- trix Chebyshev Approximation". In: Proceedings of the 34th International Conference on Machine Learning- Vol. 70. 2017, pp- 107-115 (cit. on p. 9).
[3] Mary Aprahamian, Desmond J. Higham, and Nicholas J. Higham. "Matching exponential- based and resolvent-
Networks 4.2 (2015), pp. 157-176. eprint: https: academicc oup. com/ comnet/ article-pdf 4/1////22///// 6716204/cnv01s pdf (cit. on p. 1).
centrality measures". In: Journal of Complex
[4] Erlend Aune, Jo Eidsvik, and Yvo Pokern. "Iterative numerical methods for sampling from high dimensional Gaussian distributions". In: Statistics and Computing 23 20 pp. 501-
521 (cit. on pp. 1, 6, 9).
[5] Erlend Aune, Daniel Simpson, and Jo Eidsvik. "Parameter estimation in high dimensional Gaussian distributions". In: Statistics and Computing 24 (2014), pp. 247-263 (cit. on pp. 1,
6, 9).
[6] Bernhard Beckermann and Aroo J. Kuijlaars. "Superlinear Convergence of Conjugate Gra- dients". In: SIAM Journal on Numerical Analysis 39.1 (2001). pp. 300-329 (cit. on pp. 3, 18).
[7] Vladimir Braverman, Aditya Krishnan, and Christopher Musco. "Sublinear time spectral den- sity estimation". In: Proceedings of the 54th Annual ACM SIGACT Symposium Theory of Computing (STOC 2022). 2022 (cit. on pp. 1, 1,21).
[8] Erin Carson, Jörg Liesen, and Zdenčk Strakos. "Towards understanding CG and GMRES through examples". In: Linear Algebra and its Applications 692 (2024), pp. 241-291 (cit. on pp. 2. 3, 10, 18).
[9] Tyler Chen. The Lanczos algorithm for matrix functions: handbook for scientists. 2024. arviv: 2410.11090 11090 [math. NA] (cit. on pp. 2, 10).
[10] Tyler Chen, Anne Greenbaum, Cameron Musco, and Christopher Musco. "Error Bounds for Lanczos-Based Matrix Function Approximation" In: SIAM Journal Matrix Analysis and Applications 43.2 (2022), pp. 787-811 (cit, on pp. 3, 4, 72)
[11] Tyler Chen, Anne Greenbaum, Cameron Musco, and Christopher Musco. Krylov Subspace Methods for Optimal Rational Matrix Function Approximation" In: SIAM Journal on Matrix Analysis and Applications 44.2 (2023), pp. 670-692 (cit. on pp. 10, 19, 24).
[12] Tyler Chen and Gérard Meurant. "Near-optimal convergence of the full orthogonalization method". In: Electronic Transactions Numerical Analysis 60 (2024), pp. 421-427. arXiv: 2403. 07259 [math.NA] (cit. on pp. 3, 10, 18, 24).
[13] Tyler Chen, Thomas Trogdon, and Shashanka Ubaru. "Analysis of stochastic Lanczos quadra- ture for spectrum approximation". In: Proceedings of the 37th International Conference Machine Learning (ICMI. 202 I.. 2021 (cit. on pp. 1,
[14] Tyler Chen, Thomas Trogdon, and Shashanka Ubaru. Randomized matrix-free quadrature for spectrum and spectral sum approximation. 2022. arXiv: 2204 01941 [math. NA] (cit. on pp. 21).
[15] Jane Cullum and Anne Greenbaum. "Relations between Galerkin and Norm-Minimizzin It- erative Methods for Solving Linear Systems". In: SIAM Journal on Matrix Analysis and Ap- plications 17.2 (1996), pp. 223 223-247 (cit. on p. 18).
[16] Nigel Cundy, Stefan Krieg, Guido Arnold, Andreas Frommer, Th Lippert, and Klaus Schilling. "Numerical methods for the QCD overlap operator IV: Hybrid Monte Carlo".. In: Computer physics communications 180.1 (2009), pp. 26-54 (cit. on p. 2).
11

---

[17] Kun Dong, David Eriksson, Hannes Nickisch, David Bindel, and Andrew G Wilson. "Scalable Log Determinants for Gaussian Process Kernel Learning.. In: Advances in Neural Informa- tion Processing Systems (NeurIPS 20/7). Vol. 30. 2017 (cit. on p. 1).
[18] Kun Dong, David Eriksson, Hannes Nickisch, David Bindel, and Andrew G Wilson. "Scalable Log Determinants for Gaussian Process Kernel Learning". In: Advances in Neural Informa- tion Processing Systems (NeurIPS 2017). Vol. 30. 2017 (cit. on p. I).
[19] Yuxin Dong. Tieliang Gong, Shujian Yu, and Chen Li. "Optimal Randomized Approxima- tions for Matrix- -Bsed Rényi's Entropy". In: IEEE Transactions on Information Theory 69.7 (2023), pp. 4218-423, (cit. on p. 1).
[20] Vladimir Druskin, Anne Greenbaum, and Leonid Knizhnerman. "Using Nonorthogonal Lanc- 70s Vectors in the Computation of Matrix Functions". In: SIAM Journal on Scientific Com- puting 19.1 (1998), pp. 38-54 (cit. on pp. 3, 10, 19).
[21] Vladimir Druskin and Leonid Knizhnerman. "Error Bounds in the Simple Lanczos Procedure for Computing Functions of Symmetric Matrices and Eigenvalues". In: Journal of Computa- tional Mathematics and Mathematical Physies 31.7 (1992). pp- 20-30 (cit. on p. 10).
[22] Vladimir Druskin and Leonid Knizhnerman. "Two polynomial methods of calculating func- tions of symmetric matrices". In: USSR Journal of Computational Mathematics and Mathe- matical Physics 29.6 (1989), pp. 112-121 (cit. on pp. 1, 15).
[23] Jasper van den Eshof, Andreas Frommer, Thomas Lippert, Klaus Schilling, and Henk A. van der Vorst. "Numerical methods for the QCD overlap operator. I. Sign-function and error bounds". In: Computer Physics Communications 146.2 (2002), pp. 203-224 (cit. on p. 2).
[24] Andreas Frommer, Stefan Güttel, and Marcel Schweitzer. "Efficient and Stable Arnoldi Restarts for Matrix Functions Based on Quadrature". In: SIAM Journal on Matrix Analysis and Applications 35.2 (2014), pp. 661-683 (cit. on p. 4).
[25] Andreas Frommer, Karsten Kahl, Th Lippert, and Hannah Rittich. "2-Norm Error Bounds and Estimates for Lanczos Approximations Linear Systems and Rational Matrix Functions". In: SIAM Journal on Matrix Analysis and Applications 34.3 (2013), pp. 1046-1065, eprint: https://doi.org/
1137/ / 110859749 (cit. on p. 4).
[26] Andreas Frommer, Karsten Kahl, Marcel Schweitzer, and Manuel Tsolakis. "Krylov Sub- space Restarting for Matrix Laplace Transforms", In: SIAM Journal on Matrix Analysis and Applications 44.2 (2023), pp. 693-717 (cit. on p. 2)
[27] Andreas Frommer and Marcel Schweitzer. "Error bounds and estimates for Krylov subspace approximations of Stieltjes matrix functions". In: BIT Numerical Mathematics 56.3 (2015), pp. 865-892 (cit. p. 4).
[28] Andreas Frommer and Valeria Simoncini. "Error Bounds for Lanczos Approximations of Ra- tional Functions of Matrices", In: Numerical Validation in Current Hardware Architectures. 2009, pp. 203-216 (cit. on p. 4).
[29] Roy Frostig, Cameron Musco, Christopher Musco, and Aaron Sidford, "Principal component projection without principal component analysis". In: Proceedings of the 33nd International Conference on Machine Learning (ICML 2016). 2016, pp. 2349-2357 (cit. on pp. 1, 9).
[30] E. Gallopoulos and Y. Saad, "Efficient Solution of Parabolic Equations by Krylov Approx- imation Methods". In: SIAM Journal on Scientific and Statistical Computing 13.5 (1992),
pp. 1264-1264 (cit. on p. 1).
[31] Jacob Gardner, Geoff Pleiss, Kilian Q Weinberger, David Bindel, and Andrew G Wilson. "GPyTorch: Blackbox Matrix-Matrix Gaussian Process Inference with GPU Acceleration". In: Advances in Neural Information Processing Systems (NeurIPS 2018). Vol. 31, 2018 (cit.
on p. 1).
[32] Behrooz Ghorbani, Shankar Krishnan, and Ying Xiao. "An Investigation into Neural Net Op- timization via Hessian Eigenvalue Density".. In: Proceedings of the 36th International Con- ference on Machine Learning (ICML 2019). Vol. 97. 2019, pp. 2232-2241 (cit. on pp. 1, 21).
[33] Anne Greenbaum. "Behavior of slightly perturbed Lanczos and conjugate-gradien recur- rences". In: Linear Algebra and its Applications 113 (1989), pp. 7-63 (cit. on pp. 2, 10).
12

---

[34] Anne Greenbaum. "Comparison of splittings used with the conjugate gradient algorithm", In: Numerische Mathematik 33.2 (1979), pp. 181-193 (cit. on p. 20).
[35] Anne Greenbaum. Iterative Methods for Solving Linear Systems. Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 1997 (cit. on p. 3).
[36] Stefan Güttel and Marcel Schweitzer. "A Comparison of Limited- memoryy Krylov Methods for Stieltjes Functions of Hermitian Matrices". In: SIAM Journal Matrix Analysis and Applications 42.1 (2021), pp. 83-107 (cit. on p. 9).
[37] Nicholas Hale, Nicholas J. Higham, and Lloyd N. Trefethen. "Computing A", log and Related Matrix Functions by Contour Integrals". In: SIAM Journal on Numerical Analysis 46.5 (2008), pp. 2505-2523. eprint: https: / doi.. org 10. 1137 070700607 (cit. on
pp. 4, 6).
[38]
Nathan Halko, Per-Gunnar Martinsson, and Joel A. Tropp. "Finding structure with random- ness: probabilistic algorithms for constructing approximate matrix decompositions", In: SIAM Review 53.2 (2011), pp. 217-288 (cit. on 2).
[39] Insu Han, Dmitry Malioutov, Haim Avron, and Jinwoo Shin. "Approximating Spectral Sums of Large-Scale Matrices using Stochastic Chebyshev Approximations". In: SIAM Journal on Scientific Computing 39.4 (2017), A1558-A1585 (cit. on p. 1).
[40] Clemens Hofreither. "An Algorithm for Best Rational Approximation Based on Barycentric Rational Interpolation". In: Numerical Algorithms 88.1 (2021), pp. 365-388 (cit. pp. 8, 9).
[41] Marija D. Ilce, Ian W. Turner, and Daniel P. Simpson. "A restarted Lanczos approximation to functions of symmetric matrix". In: IMA Journal of Numerical Analysis 30.4 (2009), pp. 1064-1061 (cit. on p. 4).
[42] Charles H Janson. "Measuring evolutionary constraints: a Markov model for phylogenetic transitions among seed dispersal syndromes", In: Evolution 46.1 (1992), pp. 136-158 (cit. on
p. 1).
[43] Robert A Jarrow, David Lando, and Stuart M Turnbull. "A Markov model for the term struc- ture of credit risk spreads". In: The Review of Financial Studies 10.2 (1997), pp. 481-523 (cit. on p. 1).
[44] Yujia Jin and Aaron Sidford. "Principal Component Projection and Regression in Nearly Lin- ear Time through Asymmetric SVRG". In: Advances in Neural Information Processing Sys- tems (NeurIPS 2019). Vol. 32. 2019 (cit. on pp. 1, 2, 9, 27).
[45] Sadegh Jokar and Bahman Mehri. "The best approximation of some rational functions in uniform norm". In: Applied Numerical Mathematics 55.2 (2005), pp. 204-214 (cit. on p. 18).
[46] Cornelius Lanczos. "An iteration method for the solution of the eigenvalue problem of lin- ear differential and integral operators". In: Journal of Research of the National Bureau of Standards (1950), pp. 255-282 (cit. on onp.2 2).
[47] Jörg Liesen and Zdenčk Strakoš. Krylov subspace methods: principles and analysis. 1st. Nu- merical mathematics and scientific computation. Oxford University Press, 2013 (cit. on pp. 2,
3, 18).
[48] Stefano Massei and Francesco Tudisco. "Optimizing network robustness via Krylov sub- spaces". In: ESAIM: Mathematical Modelling and Numerical Analysis 58.1 (2024). pp. 131- 155 (cit. on p. 1).
[49] Gérard Meurant and Jurjen Duintjer Tebbens. Krylov Methods for Nonsymmetric Linear Sys- tems: From Theory to Computations. Springer International Publishing. 2020 (cit. on p. 18).
[50] Gérard Meurant and Zdenčk Strakoš. "The Lanczos and conjugate gradient algorithms in finite precision arithmetic". In: Acta Numerica 15 (2006), pp. 471-542 (cit. on p. 10).
[51] Fredrik Johansson et al. mpmath: Python library for arbiraayyprecisson frac--rirccfllloatiiii floating point arith- metic (version 0. 18). /ttpa://mpmtth... org/ 2013 (cit. on p. 7).
[52] Cameron Musco, Christopher Musco, and Aaron Sidford. "Stability of the Lanezos Method for Matrix Function Approximation". In: Proceedings of the Twenty-Ninth Annual ACM-
SIAM Symposium on Discrete Algorithms (SODA 2018). 2018 pp. 1605-1624 (cit. on pp. 2,
3, 10).
13

---

[53] Cameron Musco, Praneeth Netrapalli, Aaron Sidford, Shashanka Ubaru, and David P. Woodruff. "Spectrum Approximation Beyond Fast Matrix Multiplication: Algorithms and Hardness". In: 9th Innovations in Theoretical Computer Science Conference (ITCS 2018). 2018 (cit. on p. 9).
[54] Lorenzo Orecchia, Sushant Sachdeva, and Nisheeth K. Vishnoi. "Approximating the expo- nential, the Lanczos method and an O(m)-time spectral algorithm for balanced separator". In: Proceedings of the 44th symposium on Theory of Computing (STOC 2012). 2012 (cit. on p. 3).
[55] Christopher Conway Paige. "Accuracy and effectiveness of the Lanczos algorithm for the symmetric eigenproblem". In: Linear Algebra and its Applications 34 (1980), pp. 235-258 (cit. on p. 10).
[56] Christopher Conway Paige. "Error Analysis of the Lanczos Algorithm for Tridiagonalizing a Symmetric Matrix". In: IMA Journal of Applied Mathematics 18.3 (1976), pp. 341-349 (cit.
on p. 10).
[57] Christopher Conway Paige. "The computation of eigenvalues and eigenvectors of very large sparse matrices." PhD thesis. University of London, 1971 (cit. on p. 10).
[58]
Vardan Papyan. The Full Spectrum of Deepnet Hessians at Scale: Dynamics with SGD Train- ing and Sample Size. 2019, arXiv: ..00620722 [cs LG] (cit. on pp. 1, 21).
[59] Geoff Pleiss, Martin Jankowiak, David Eriksson, Anil Damle, and Jacob R. Gardner. "Fast matrix square roots with applications to gaussian processes and Bayesian optimization". In: Proceedings of the 34th International Conference on Neural Information Processing Systems (NeurIPS 2020). 2020 (cit. on pp. 1, 2, 6, 9, 21).
[60] Theodore Rivlin. An introduction to the approximation of functions. Dover books on ad- vanced mathematics. Dover, 1981 (cit. on p. 20).
[61] Yousef Saad. "Analysis of Some Krylov Subspace Approximations to the Matrix Exponential Operator". In: SIAM Journal on Numerical Analysis 29.1 (1992), pp. 209-228. eprint: https: / //doi. org/ 10. 1137 00729014 14 (cit. on pp. 3, 15, 18).
[62] Yousef Saad. Iterative Methods for Sparse Linear Systems. Society for Industrial and Applied Mathematics, 2003 (cit. on p. 2).
[63] E. B. Saff, A. Schönhage, and R. S. Varga. "Geometric convergence to e by rational func- tions with real poles". In: Numerische Mathematik 25.3 (1975), pp. 307-322 (cit. on pp. 6, 8).
[64] Lloyd N. Trefethen. Approximation Theory and Approximation Practice, Extended Edition. Society for Industrial and Applied Mathematics, 2019 (cit. on pp. 6, 20).
[65] Joel A. Tropp and Robert J. Webber. Randomized algorithms for low-rank matrix approxi- mation: Design, analysis, and applications. 2023. arXiv: 2306 12418 [math.NA NA] (cit. on
p. 2).
[66] Roman Vershynin. Tligh-Dimensional Probubility: An Introduction with Applications in Data Science. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge Univer- sity Press, 2018 (cit. p. 21).
[67] Sheng Wang. Yuan Sun, Christopher Musco, and Zhifeng Bao. "Public Transport Planning: When Transit Network Connectivity Meets Commuting Demand". In: Proceedings of the 202/ International Conference on Management of Data. SIGMOD '21. 2021, pp. 1906-1919 (cit. on p. 1).
[68] Frederick V Waugh and Martin E Abel. "On fractional powers of matrix". In: Journal of the American Statistical Association 62.319 (1967), pp. 1018-1021 (cit. on p. 1).
[69] Zhewei Yao, Amir Gholami, Kurt Keutzer, and Michael Mahoney. PyHessian: Neural Net- works Through the Lens of the Hessian. 2020. arXiv: 1912. 07145 [cs. LG] (cit. pp. 1, 21).
14

---

A Proof of Theorem 4
A.1 Notation
We first introduce notation used throughout. Given rational function r(x) n(x)/m(x) with numerator degree p denominator degree q, for < sisj q, we define:
mi, (x) == 1
(x x zR).
R=1
We adopt the convention that m )+1,(xx) == 1. Note that mi.j is (j - -i 1)-degree polynomial. Define also
ry(x) n(x) m(x)//1.,(x))
Note that m1.4 = 80 rg r(x). Recall that for any function fx), the Lanczos-FA approximation ((A) is (f:,( Q Define
err, (f) f(A)b lan (f: b).
Proving Theorem 4 amounts proving upper bound |lerr=(r))l Finally, for any symmetric positive definite matrix M, define
opt, (; == argmin (If(A)b x||M- XE (A.b)
For brevity, in the analysis below will write lan (f) and opt, (f.M) in place of lanx (f: b) and opt, (f. M; A, b), since A and are fixed throughout the analysis.
A.2 Simplifying lane(r) (ry) and opt, (ry, Aj)
We begin with few standard results that will be useful for the proof of Theorem 4.
Lemma 8 ([22, 61]). For a polynomial n(x), (n) n(A)b if k deg(n).
Proof. By definition the Krylov subspace, A/b Kg (A, for all k. Since QQ is the 2- norm orthogonal projector onto Kg (A, b), then QQ A/b A/b for all k. Iteratively applying this fact,
A/b QQ A/b AA / AA
= 20 AQQ A A AQ ab
Using that QTAQ = T, we find that A/b QT and thus, by linearity, n(A)b Qu(T)Q lanx(n),
Lemma 9. Let A₁ be in (5). Then, for any opt, (ry, A,) Q(T -
Proof. The A -nor projector onto K2 (A, b) is 1(0 so
opt. (ry, A,) 1100
Recall that A, #(A-2)1), Since TA and Q== we have ("==,,( =(T 1)) Noting also that Ayr, (A) (A) and plugging into (9) gives the result.
A.3 telescoping sum for the error in terms of opt, Ar)
Our first main result is decomposition for the Lanczos-Fa emror
Lemma 10. Suppose & deg(n). Then,
Σ
i' i=j+1
9 Q(T 2, -
(r)( (r,
r(A)b - labb-lan (r)
where we adopt the convention that 17
B, for any set of matrices (B.).
15

---

Proof. We can decompose err; (T))
(rr) =r)(A)b Qr -
[r)( opt, (r) -0000)0))0000)))))00000000000000
(10)
Focusing on the second term and using Lemma 9,
opt (rj, Aj)
Qr,(T)Q
Q(T Q 1)=1(A)b-
Q(T
2, z,
e
[err((;-11]
Substituting this into (10), we have
(rr)(r)) = (A)b opt, ()), 2 (
(11)
Next, notice that opt, (/1, A) is exactly the Lanczos approximation to r1 (A)b. Indeed, since QTQ I, Lemmas 8 and 9 imply
lan Q(T - -
QQTT Q
Q Q nA)b opt, A1).
Using this base case, repeatedly apply our decomposition of err, (ry) in (11) obtain
err (r) = (Tg) Σ
QT O (r,(A) - opt, (1)(A)b-opt((r)) (r,A,)
1=1+
A.4 Bounding each term the telescoping sum and combining
The following lemma allows us to relate the optimality of the functions rj measured in the Ay: norm that of rg = as measured in the 2-norm. Lemma 11. For any k > -.. j. [A)j opt, (r.,A))l
x(A)) 1/2 [IM]+llg(A)ll (e((p)<k((g-)))
min
||[(A)b p(A)bll2.
Proof. Recall that opt, (r, Aj) is the optimal approximation to r, (A)b b in the A -norm, can be related the optimal approximation the 2-norm.
(A) opt, (r,, Al2 |A, A)
1/ (r) opt. (r)
S ||A] 1/2 |Iry (A)b opt, (r) (, Ay)IA
= ||A] 11A//11/2 deg(p)<k
min
[|r] p(A)blla,
< ||A, 11 1/2 || 1/2 min
[A) ((Abll22
deg(p)<c
= x(A min deg(p)<k
(A) p(A)bll.
We now relate the error of approximating (A) to that of approximating (A)b r(A)b: (A)b p(A)bll2
min deg(p)<c
min deg(p)<k (MM(LLg(A)r (A)b p(A)bl22
min [| (A) (A)b m m+ m (A)b
dog(p)<k-19-1)
$ [llmj+l.q deg(p)ck-(g-1)
min
[|r(A]b p(A)bll2
16

---

Combining these two steps proves the lemma.
With the above results place, we can prove Theorem We will apply the triangle inequality to the telescoping sum for err, (r) (Lemma 10) and bound each term by the error of the optimal approximation r r(A)b (Lemma 11).
Proof of Theorem 4. Focus on a single term in the sum in Lemma 10. First we will get rid of Q and Q For $ j < 4, using that Q = I we have that
Q(T Q(T-==))-- Q = (T) (T) s 1 (T) (T)
(12)
inf+1
Also note that by the convention adopted in Lemma 10, for q. 11 (xT=2xx)--0 =1= m (T). Next we claim that
[|mj+11q (T) z [11 112-
(13)
To see why this is the case, note that because T ATAQ, the eigenvalues of T are contained in
= [mmin (A), Amax (A)]. By assumption, the roots of mi1,, (x) are all real and lie outside of J. Since there is at most one critical point between distinct roots, there can only be one critical point of
m | (x) in I. Thus, there can be no local minima of [m)+1.g(x)) in the interior of I, Rather, the Amax-
minimum of [m,+1.g(2) (x) over I must be attained at the boundary; i.. at Amin
We now apply the triangle inequality to the telescoping sum of Lemma 10:
(r)lr Σ
q Q(T 21) -
||r((A)b (r,
E++1
Applying (12) and (13) we then have
(err(()) &
Ilm (1))+1((A) (A) ||)( -ptt
Finally, using Lemma 11 we find
min [r)l2 12 ((m)+1.(())) deg(p)<k-(q-))
[lp(A)b r(A)bll2. (14)
We can simplify by noting that
min (kkp))k-((-J))
min IIp(A)b- (k(p))kk(g--))
Ilp(A)b
and we can combine all the condition number factors using
Σ
2 1
x(A)) 12 K(mj+1.g (A)) s
x(A))
<(A,) S x(Ai).
1=1
1=1
Ej+1
i=1
The result follows by plugging into (14).
A.5 Proof of Lemma
We can bound the error of Lanczos-FA on f(x) using triangle inequality follows:
|l((A) - Ian (/:,()AA (A) (aaf/((Abrr(a ( (AA)bl2 [|r(A)) lan b)
+Illan (r:; A, b) lan A,
(15)
The first and third terms of (15) are controlled by the maximum error of approximating f(x) with r(x) over I. Specifically, if we let ||r - filr == max [r(x) f(x).
IlS(A)b r(A)bl22 |lb]l (I((A) r(A)ll2 ||b]] ||r fll
17

---

and similarly, using that A(T) C I.
| (r: (fa A. b)ll The second term of (15) controlled using (7) and a triangle incquality:
(f |lbll2 ||r fll..
[|r(A)b lan; Ar; s C min deg(p)<o,k
[|r(A)b p(A)bll2
s C, deg(p)<c,k
min
(IIr(A)b- f(A)bll2 (IS(A)) ((A)bll2)
- C ||r fll + Cr min
II((A)b - p(A)bll2.
deg(p)<<ck
Combining, we obtain the bound
(If( - lan, (f: b)ll2
s min (C, + 2) ||r- fllz Cr min
II/(A)) 1100--((Abbii))) p(A)bllz).
(16)
deglp)<c,k
B Comparison to Prior Work
B.1 Details of existing near-optimality guarantees for Lanczos-FA
In this section, we review in detail the prior analyses of Lanczos-FA cited in Section 1.3. These are the only near-optimality guarantees for Lanczos-FA of which we are aware. Instance optimal- ity trivially holds when f(x) is a polynomial of degree k. In this case, it is well known that lany (f: A, b) exactly applies f(A)b; i.c. (If lan; (f:A, b)ll2 = 0 [61]. When f() 1/x and A is positive definite, the Lanczos-FA -FA algorithm is mathematically equivalent the well-known conjugate gradient algorithm for solving a system Ax = b. This implies that Lanczos-FA is the op- timal approximation in the Krylov subspace with respect the A-norm; that is,
lan, (1/x; A, b) = p(A)b,
p = argmin ||A ||A p(A)bllA-
dc
This immediately yields near instance optimality in the Euclidean norm:
|| - lan, (1/
vk(A) min
||A p(A)bll2. deg(p)<k
(17)
where x(A) is the condition number Amax/Amin of A. That is, Definition 2 is satisfied with C vk(A) and 1.
Note that the best polynomial approximation to 1/x on I [mmin mmxx] already converges at geometric rate [45]:
8r1
min max 11/x - p(x) [= (12 1 (A Amin)
deg(p)<ck xex
1 Vx(A)
That is, Fact suffices to prove the exponential convergence of the conjugate gradient method. However, conjugate gradient (and equivalently Lanczos-FA) often converges even faster, which is explained theoretically by the stronger instance- and spectrum -optttulmpppiialityy guarantees that the method satisfies [6, 47, 8).
If A is positive definite, T may have eigenvalue at FA approximation to "b can be arbitrarily large. In fact, the same is true for any function f(x) which is much larger on I [Amin Amax] than on A. the of eigenvalues of A. However, while the Lanczos-FA FA iterates may be bad at some iterations, the overall convergence of Lanczos-FA with f(x) 1/x is actually good [12]. In particular, the convergence of Lanczos-FA can be related the MINRES algorithm [15, 49], which produces an optimal approximation to A - with respect to the A -norm. This allows certain optimality guarantees for MINRES be transferred Lanczos-Fit even on indefinite systems. In particular, [12] asserts that for every k, there exists k < k such that'
and the error of the Lanczos-
|| lb - lanx (1/x; A, b) s v(A) min deg(p)<<k
||A p(A)bll2.
(18)
For indefinite matrices, x(A) denotes the ratio of the largest eigenvalue magnitude to the smallest eigen- value magnitude.
18

---

k), we say that a Krylov While (18) does not quite fit the form of Definition 2 because of the k* on the left side, it is similar in spirit. Also note that the dependence on k in the prefactor V + is not of great significance. Indeed, the convergence of the optimal polynomial approximation to 1/x on two intervals bounded away from zero is geometric.
Finally, guarantee for the matrix exponential is proved in [20, (45)). They show that the Lanczos- FA iterate satisfies the guarantee
exp(-1A)b (epp(ex-- b)ll2
s 3||A||
max 0gsar
min deg(p)<k-2
||ep(=ss)) (0000-((A)bi) p(A)bll2)
Again, this bound does not quite fit Definition due to the maximization over The authors of [20] state, "It is not known whether the maximum over E [0,r] in (45) can be omitted and set equal in the right-hand side."
B.2 Comparison to Lanczos-OR
In [11] an algorithm called the Lanezos method for optimal rational function approximation (Lanczos-OR) was developed. For rational functions of the form (4), when run for k > deg(n) iterations, Lanczos-OR produces iterates
lrnOR == argmin [|{(A) (()))()))
x
This implies a 2-norm instance-optimality guarantee (Definition 2):
[|r)) - I A, b)ll2 S vk(r(A)) min deg(p)<k-l9/2)
((A)b=p(A)bll>>
(19)
This guarantee is similar to and often somewhat better than Theorem 4. For example, when r(x) 1/m(x), (x(r(A)) vx(A,) x(A))...x(A)) x(Ag). In this case (19) improves on Theorem 4 by factor and a square root. However, the bound for a different algorithm; it does not extend to the ubiquitous Lanczos-FA method.
Moreover, despite this bound, Lanczos-OP OR usually performs worse than Lanezos-FA in practice enro is measured in the 2-norm. It is not hard to find examples where (19) is essentially sharp: see for instance Figure 8, On the other hand, we have been unable to find any numerical examples where Theorem 4 is sharp, suggesting that it may be possible to prove tighter bound for Lanczos-FA.
C Instance, spectrum, and FOV optimality
The convergence of Lanczos-Fil entirely determined by spectral properties of A. Therefore, fol- lowing the classical analysis of f(x) = 1/x, we can relax the notion of instance optimality (Defini- tion 2) by removing dependence on b:
Definition (Near Spectrum Optimality). For an input instance (f.A. k), we say that a Krylov method is nearly spectrum optimal with parameters C, c.c.l if
(lf(A)b alga (f:AAfAA)ll. Cl[blI2
min deg(p)<ck
IIf(A) (PA))-(A))l2.
Note that Definition 12 depends only the spectrum of A, and not the interaction of the eigen- vectors of A with the starting vector Hence, we call it "spectrum optimality". Note also that the guarantec of Theorem 6 fits this form. We can further relax this notion of optimality by weakening its dependence on A:
Definition 13 (Near FOV Optimality). For an input instance method is nearly FOV optimal with parameters C. if
1f(x) (f(x)-p(x)) p(s))
[A - (aa))ball,.4 (f:A, s Cl|bil2
min
max
deg(p) <ok (
19

---

Note that Fact 3 fits this form. Each type of optimality tightly relaxes the previous one. Definition implies Definition 12 and, as we show in Theorem 14, for any f. and A, there exists "worst-case" choice of b for which their bounds coincide; likewise, Definition 12 implies Definition 13, and for any fxx) and there exists a "worst-case" choice of and b for which their bounds coincide. Theorem 14.
1. Given any d, discrete set A of most d real points, f(x), and k < d, there exists symmetric matrix BBddd and vector E Rd such that the spectrum of A is A and min [If((Abb (000-p((Abl///llll)l)) deglp)<k
min max ( deg(p)<k XEA
2. Given any d, Amin, Amax, f(x) continuous on [Amins Amax], and k d, there exists sym- metric matrix Raxd and vector b Bd such that the smallest eigenvalue of A is Amin- the largest eigenvalue of A is max, and
min IIf(A)b p(A))ll2/(lblle> dk
min max If(x) (x)xxxx)) p(x)l.
deg(p)<k xe1
Our proof follows the approach of [34, Theorem I] closely.
Proof. For Item note that if |A] k, then both sides of the expression to consider the case |A| 2 k+ + 1. On a discrete set A with at least & + distinct points, f(xx has best (on A) polynomial approximation p" of degree k that equioscillates f(x) at k points [60, Theorem 1.11]. That is, there exist values 21 < < Ax+1 contained in A such that
Thus, it suffices
f(A) - p" (2) = (-1)
i=
where
: max [f(x) px()x
Similarly, for Item if fx) is continuous on [Aminn Amax], then f(x has a unique best poly- nomial approximation p' of degree k - which equioseillates at k points dett in I (see for instance [60, 64]), and we analogously define e as
E = max [f(x p (x)]
xeX
Set
A diag (1, ... det, d 4+1). [bi..... b 0, OJ.
d-k-1 times
d-k-1 times
Then,
k+1 B2U(()))
k-1
min deg(p)<c
(l((A)b-p(A)bl))
min agg. mg-1
Σ
p()) = Eajal
1=0
We would like to find b so that f(1:) p(d:) (-1) i.e. so that p == p". When b fixed, for each j 0, 1,....
the solution to this least squares problem must satisfy
d
- (1)))-((a))=-22) Σ (F())) p(2.)
0
....k....kk
dæz 1=1
When p p", this gives the conditions
0= 0= ----)--
L......kk
1=1
Without loss of generality, we can assume linear system
||b||2 so that b7 ..... b2
Thus, we obtain
1
62
-2e
2ed9
2(-1) k+1 10
by
E+
-22
201h-1
2(-1) k+l 2
b2
k+1
20

---

We can rewrite this system as
1
-1
(-1)*)
22
12
b
R+1
14-
1--1
1--1
(-1)
k+
This system can be solved analytically via Cramer's ' rule and has solution
k+1
(2
bi
a1 2n)
(1,
10)
1=+1 is!
m==l 1=1 j=1+1
jut
iim jum
Clearly (Aj- - Ai) is positive for j E (i l.... ..... + 1}. Thus the numerator and denominator above positive, the entries of b well-defined and real. This implies that, for these choices of A and b,
k+1 by:²²
min (If(A)b- ((A)bII dog(p)<k
min am.... 4.g-
((4))))((4)))--
--pee 2.
Taking the square root of both sides gives the result.
C.1 Relation between spectrum and instance optimality for random vectors
We also note that, with some additional assumptions on b, spectrum optimality implies instance optimality:
Lemma 15. Lei (ui..... Ud} be the eigenvectors of A. Given a problem instance (S, k), if an algorithm is nearly spectrum optimal with parameters C and then it is nearly instance optimal
with parameters C [|b]] and
min; m
The bound is useful e.g., when b is a random vector. For example, when has independent and iden- [bl],
tically distributed Gaussian entries, min, bu]b] Odd3/2) with high probability. Such random vectors
occur when sampling Gaussians [59] and are extremely common in machine learning applications related to trace and spectrum approximation [58, 32, 69, 7, 13, 14].
Let V be the matrix whose ith column is V/. Let w/ v.b and w = v be the vector whose ith
entry is wy. Then,
[If(A) P(A)bI -((2)-
> min Σ ((())
1=1
Clearly (f(a)) ( max,ea 1f(x) pxx] substituting into Definition 12 proves the result.
Taking the square root of both sides and
For completeness, now show that if b (O., Lz), then with high probability,
|lbll
= Od)33))
min, by]b]
First, with high probability, [[b]]2 = O(V)); see, for example, [66, p. 3.1.2]. Second, we show that with high probability |v] ((11d) for all j. Since is orthonormal, = vTb (0, L4). By anti-concentration of the normal distribution, the probability that [wyl > e is at least 0.4e for any By union bound, this holds simultaneously for all with probability at least - 0.4ed. Setting = O(1/d finishes the argument. Finally, by another union bound, our bounds on the numerator and the denominator hold simultancously with high probability.
21

---

D Proofs of Theorem 6 and Theorem
We begin by quoting bounds from [10]:
Lemma 16. For all 2 1, the Lanczos-Pa FA iterate satisfies the bounds
A
1. | 1aa-b-landde (x (x
||A- llan
2R3/2
Amax
2. ||A 1 blan (x-1224 b)ll2 v
πk ||A
lan (1/x; b)ll2.
Proof. For Part 1, see Example 4.1 in [10]. For part Part 2, we use the same proof again but sub- stituting the inverse square root for the square root. For both parts, we have simplified the resulting bounds using the fact that for k 1.
r(k 1/2) [(k 1)
V
r(k ((2)1/2) T(k +
$
S
k3/2
VK
Next, we bound the error of Lanczos-FA on linear systems by the error of the optimal polynomial approximation to x 1/2, Lemma 17. For all 1, the Lanczos-Fa FA iterate satisfies the bounds
3 v<(A)||blz
max xEA Vxx p(.)
| b (1/ s
min (pg(p)<kp2
Va
Proof. Let
p'(x) argmin max - p(x) dk/(p)<</2 XEA VX
and note that p2x)2 is polynomial of degree less than k. Therefore,
| p(A)bl22 = |lbll2 min max p(x) s [lblla max ! x
p (x)2
min dc
deg p<ck x x
x
Thus, since 11/x-p'(xx 11/va+p'<>) 11/V5 (x)) we plug the above into (17) to
obtain a bound
p (()))
- -0---((())
1
P'(o))
| IA- = |x(A)||bll
max xEλ
max xea
Vx
Vx
Now use the optimality of pr'ppxxxxx (x)in approximating 1/ 1/v on A to bound:
2
1
max
+x((x)
max
+ p px)
VX
xEA
Vx
Vx
S max
max
p p'(x)
xεA
Vx
3
max
0(x)
5
Va
X
vr
Vd
Substituting this inequality above proves the lemma.
The main proofs now follow directly.
Proof of Theorem 6. Substitute Lemma 17 into Lemma 16, part (b).
22

---

10a
10
10
10-13
10-11
10
10-19
42
10:21
10"54
10:36
Fast
10-4
Taconem 6 (lef)
Tacopem 7 (right)
10
10-51
LanezcoFAA
10
Opiima)
10
20
30
40
50
60
Namber of iterations (K)
Figure 7: The bounds of Theorems and capture the convergence behavior of Lanczos-FA for b. In particular, they can be tighter than the FOV optimality bound of Fact 3. The predicted rate of convergence is about half that observed for Lanczos-FA this example; this is due to the parameter 1/2 in both bounds.
Proof of Theorem 7. Let
p(x) p'(x) (0 (O)))
p (x) =
argmin dkg(p)<</2+1 xEAU(D)
max | - p(x)
and note that p polynomial of degree less than k/2. Then
p(o)
min (k)
max XEA
max
p(x)
VX
XEA
VF
max x
(11/xl - [xxx((xx]]
max - - p0(0)
Anin
max (lVi ([VT-p'(0r p0(x) -v0-p'(0))
Amin xEλ
max Amin XEAU(O)
- px)]
Combining Lemma 16 part (a), Lemma 17, and the above proves the theorem.
E
Additional Experiments
E.1 Validating Theorems and 7
To understand Theorems 6 and 7, we repeat the experiment using the inverse square root and square root functions. We reuse the second and third matrix and b vectors from the previous experiment. As Figure shows, our bounds are tighter than Fact 3. They closely resemble the true convergence of Lanczos-FA, but stretched horizontally by factor minimization in both bounds is k/2, while in these examples Lanczos-FA performs nearly instance optimally (that is, as well the best degree polynomial). Still, our bounds capture the correct qualitative behavior of the algorithm.
This is because the degree of the polynomial
E.2 Rational functions with poles between eigenvalues
In Theorem 4, we assume that the poles of the rational function r(x) lie outside the interval I [Amin ] containing A eigenvalues. If the poles of r(x) live in I, then the Lanczos-FA iterate
23

---

10:0
10:43
10
10
Optimalty Rutio
10:35
10
10
10
10FT
10t5
10ss
10
10' 10 10" 10P
The maximum observed ratio between the error of Lanczos-OR and the optimal error choice of right hand side b when approximating A for matrices with varying condition number Each point shows the optimality ratio for a different pair of к and q. Points with the same color correspond to the same value of K. On the left, the dotted line plots g(q) x for the maximum considered (10%). On the right, the dotted line plots g(x) x4/2 for the maximum considered (26 = 64). Overall, the optimality ratio appears to grow as Q(x4/2), Contrast with Figure 3.
Figure 8:
can arbitrarily far from the optimal iterate. Indeed, an eigenvalue of T might be very close to pole, causing S(T) to be poorly behaved. This behavior is easily observed e.g. for f(A) A and rules out a direct analog of Theorem when r(x) has poles in I. However, as noted in (18), for 1/x there exists bounds for Lanczos-FA FA in terms of the best possible KSM [12].
However, as discussed Appendix B.1, the "overall" convergence of Lanczos-FA still seems to closely follow the optimal approximation for many matrix functions. Specifically, even if the Lanczos-FA iterate deviates significantly from optimal on some iterations, there are other iterations where it is very close to optimal. We illustrate this phenomena in Figure 6 for two rational func- tions with poles in I = [Amin- Amax], as well as the function sign(x), which has discontinuity in I and would typically be approximated by rational function with conjugate pairs of imaginary poles on the imaginary axis which become increasingly close to the the interval as the rational function degree increases."
E.3 Comparison with Lanczos-OR
As discussed in Appendix B.2, the Lanezos-On OR method is an alternative method for rational matrix function approximation for which a tighter guarantee than Theorem 4 is currently known (the iter- ates are optimal in a certain norm). However, Lanczos-FA usually outperforms Lanezos-OR when measured in the 2-norm. For example, as discussed in Section 4.1, when approximating A bb,,, the worst optimality ratio of Lanczos-FA that were able observe was O(vq x(A)), much lower than Theorem& 4 would suggest). Repeating the same experiment with Lanczos-OR shows that the method's optimality ratio appears to grow as ((((A)()) (see Figure 8).
We remark that [11] shows that applying their Lanczos-O0 OR method to each term in the partial fraction decomposition of the approximating rational function yields an approximation to the sign function that seems avoid the oscillations of Lanczos-FA.
24

---

NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: Our abstract refers the main theoretical result, Theorem 4, and the ex- perimental results of Section 4, particularly Figure 2.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made in the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question:
Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: Each of our theorems clearly states all assumptions. Section 5 discusses the limitations of our current results and the open questions that remain.
Guidelines:
The answer NA means that the paper has limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create separate "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The au- thors should reflect on how these assumptions might be violated practice and what the implications would be.
The authors should reflect the scope of the claims made, e.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the ap- proach. For example, facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech- speech-to-tex system might not be used reliably provide closed captions for online lectures because fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to ad- dress problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community, Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory Assumptions and Proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
25

---

Answer:
[Yes]
Justification: Formal theorem statements include all assumptions. All proofs appear in Appendices A, C, and D.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced in the statement of any theo- rems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be comple- mented by formal proofs provided in appendix supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental Result Reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main experimental results of the paper to the extent that it affects the main claims and/or conclu- sions of the paper (regardless of whether the code and data provided not)?
Answer: [Yes]
Justification: We release our code in the supplementary material and on GitHub. In addi- tion, details of the experiments are provided in the text.
Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided not.
If the contribution is a dataset and/or model, the authors should describe the steps taken make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is a specific model and empirical evaluation, may be necessary to either make it possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how replicate the results, access a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate the research performed.
While NeurlPS does not require releasing code, the conference does require all sub- missions to provide some reasonable avenue for reproducibility, which may depend on the nature the contribution. For example
(a) the contribution primarily a new algorithm, the paper should make it clear how to reproduce that algorithm.
(b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results way produce the model (e.g., with an open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case au- thors are welcome to describe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should possible for other researchers to have some path to reproducing or verifying the results.
26

---

5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, described in supplemental material?
Answer: [Yes]
Justification: We provide our code in the supplementary material and on GitHub. The code to produce Figure is not released, as it was shared with us privately by the authors of the paper in which it first appeared [44]. Standard Python libraries are used for reproducibilityy and a README is provided with instructions for how to generate the figures.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (attps://nips) cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https:
cttpsccccpbbbcicc httppppullllupeeeeeeeeppppppaeae guw.es/ /eaidddd///oodbbnnne CooeeSubeeiieeeecc ssi onPolicy) for more details.
The authors should provide instructions on data access and preparation, including how access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing much information possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental Setting/Details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: These are described in the text and the code included in the supplementary materials.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment Statistical Significance
Question: Does the paper report error bars suitably and correctly defined or other appropri- ate information about the statistical significance of the experiments?
Answer:
[No]
Justification:
The algorithms we test are random and no sampling is performed.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
27

---

The factors of variability that the error bars are capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall run with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
is OK to report I-sigma error bars, but one should state it. The authors should prefer- ably report a 2-sigma error bar than state that they have a 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful not show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain in the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments Compute Resources Question:
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory. time of execution) needed reproduce the experiments?
Answer: [No]
Justification: No special computing resources used. All experiments were performed on a standard laptop in a few minutes.
Guidelines:
The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount of compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported the paper (e.g.. preliminary or failed experiments that didn't t make it into the paper).
9. Code Of Ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics https:
Answer:
[Yes]
Justification:
No human subjects data was used. No clear potential exists for harmful consequences.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader Impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [No]
Justification: This paper provides new theoretical analysis of an old and widely used al- gorithm, the Lanczos method. The goal to achieve better scientific understanding of the method. The Lanczos method has been used in huge variety of applications across machine learning and scientific computing. and thus has societal impact, but the results and focus of this paper do not directly relate this impact.
28

---

Guidelines:
The answer NA means that there is no societal impact of the work performed.
If the authors answer NA or No. they should explain why their work has no societal impact or why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact spe- cific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied partieular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point that improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitiga- tion strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of MI.).
11. Safeguards
Does the paper deseribe safeguards that have been put in place for responsible release of data or models that have high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
Question:
Answer: [NA]
Justification: See above.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by re- quiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification:
We cite third party software libraries and use them properly.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package or dataset.
The authors should state which version of the asset used and, if possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
29

---

For scraped data from a particular source (e.g., website), the copyright and terms of service of that souree should be provided.
assets are released, the license, copyright information, and terms of use the pack- age should be provided. For popular datasets, paperswi tacodhcodee.. com/ datasets has curated licenses for some datasets. Their licensing guide can help determine the li- cense of dataset.
For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided.
If this information is not available online, the authors encouraged to reach out to the asset's ' creators.
13. New Assets
Question: Are new assets introduced in the paper well documented and the documenta- tion provided alongside the assets?
Answer: [Yes]
Justification: Our only asset is code to perform our experiments. This is documented in the code and README.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details of the dataset/code/mmoel part of their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset is used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and Research with Human Subjects
For crowdsourcing experiments and research with human subjects, does the pa- per include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Question:
Answer: [NA]
Justification: No human subjects were used.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contri- bution of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurlPS Code of Ethics, workers involved in data collection, cura- tion, or other labor should paid at least the minimum wage the country of the data collector.
15. Institutional Review Board (IRB) Approvals or Equivalent for Research with Human Subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or an equivalent approval/review based on the requirements of your country or institution) were obtained?
Answer: [NA]
Justification:
No IRB approval was required.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
30

---

Depending on the country in which research is conducted, IRB approval (or equiva- lent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
31


Reviews:
[{"review_id": "8djK5ZkBD8", "reviewer": "Reviewer_cS79", "summary": "The goal of this paper is to better understand the performance of using Lanczos for matrix-function times vector computations. Using Lanczos for this task is the de-facto standard, since it works very well (best vs competitors) in practice. Previous theory is very weak: there is a big gap between the bounds it specifies and the actual performance. The goal of this paper is shrink and where possible eliminate this gap.\n\nA stronger bound for Lanczos is established.  Extensive experiments are reported.", "questions": "- How hard is it to build rational approximation with poles outside the interval of eigenvalues? Can you give a reference to a general result?\n- Did you conduct experiments with finite precision to see how much the theory carries over?", "limitations": "Nothing to add.", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "Very well written paper.\n\nThe main value of the paper is that it takes a well established method and de-facto standard (using Lanczos for matrix-function times vector computations), but with weak theory, and establishes much stronger theory. The reason this is important is that the weak bound motivated researchers to find better algorithms. Such algorithms had better theoretical bounds, but in practice were not as good as Lanczos. This paper closes the gap, giving stronger theory for Lanczos.\n\nThe experimental evaluation is robust. Appendix also contains additional bounds specific for matrix square root and inverse square root.", "weaknesses": "- Bounds are theoretic, and only explain what an existing algorithm achieves. There is no new algorithm in the paper.\n- The assumption that the poles are outside the gap that contains the eigenvalues seem to imply that the bound works only for a very restricted subset of rational functions. How common is it?\n- Bounds are only for exact arithmetic. The picture for finite-precision arithmetic can be very different. Even though I can understand why bound for finite-precision arithmetic are outside the scope for this paper, the authors could have conducted experiments with them to see how much the bounds carry to finite-precision.\n- Even though the problem of computing f(A)b has many ML applications, at the core the paper is focused on the quality of a NLA problem (solving f(A)b). It is more a NLA paper than a ML paper.\n- In some experiments, in order to show significant gap between the new bound and the old bound (Fact 1), a large number of iteration are used so the error becomes very very small. For example, in Figure 2 (center) the error starts at 10^-5 and goes to 10^-69. In Figure 2 (right) errors start at 10^-10!  In applications you will rarely work in these regimes, and will do with much smaller error. It is unclear whether the gap between Fact 1 and the new bound is significant for that regime.\n- Lemma 1 seems a bit weak: the multiplicative factor or norm(b) is C_r which can be huge. Sure, it is present in (7), but there the minimum it multiplies presumed to be very small. In Lemma 1 it appears on norm(b), which is constant.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thorough and thoughtful comments! This review raises many valuable points:\n\nRegarding the difference between exact arithmetic and finite precision, please see the global response.\n\nThe review asks if it is realistic to build a rational approximation with poles outside the interval of eigenvalues. For many functions (e.g. exponential, $x^{\\alpha}$) explicit rational approximations with poles outside of the interval of eigenvalues are known. (See for instance, reference 36 and 62.) For other “nice” (e.g. continuous) functions, one could apply the Remez algorithm to find the best rational approximation with real poles (see reference 63, chapter 24-5). The resulting approximation will not have any pole inside the interval of approximation, because otherwise the approximation would have infinite error near that pole. For less nice functions (e.g. absolute value which has a discontinuous derivative), the best rational approximations actually have *complex* poles which are clustered about discontinuities in the (higher order) derivatives. For instance, for the absolute value, the best rational approximation has poles on the imaginary axis. Our theory does not apply to this case, but this would be one of the most natural classes of functions to explore next.\n\nFor experiments, we chose to run Lanczos in extended precision because it allows us to highlight the qualitative behavior of our bound: Theorem 1 captures the correct shape of the convergence curve better than Fact 1. Nevertheless, several of our plots (e.g., the left and center panels of Figure 2) show that our analysis can still significantly outperform Fact 1, even when the desired accuracy or number of iterations is low. For some problems, this is not the case, mainly because of the large prefactor in Theorem 1, which depends on the condition number of $A_1, \\ldots, A_q$. Such a prefactor does not appear in Fact 1, so our new bound is only better for a large number of iterations. That said, as discussed in Section 4, we believe the prefactor can likely be reduced with a better theoretical analysis, in which case the improvement of Theorem 1 would kick in sooner.\n\nThe review also asks about the multiplicative factor in Lemma 1. It is correct to note that in Equation 8, $\\||b\\||_2$ is a constant. However, in this term, $C_r$ is *also* multiplied by $\\||f - r\\||_I$, which denotes the approximation error in the max-norm for the rational function approximation $r$. As the degree of $r$ grows, this goes to zero, often much more quickly than the error of the best polynomial approximation to $f$ on the interval $I$.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bfUJa9kLD3", "reviewer": "Reviewer_CAYW", "summary": "The authors study computation of matvecs of a matrix function via the Lanczos method, an in particular try to answer the question of why the basic Lanczos method is competitive or even superior to sophisticated implementations targeting specific matrix functions. They provide theoretical bounds on the error of Lanczos methods for rational functions which they use to produce a bound on generic functions based on distance to the nearest rational, and demonstrate this in numerical experiments.", "questions": "Congratulations on this paper.", "limitations": "The theoretical results rely on the poles of the rational function lying outside of the spectrum of the operator which might not always be the case, but for many practical situations it will be, and the authors' assertion that this case be considered future work is convincing to me.", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The first major contribution of this article lies in the introduction section, which is establishment of the fact that Lanczos methods outperform sophisticated, problem specific peers in some instances. This \"theory-practice gap\" in and of itself is an important and consequential fact to establish. The numerical experiments they do to demonstrate this are convincing and well communicated. \n\nSubsequently, they offer important theorems governing the ability of Lanczos methods to approximate rational functions, and use this to find a bound based on the l infinity distance between a given function and the nearest rational function of low degree. Numerical experiments demonstrate these results on important practical cases. \n\nThe way that the main theorem was sketched as special case in Section 2.1 is an exceptional move from a mathematical communication perspective. I was easily able to follow that specific proof, and it gave me intuition as to why the general case might be true without going through the more extensive proof in the appendix. I wish more articles did this, and I will try to do this in my future articles.", "weaknesses": "In practice, we care about the effectiveness of this method in floating point arithmetic. But I think it makes sense to partition that off as a separate project, as the contents of the article as it stands are convincing on their own. \n\nHere are some grammar issues I noticed (did not affect my scoring):\n33) functions repeated.\n72) an problem instance -> a problem instance\n146) the \"and <=\" may be missing a term.\n147 and 148) May be worth clarifying that it is the right hand side of (6) that is being referred to.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thorough review! The point about floating point arithmetic is well-taken; see the global response for more. Thank you as well for the grammar corrections. These have all been corrected in our revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bZrrzdt0oz", "reviewer": "Reviewer_nxLf", "summary": "The submission discusses the optimality (in terms of approximation error) of approximation $(A, b) \\rightarrow f(A)b$ with Lanczos' algorithm.\nThe main result (Theorem 1) states that the Lanczos iteration is \"near instance optimal\". Near-instance-optimality relates to the optimal reconstruction of $r(A)b$ for rational functions $r$. \nVia the triangle inequality, one can deduce similar statements for general $f$ provided $f$ can be approximated well by rational functions.\nThe submission proves Theorem 1 by repeatedly applying a similar optimality result for $A^{-1} b$ (i.e. $f=1/x$), which is why the final bound comes with scalar $\\kappa(A)^q$ where $\\kappa$ is the condition number of $A$ and $q$ the order of the denominator of the rational function.", "questions": "What follows is not directly a question, but the submission's format appears to deviate from the Neurips style in two ways:\n\n1. The abstract contains two paragraphs, even though the template states, \"The abstract must be limited to one paragraph.\"\n2. It seems that the mathematics font has been changed from the original Neurips template. \n\nI am unsure what Neurips' policy is here. Should these two things be changed?", "limitations": "Limitations are acknowledged. (I discuss this under \"Weaknesses\" above.)", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The paper provides an interesting perspective on the performance of Lanczos' algorithm for computing matrix functions.\nThe main strength is clarity, partly because the result is relatively strong and partly because the presentation is easy to follow.\nI appreciate the proof sketch in Section 2.1 and the comparison to prior work in Appendix B (especially in B.2, which I was wondering about while reading the other parts of the manuscript).\nI also appreciate the code submission.\nOverall, I think this is a nice paper.", "weaknesses": "The paper's main contribution is Theorem 1: near-instance optimality of Lanczos' method for computing $f(A)b$.\nThe weaknesses of this result are twofold:\n1. The analysis is limited to exact arithmetic. It is well-known that implementations of Lanczos' method in exact versus finite-precision arithmetic may not match. Focussing on exact arithmetic somewhat limits the practical applicability of Theorem 1.  \n2. The bound in Theorem 1 includes the constant $\\kappa(A)^q$. Similar works (including those discussed in Appendix B) achieve smaller constants, and in practical applications, $\\kappa(A)$ can be gigantic. Since the largest $\\kappa$ in the experiments seems to be $10^6$ (unless I have missed something), it is slightly unclear whether the bounds actually do explain the superiority of Lanczos' algorithm in practice (where $\\kappa$ can exceed $10^6$ by a large margin; see, for instance,  the condition number of Gaussian process covariance matrices or discretisations of PDEs).\n\nThat said, both weaknesses are openly discussed as limitations.\n\nFurthermore, Krylov methods have gained popularity in machine learning in recent years, as evidenced by the growing number of papers/submissions on Lanczos, Arnoldi, CG, GMRES, etc. \nNevertheless, the submission's connection to existing machine learning literature would be strengthened if at least some of the matrix functions mentioned in the introduction would reappear in the experiments.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thorough and thoughtful comments. This review raises several valuable points:\n\nOn the issues of exact vs. floating point arithmetic, please see the global response.\n\nRegarding the leading constant in Theorem 1, it is true that our result is weaker if the condition number of each $A_i$ is large. This means that, as Figure 2 shows, it is sometimes better than existing bounds (Fact 1) and sometimes worse. See the discussion below Theorem 1, at the bottom of page 4. However, as discussed in Section 4, we believe that the dependence on the condition number in the leading constant $C$ can likely be improved significantly with a tighter theoretical analysis, bringing the bound closer to more specialized results (for specific matrix functions), like those discussed in Appendix B.1. \n\nThe review notes that it would be valuable for the paper to touch on the matrix functions mentioned in the introduction. While our main focus is rational functions, we do study several of the functions most relevant for machine learning both theoretically and empirically:\n- Matrix square root and inverse square root: Analyzed theoretically in Appendix D. Experiments in Appendix E.1. Also discussed below Lemma 1.\n- Matrix exponential: For theoretical analysis, see the discussion below Lemma 1. Experiments are shown in the right panel of Figure 1. A rational approximation to it is studied in Figure 2. See the discussion accompanying these figures at the beginning of Section 3.\n- Matrix log: Experiments are shown in the right panel of Figure 1. A rational approximation to it is studied in Figure 2. See the discussion accompanying these figures at the beginning of Section 3.\n- Matrix sign: This is discussed at the end of Section 3.2. Results are shown in Figure 5 and Figure 6\n\nThank you for spotting the deviations from the NeurIPS style. These were oversights and we are more than happy to correct them for the revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "jCKyCO8ZCD", "reviewer": "Reviewer_uopa", "summary": "A matrix function f(A) for a real, symmetric matrix A and a univariate function f is given as \\sum f(\\lambda_i) u_i u_i^\\top, where \\lambda_i are the eigenvalues of A and u_i the corresponding eigenvectors. In practice, one is often interested in matrix-vector products f(A)b for some vector b. In the paper, the approximation quality of the Lanczos-FA method for approximating f(A)b for rational functions f is analyzed theoretically and experimentally. In particular, a known bound that is uniform in the A and b is improved to instance specific bounds using the notion of near instance optimality. The instance specific bounds are not tight. However, they can better explain better the experimentally observed good approximation quality of the Lanczos-FA method that often outperforms more advanced Krylov subspace methods.", "questions": "1. I like the proof sketch in Section 2.1, which was easy to follow for a non-expert like me. The only thing that I did not see immediately was the equality below Line 145. It is probably an elementary result, but maybe you can provide a few more insights why this equality holds true.\n\nIs well addressed in the author's rebuttal.\n\n2. In the second paragraph of Section 3 (Experiments): Why did you choose the matrices (or the spectra of the matrices) as you did? How do practical spectra look like? What are extreme/challenging spectra?\n\nIs also well addressed in the author's rebuttal.", "limitations": "As far as I can tell, there are no limitations. The bounds are not tight, but I do not consider this a limitation of the paper.", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "Originality: The paper provides new theoretical insights into the observed practical behavior of the Lanczos-FA method. It contributes to a better understanding of Krylov subspace methods. The results are non-trivial.\n\nClarity:  The paper is well is structured and, in general, easy to read. I did not read the full proof of Theorem 1, but appreciate the proof sketch in Section 2.1.  \n\nSignificance: As far as I can tell, the paper addresses an interesting, practically relevant problem. The paper does not only present upper bounds but also aims to find instances for the Lanczos-FA.", "weaknesses": "As far as I can tell, there are no major weaknesses. \n\nMinor typos and suggestions for improvement:\n\nIntroduction: Bayseian\n\nSection 1.4: ... functions of interest(ing) ...\n\nFootnote 2: Should \\|x\\|_A be \\|b\\|_a?\n\nLine 148: The (second factor) of the second term ...\n\nFigure 2: The axis labels are hardly legible.\n\nFigure 5: matvec -> matrix-vector product", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for the thoughtful feedback! We will address the two questions in turn:\n\nThe first asks for clarification about Line 145. The projection of a vector $x$ onto a linear subspace $S$ is defined as $\\underset{y \\in S}{\\operatorname{argmin}} \\||y - x\\||_2$. The $A$-norm projection is defined as $\\underset{y \\in S}{\\operatorname{argmin}} \\||y - x\\||_A = \\underset{y \\in S}{\\operatorname{argmin}} \\||A^{1/2} y - A^{1/2} x\\||_2$. In our case, the linear subspace is the Krylov subspace $K_k(A,b)$. By definition, $Q$ is a basis for this subspace. Therefore, any $y \\in K_k(A,b)$ can be written as $Qc$ for some $c$. Thus, the $A$-norm projection can be written $\\underset{c}{\\operatorname{argmin}} \\||A^{1/2} Q c - A^{1/2} x\\||_2$. Solving this least squares problem using the normal equations, we find that $c = (Q^\\top A Q)^{-1} Q^\\top A x = T^{-1} Q^\\top A x$. Therefore, the projection is $y = Qc = Q T^{-1} Q^\\top A x$. In other words, $\\||x - Q T^{-1} Q^\\top A x \\||_2 = \\underset{y \\in K_k}{\\operatorname{min}} \\||x - y\\||_2$. Finally, by definition, any vector in the Krylov subspace can be written as $p(A)b$ for some polynomial $p$ with degree $< k$, so we have $\\underset{y \\in K_k(A, b)}{\\operatorname{min}} \\||x - y\\||_2 = \\underset{p}{\\operatorname{min}} \\||x - p(A)x\\||_2$. To finish the argument, replace $x$ with $A^{-2}b$. We will modify the proof sketch to clarify this step.\n\nThe second asks about the choice of matrices for our experiments. First, note that all Krylov subspace methods are equivariant to orthogonal transformations, and the errors $\\||f(A)b - \\mathrm{alg}\\||$ are invariant to such transforms (in the Euclidean norm). So, without loss of generality, we choose $A$ to be a diagonal matrix, as is standard in the literature. \n\nIt is well understood in the literature that the eigenvalue distribution has a large influence on the performance of Krylov subspace methods, often in subtle ways. For instance, see reference [8], “Towards understanding CG and GMRES through examples”. We used uniform, skewed, and clustered eigenvalue distributions to get some variety in the convergence behavior (as reflected in the different shapes of the curves in Figure 2). In particular, eigenvalue distributions that are clustered at the ends of the spectrum often have interesting behavior. A non-uniform distribution can emphasize the difference between the uniform optimality of Fact 1 and the instance optimality of our bounds. As we report in Section 3.1, distributions with a single outlying eigenvalue are among the most challenging for Lanczos-FA.\n\nWhat practical spectra look like depends a lot on the specific problem, but they are often characterized by the presence of many small eigenvalues and a smaller number of large outlying and clustered eigenvalues (hence, why we consider skewed and clustered eigenvalue distributions). \n\nWe also thank you for spotting these typos. They will all be corrected in the revision.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "IdDcDI4Hr5", "reviewer": "Reviewer_dxxb", "summary": "The authors investigate the convergence of the Lanczos method for computing $f(A)v$ for a matrix $A$ and a vector $v$. They show a near instance optimality bound that involves the product of the condition numbers of matrices related to $A$ and empirically investigate the bound.", "questions": "- Can you generalize your results to other types of Krylov subspace methods such as rational Krylov methods?\n- In the middle and the right figures in Figure 1, the relative error seems very small. Did you use higher precision floating point format than double precision?", "limitations": "The authors discuss limitations in Section 4.", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "The authors investigate the convergence of the Lanczos method for computing matrix functions theoretically and numerically. The results are also meaningful for the machine learning community since matrix functions appear in many machine learning algorithms.", "weaknesses": "I think the authors should explain more about how we can interprete the proposed bound. The authors insist the bound captures the convergence behavior of the Lanczos method. The proposed bound involves the condition numbers of the matrices appearing in the denominator of the rational function. I think it is natural, but I couldn't clearly understand why the proposed bound can explain the convergence behavior better than the existing bound in Fact 1,  since if the condition numbers of these matrices are large, than the magnitude of $f$ in the bound in Fact 1 can be large, and the factor involving $f$ can also become large.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful feedback. We will respond to each of the questions in turn:\n\n> “Can you generalize your results to other types of Krylov subspace methods such as rational Krylov methods?”\n\nIt would be interesting to see whether bounds of the flavor we describe can be extended to rational Krylov subspace methods. However, the class of functions we are studying are already rational functions. If one were to use a rational Krylov subspace method (with a suitably chosen set of poles), these functions could be applied exactly. So it is somewhat less clear exactly what the generalization would look like.\n\n> “In the middle and the right figures in Figure 1, the relative error seems very small. Did you use higher precision floating point format than double precision?”\n\nYes, we used higher precision arithmetic. The precise setup is described in the first paragraph of Section 3.\n\n> “I couldn't clearly understand why the proposed bound can explain the convergence behavior better than the existing bound in Fact 1”\n\nIt is true that our main bound (Theorem 1) is weaker if the condition number of each $A_i$ is large. This means that, as Figure 2 shows, it is sometimes better than Fact 1 and sometimes worse. See the discussion below Theorem 1, at the bottom of page 4. However, importantly, while it depends on the condition numbers, our leading constant factor $C$ does *not* depend on the number of Lanczos iterations. On the other hand, the term $\\min_{\\mathrm{deg}(p) < k-q+1} \\||r(A)b - p(A)b\\||_2$ in Theorem 1 does: it decreases with the number of iterations, $k$, typically at a much faster rate than Fact 1. For this reason, Theorem 1 usually beats Fact 1 eventually as the number of iterations increases. If our current $C$ factor can be improved (we suspect it can be), then Theorem 1 would beat Fact 1 earlier.\n\nAs an extreme example, it might help to consider the simplified case of a 2x2 matrix with only two eigenvalues: $\\lambda_{\\min}$ and $\\lambda_{\\max}$. Let’s compare Definition 2, which is the form of our main theorem, to Fact 1. Both of them upper bound the error by a minimization problem over polynomials of a certain degree. To get a good guarantee from Fact 1, we would have to find a degree-$k$ polynomial that closely matches $f$ on the entire range $[\\lambda_{\\min}, \\lambda_{\\max}]$. When the condition number of $A$ is large, this range will be large, and so this will be impossible. However, to get a good guarantee from Definition 2, we only need to find a polynomial that matches $f$ at two discrete points: $\\lambda_{\\min}$, and $\\lambda_{\\max}$. This is much easier to do. In fact, if $k \\geq 1$, then our bound guarantees that we get zero error for any $\\lambda_{\\min}$, $\\lambda_{\\max}$, and any $f$. Fact 1 does not capture this behavior.", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:
We thank the reviewers for their thoughtful comments. We were pleased with the quality and number of reviews we received, and have taken the feedback into account in order to improve the paper. In our global response we address some common points raised by the reviewers. 

Several of the reviewers raised very reasonable questions about the applicability of our bounds to finite precision arithmetic. Generally speaking, the behavior of Lanczos in finite precision arithmetic is significantly more complicated than exact arithmetic. Therefore, bounds for exact arithmetic, such as the ones in this paper, serve as a starting point for a more complete understanding of the algorithm in finite precision arithmetic. We do observe that, typically, Lanczos-FA continues to perform near-optimally, even when run with the standard 16 digits of precision. Concretely, plots like those in Figure 1 and Figure 2 can be reproduced in finite precision. We will consider adding such plots to our revised version of the paper (as suggested by Reviewer cS79) to highlight this point.

We make three additional comments:

- First, we note that reorthogonalization is often used in practice. In this case, the behavior of Lanczos in finite precision is very similar to exact arithmetic and we expect our bounds would carry over with a little work.
- Second, even without reorthogonalization, there is some hope our theory can be applied. In particular, the analysis of Greenbaum (ref [32]) guarantees that the behavior of Lanczos in finite precision arithmetic is equivalent (in a precise way) to the behavior of Lanczos in exact arithmetic on a certain related problem. This allows exact arithmetic theory (e.g., our bounds) to be transferred to finite precision arithmetic. We discuss this in the final paragraph of Section 4. However, the precise equivalence is somewhat complicated. This would be an interesting direction for further work.
- Finally, we note that the analogous optimality guarantees for conjugate gradient in exact arithmetic are generally regarded as very important (even if they do not hold precisely in finite precision), because they provide intuition about how the algorithm works and are predictive of its excellent performance in practice.

Meta-review: The presented paper studies the approximation of the matrix function using a Lanczos-based method. It shows that, for a natural class of rational functions, Lanczos matches the error of the best possible Krylov subspace method, up to a multiplicative approximation factor.

Five reviewers evaluated the paper, and their overall assessment is positive. I agree with their evaluation and believe the paper makes a significant contribution with compelling results. In particular, I think the presentation of the paper is truly exceptional. Some reviewers raised concerns about the applicability of the proposed bounds in finite precision arithmetic, which the authors addressed satisfactorily.
Decision: Accept (spotlight)