A posteriori error bounds for the block-Lanczos method for matrix function approximation

Abstract

We extend the error bounds from Chen et al. (SIAM J. Matrix Anal. Appl 43(2):787–811, 2022) for the Lanczos method for matrix function approximation to the block algorithm. Numerical experiments suggest that our bounds are fairly robust to changing block size and have the potential for use as a practical stopping criterion. Further experiments work toward a better understanding of how certain hyperparameters should be chosen in order to maximize the quality of the error bounds, even in the previously studied block size one case.

Appendices

Appendix 1. Estimates for the block-CG error norm

In this section, we discuss how to obtain an error estimate for the $(\textbf{H}-w\textbf{I})$-norm of the error of the block-Lanczos algorithm used to approximate $(\textbf{H}-w\textbf{I})^{-1}\textbf{V}$; i.e. with $f(x)=(x-w)^{-1}$. For simplicity, without loss of generality, we will assume $w=0$. Closely related approaches are widely studied for the block size one case [13, 28,29,30, 44]. Bounds for the block algorithm have also been studied [38].

Fig. 9

Illustration of the quality of Fig. 9 for varying d. The block-Lanczos-FA error is measured using $\Vert \cdot \Vert _{\textbf{H}-w\textbf{I}}$ with $w = 0$

By a simple triangle inequality, we obtain a bound

$$\begin{aligned} \Vert \textsf {err} _k(w) \Vert _{\textbf{H}}&= \Vert \textbf{H}^{-1} \textbf{V} - \textsf {lan} _k(f) \Vert _{\textbf{H}} \\ {}&\le \Vert \textbf{H}^{-1}\textbf{V} - \textsf {lan} _{k+d}(f) \Vert _{\textbf{H}} + \Vert \textsf {lan} _{k+d}(f) - \textsf {lan} _k(f)\Vert _{\textbf{H}} . \end{aligned}$$

If we assume that d is chosen sufficiently large so that

$$\begin{aligned} \Vert \textbf{H}^{-1}\textbf{V} - \textsf {lan} _{k+d}(f) \Vert _{\textbf{H}} \ll \Vert \textsf {lan} _{k+d}(f) - \textsf {lan} _k(f)\Vert _{\textbf{H}}, \end{aligned}$$

then we obtain an approximate error bound (error estimate)

$$\begin{aligned} \Vert \textsf {err} _k(w) \Vert _{\textbf{H}} \lesssim \Vert \textsf {lan} _{k+d}(f) - \textsf {lan} _k(f)\Vert _{\textbf{H}}. \end{aligned}$$

(18)

It is natural to ask why we don’t simply apply this approach to the block-Lanczos-FA algorithm for any function f. Such an approach was suggested for the Lanczos-FA algorithm in [12], and in many cases can work well. However, if the convergence of the block-Lanczos-FA algorithm is not monotonic, then the critical assumption that $\textsf {lan} _{k+d}(f)$ is a much better approximation to $f(\textbf{H})\textbf{V}$ than $\textsf {lan} _k(f)$ may not hold. Perhaps more importantly, the decomposition in Theorem 2.1 also allows intuition about how the distribution of eigenvalues of $\textbf{H}$ impact the convergence of block-CG to be extended to the block-Lanczos-FA algorithm.

In Fig. 9 we show a numerical experiment illustrating the quality of (18) for several values of d. As expected, increasing d improves the quality of the estimate. However, even for small d, Fig. 9 appears suitable for use as a practical stopping criterion on this problem. Of course, on problems where the convergence of block-CG stagnates, such an error estimate will require d fairly large. Further study of this topic is outside the scope of the present paper.

Appendix 2. Proofs of known results

Proof of (3)

Owing to linearity, in order to show $p(\textbf{H})\textbf{V} = \textsf {lan} _k(p)$ for any polynomial p with $\deg (p)<k$, it suffices to consider the case $p(x) = x^j$, for $j=0, \ldots , k-1$. Left multiplying (2) by $\textbf{H}^{j-1}$ and then repeatedly applying (2) we find

