Code to replicate the figures in the corresponding paper can be found on Github.
Introduction
The Lanczos method for matrix function approximation (Lanczos-FA) can be used to approximate f(A)b (read more about this task here), and in the case that f(x)=1/x and A is positive definite, this approximation is optimal over Krylov subspace.
This case is very well studied and a range of error bounds and estimates exist.
However, for other functions, the standard bounds are often too pessimistic as they do not take into account fine grained information about the spectrum of A such as outlaying or clustered eigenvalues.
This makes it difficult to know when Lanczos-FA has reached a suitable accuracy for a given problem.
In this paper we show how to reduce the error of approximating f(A)b with Lanczos-FA to the error of solving a certain linear system with the Lanczos-FA.
This allows us to leverage the range of existing bounds for the convergence of Lanczos-FA on linear systems to easily obtain a priori and a posteriori bounds for other matrix functions including piecewise analytic functions such as the sign function.
Our a posteriori bounds are highly accurate and can be used as practical stopping criteria.
The basic idea
The k-th Lanczos-FA approximation to f(A)b is defined as
lank(f):=Qf(T)QTb,
where Q and T are produced by the Lanczos method run for k steps on (A,b).
Suppose A is a Hermitian matrix and If f is analytic in a neighborhood of the eigenvalues of A and Γ is a contour in this neighborhood containing these eigenvalues,
f(A)b=−2πi1∮Γf(z)(A−zI)−1bdz.
If the eigenvalues of T are also contained in Γ, we similarly have
Qf(T)QTb=−2πi1∮Γf(z)Q(T−zI)−1QTbdz.
Combining, these, we see that the Lanczos-FA error can be written as
f(A)b−Qf(T)QTb=−2πi1∮Γf(z)errk(z)dz.
For z∈C, define the k-th Lanczos-FA error and residual for the linear system (A−zI)x=b as,
errk(z,A,b)resk(z,A,b):=(A−zI)−1b−Q(T−zI)−1QTb,:=b−(A−zI)Q(T−zI)−1QTb.
It is a well-known fact that
resk(z)=(det(T−zI)(−1)kj=1∏kβj)∥b∥2qk+1.
Consequently, for all z,w∈C, where A−zI and A−wI are both invertible,
errk(z)resk(z)=det(hw,z(T))hw,z(A)errk(w)=det(hw,z(T))resk(w),
where
hw,z(x)=x−zx−w.
Thus, if Γ is a simple closed curve or union of simple closed curves inside this neighborhood and enclosing the eigenvalues of A and T and w a point not in Λ(T)∪Λ(A),
f(A)b−lank(f)=(−2πi1∮Γf(z)det(hw,z(T))hw,z(A)dz)errk(w).
Applying the triangle inequality and the submultiplicitivity of matrix-norms, we can then obtain a bound ∥f(A)b−lank(f)∥≤integral term(2π1∮Γ∣f(z)∣(i=1∏k∥hw,z∥Si)∥hw,z∥S0∣dz∣)linear system error∥errk(w)∥,
where Si are some suitably chosen sets and ∥g∥S:=maxx∈S∣g(x)∣.
Note that the integral term and linear system error term in the theorem are entirely decoupled!
Thus, once the integral term is computed, bounding the error of Lanczos-FA for f(A)b is reduced to bounding ∥errk(w)∥, and if the integral term can be bounded independently of k, implies that, up to a constant factor, the Lanczos-FA approximation to f(A)b converges at least as fast as ∥errk(w)∥.
Much of the paper is focused on practical aspects regarding the use of such a bound.
In particular:
we provide a detailed discussion and analysis of the validity of our bounds when the Lanczos iterate was computed using finite precision arithmetic.
we provide several analytic examples where the integral term is computed directly
we use numerical experiments to demonstrate the effectiveness of our bounds on a variety of functions such as the square root, inverse square root, and sign functions, and
we derive similar bounds for quadratic forms.
It’s worth noting that similar ideas have been previously used to get error bounds for Stieltjes functions [1–3].
As discussed in Section 2.2 of the paper, our bounds differ in a number of key ways.
One key difference is that we reduce error bounds for Lanczos-FA a general function to error bounds for Lanczos-FA on a fixed linear system.
This allows intuition about the convergence of Lanczos-FA on linear systems to be transferred to other functions.
In addition, our analysis allows bounds for piecewise continuous functions like the sign function.
References
1. Ilic, M.D.; Turner, I.W.; Simpson, D.P. A Restarted Lanczos Approximation to Functions of a Symmetric Matrix. IMA Journal of Numerical Analysis2009, 30, 1044–1061, doi:10.1093/imanum/drp003.
2. Frommer, A.; Güttel, S.; Schweitzer, M. Efficient and Stable Arnoldi Restarts for Matrix Functions Based on Quadrature. SIAM Journal on Matrix Analysis and Applications2014, 35, 661–683, doi:10.1137/13093491x.
3. Frommer, A.; Schweitzer, M. Error Bounds and Estimates for Krylov Subspace Approximations of Stieltjes Matrix Functions. BIT Numerical Mathematics2015, 56, 865–892, doi:10.1007/s10543-015-0596-3.