Nearly Optimal Approximation of Matrix Functions by the Lanczos Method

Thank you for your thoughtful feedback. We will respond to each of the questions in turn:

“Can you generalize your results to other types of Krylov subspace methods such as rational Krylov methods?”

It 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.

“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?”

Yes, we used higher precision arithmetic. The precise setup is described in the first paragraph of Section 3.

“I couldn't clearly understand why the proposed bound can explain the convergence behavior better than the existing bound in Fact 1”

It 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.

As 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.

Thank you for the thoughtful feedback! We will address the two questions in turn:

The 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.

The 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.

It 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.

What 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).

We also thank you for spotting these typos. They will all be corrected in the revision.

Thank you for the thorough and thoughtful comments. This review raises several valuable points:

On the issues of exact vs. floating point arithmetic, please see the global response.

Regarding 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.

The 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:

• Matrix square root and inverse square root: Analyzed theoretically in Appendix D. Experiments in Appendix E.1. Also discussed below Lemma 1.
• 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.
• 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.
• Matrix sign: This is discussed at the end of Section 3.2. Results are shown in Figure 5 and Figure 6

Thank you for spotting the deviations from the NeurIPS style. These were oversights and we are more than happy to correct them for the revision.

Thank you for the thorough and thoughtful comments! This review raises many valuable points:

Regarding the difference between exact arithmetic and finite precision, please see the global response.

The 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.

For 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.

The 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$.

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.