$$\begin{aligned} \textbf{H}^j \textbf{Q}_k&= \textbf{H}^{j-1} \textbf{Q}_k \textbf{T}_k + \textbf{H}^{j-1} \overline{\textbf{Q}}_{k+1} \textbf{B}_k \textbf{E}_k^*\\&= \textbf{H}^{j-2} (\textbf{Q}_k \textbf{T}_k + \overline{\textbf{Q}}_{k+1} \textbf{B}_k \textbf{E}_k^*) \textbf{T}_k + \textbf{H}^{j-1} \overline{\textbf{Q}}_{k+1} \textbf{B}_k \textbf{E}_k^*\\&\vdots \\&= \textbf{Q}_k \textbf{T}_k^{j} + \sum _{i=1}^{j} \textbf{H}^{j-i} \overline{\textbf{Q}}_{k+1} \textbf{B}_k \textbf{E}_k^* \textbf{T}_k^{i-1}. \end{aligned}$$

Thus, using that $\textbf{V} = \textbf{Q}_k \textbf{E}_1\textbf{B}_0$,

$$\begin{aligned} \textbf{H}^j \textbf{V} = \textbf{H}^j \textbf{Q}_k \textbf{E}_1\textbf{B}_0 = \textbf{Q}_k \textbf{T}_k^{j}\textbf{E}_1\textbf{B}_0 + \sum _{i=1}^{j} \textbf{H}^{j-i} \overline{\textbf{Q}}_{k+1} \textbf{B}_k \textbf{E}_k^* \textbf{T}_k^{i-1}\textbf{E}_1\textbf{B}_0. \end{aligned}$$

Since $\textbf{T}_k$ is banded with half bandwidths b, $\textbf{T}_k^{i-1}$ is banded with half bandwidth $(i-1)b$. Thus, the bottom left $b\times b$ block of $\textbf{T}_k^{i-1}$ is all zero provided $i-1 \le k-2$; i.e., $\textbf{E}_k^* \textbf{T}_k^{i-1}\textbf{E}_1\textbf{B}_0 = \textbf{0}$ provided $i<k$. We therefore find that for all $j<k$,

$$\begin{aligned} \textbf{H}^j \textbf{V} =\textbf{Q}_k \textbf{T}_k^{j}\textbf{E}_1\textbf{B}_0, \end{aligned}$$

from which we infer $p(\textbf{H})\textbf{V} = \textsf {lan} _k(p)$ for any p with $\deg (p)<k$. This is a well-known result; [10, 17, 18, 37].

Therefore, using the triangle inequality, for any polynomial p with $\deg (p)<k$,

$$\begin{aligned} \Vert f(\textbf{H}) \textbf{V} - \textsf {lan} _k(f) \Vert _2&= \Vert f(\textbf{H}) \textbf{V} - \textbf{Q}_kf(\textbf{T}_k) \textbf{E}_1 \textbf{B}_0 \Vert _2 \\ {}&\le \Vert f(\textbf{H}) \textbf{V} - p(\textbf{H}) \textbf{V} \Vert + \Vert \textbf{Q}_k p(\textbf{T}_k) \textbf{E}_1 \textbf{B}_0 - \textbf{Q}f(\textbf{T}_k) \textbf{E}_1 \textbf{B}_0 \Vert _2 \\ {}&\!\le \! \Vert f(\textbf{H}) \!-\! p(\textbf{H}) \Vert _2 \Vert \textbf{V} \Vert _2 \!+\! \Vert \textbf{Q}_k \Vert _2 \Vert p(\textbf{T}_k) \!-\! f(\textbf{T}_k) \Vert \Vert \textbf{E}_1 \textbf{B}_0 \Vert _2. \\&\le \Vert \textbf{V} \Vert _2 \left( \Vert f(\textbf{H}) - p(\textbf{H}) \Vert _2 + \Vert p(\textbf{T}_k) - f(\textbf{T}_k) \Vert _2 \right) . \end{aligned}$$

In the final inequality we have use that $\Vert \textbf{Q}_k\Vert \le 1$ since $\textbf{Q}_k$ has orthonormal columns and that $\Vert \textbf{E}_1 \textbf{B}_0\Vert = \Vert \textbf{Q}_k^*\textbf{V}\Vert \le \Vert \textbf{V}\Vert$. Next, note that

$$\begin{aligned} \Vert f(\textbf{H}) - p(\textbf{H}) \Vert _2 = \max _{x\in \Lambda (\textbf{H})} | f(x) - p(x) | ,\qquad \Vert p(\textbf{T}_k) - f(\textbf{T}_k) \Vert _2 = \max _{x\in \Lambda (\textbf{T}_k)} | f(x) - p(x) | \end{aligned}$$

so, since each of $\Lambda (\textbf{H})$ and $\Lambda (\textbf{T}_k)$ are contained in $[\lambda _{min }(\textbf{H}),\lambda _{max }(\textbf{H})]$, we find

$$\begin{aligned} \Vert f(\textbf{H}) \textbf{V} - \textsf {lan} _k(f) \Vert _2&\le 2 \Vert \textbf{V} \Vert \max _{x\in [\lambda _{min }(\textbf{H}),\lambda _{max }(\textbf{H})]} | f(x) - p(x) |. \end{aligned}$$

Since p was arbitrary, we can optimize over all polynomials p with $\deg (p)<k$.$\square$

Proof of Lemma 2.2

From Definition 1.2,

$$\begin{aligned} \textsf {res} _k(z) := \textbf{V} - (\textbf{H}-z\textbf{I}) \textbf{Q}_k (\textbf{T}_k - z\textbf{I})^{-1} \textbf{E}_1\textbf{B}_0. \end{aligned}$$

Using the fact that the Lanczos factorization (2) can be shifted,

$$\begin{aligned} (\textbf{H}-z\textbf{I}) \textbf{Q}_k = \textbf{Q}_k(\textbf{T}_k-z\textbf{I})+\overline{\textbf{Q}}_{k+1}\textbf{B}_k\textbf{E}_k^*. \end{aligned}$$

Thus, by substitution, we have

$$\begin{aligned} \textsf {res} _k(z)&= \textbf{V} - (\textbf{Q}_k(\textbf{T}_k-z\textbf{I})+\overline{\textbf{Q}}_{k+1}\textbf{B}_k\textbf{E}_k^*) (\textbf{T}_k - z\textbf{I})^{-1} \textbf{E}_1\textbf{B}_0\\&= \textbf{V} - \textbf{Q}_k\textbf{E}_1\textbf{B}_0-\overline{\textbf{Q}}_{k+1}\textbf{B}_k\textbf{E}_k^*(\textbf{T}_k - z\textbf{I})^{-1} \textbf{E}_1\textbf{B}_0 \end{aligned}$$

Since $\textbf{E}_1 = \left[ \textbf{I}_{b}, \textbf{0},..., \textbf{0 }\right] ^*$, where $\textbf{I}_{b}$ is the $b\times b$ identity matrix, we know that $\textbf{Q}\textbf{E}_1\textbf{B}_0 = \textbf{Q}_1\textbf{B}_0 = \textbf{V}$. Therefore,

$$\begin{aligned} \textsf {res} _k(z)&= \textbf{V} - \textbf{Q}_k\textbf{E}_1\textbf{B}_0-\overline{\textbf{Q}}_{k+1}\textbf{B}_k\textbf{E}_k^*(\textbf{T}_k - z\textbf{I})^{-1} \textbf{E}_1\textbf{B}_0\\&= -\overline{\textbf{Q}}_{k+1}\textbf{B}_k\textbf{E}_k^*(\textbf{T}_k - z\textbf{I})^{-1} \textbf{E}_1\textbf{B}_0. \end{aligned}$$

The result follows from the definition of $\textbf{C}_k(z)$. $\square$

Proof of Corollary 2.1

We have shown in Lemma 2.2 that $\textsf {res} _k(z) = \overline{\textbf{Q}}_{k+1}\textbf{B}_k \textbf{C}_k(z)$. Thus,

$$\begin{aligned} \textsf {res} _k(z) \!=\! \overline{\textbf{Q}}_{k+1}\textbf{B}_k \textbf{C}_k(z) \!=\! \overline{\textbf{Q}}_{k+1}\textbf{B}_k \textbf{C}_k(w)\textbf{C}_k(w)^{-1} \textbf{C}_k(z) \!=\! \textsf {res} _k(w) \textbf{C}_k(w)^{-1} \textbf{C}_k(z). \end{aligned}$$

Next, notice that

$$\begin{aligned} \textsf {err} _k(z) = (\textbf{H}-z\textbf{I})^{-1}\textsf {res} _k(z)&=(\textbf{H}-z\textbf{I})^{-1}\textsf {res} _k(w)\textbf{C}_k(w)^{-1} \textbf{C}_k(z)\\&= h_{w, z}(\textbf{H})(\textbf{H}-w\textbf{I})^{-1}\textsf {res} _k(w)\textbf{C}_k(w)^{-1} \textbf{C}_k(z)\\&=h_{w, z}(\textbf{H})\textsf {err} _k(w)\textbf{C}_k(w)^{-1} \textbf{C}_k(z). \end{aligned}$$

The result is established.$\square$