PDF
\(\newcommand{\footnotename}{footnote}\) \(\def \LWRfootnote {1}\) \(\newcommand {\footnote }[2][\LWRfootnote ]{{}^{\mathrm {#1}}}\) \(\newcommand {\footnotemark }[1][\LWRfootnote ]{{}^{\mathrm {#1}}}\) \(\let \LWRorighspace \hspace \) \(\renewcommand {\hspace }{\ifstar \LWRorighspace \LWRorighspace }\) \(\newcommand {\mathnormal }[1]{{#1}}\) \(\newcommand \ensuremath [1]{#1}\) \(\newcommand {\LWRframebox }[2][]{\fbox {#2}} \newcommand {\framebox }[1][]{\LWRframebox } \) \(\newcommand {\setlength }[2]{}\) \(\newcommand {\addtolength }[2]{}\) \(\newcommand {\setcounter }[2]{}\) \(\newcommand {\addtocounter }[2]{}\) \(\newcommand {\arabic }[1]{}\) \(\newcommand {\number }[1]{}\) \(\newcommand {\noalign }[1]{\text {#1}\notag \\}\) \(\newcommand {\cline }[1]{}\) \(\newcommand {\directlua }[1]{\text {(directlua)}}\) \(\newcommand {\luatexdirectlua }[1]{\text {(directlua)}}\) \(\newcommand {\protect }{}\) \(\def \LWRabsorbnumber #1 {}\) \(\def \LWRabsorbquotenumber "#1 {}\) \(\newcommand {\LWRabsorboption }[1][]{}\) \(\newcommand {\LWRabsorbtwooptions }[1][]{\LWRabsorboption }\) \(\def \mathchar {\ifnextchar "\LWRabsorbquotenumber \LWRabsorbnumber }\) \(\def \mathcode #1={\mathchar }\) \(\let \delcode \mathcode \) \(\let \delimiter \mathchar \) \(\def \oe {\unicode {x0153}}\) \(\def \OE {\unicode {x0152}}\) \(\def \ae {\unicode {x00E6}}\) \(\def \AE {\unicode {x00C6}}\) \(\def \aa {\unicode {x00E5}}\) \(\def \AA {\unicode {x00C5}}\) \(\def \o {\unicode {x00F8}}\) \(\def \O {\unicode {x00D8}}\) \(\def \l {\unicode {x0142}}\) \(\def \L {\unicode {x0141}}\) \(\def \ss {\unicode {x00DF}}\) \(\def \SS {\unicode {x1E9E}}\) \(\def \dag {\unicode {x2020}}\) \(\def \ddag {\unicode {x2021}}\) \(\def \P {\unicode {x00B6}}\) \(\def \copyright {\unicode {x00A9}}\) \(\def \pounds {\unicode {x00A3}}\) \(\let \LWRref \ref \) \(\renewcommand {\ref }{\ifstar \LWRref \LWRref }\) \( \newcommand {\multicolumn }[3]{#3}\) \(\require {textcomp}\) \(\newcommand {\intertext }[1]{\text {#1}\notag \\}\) \(\let \Hat \hat \) \(\let \Check \check \) \(\let \Tilde \tilde \) \(\let \Acute \acute \) \(\let \Grave \grave \) \(\let \Dot \dot \) \(\let \Ddot \ddot \) \(\let \Breve \breve \) \(\let \Bar \bar \) \(\let \Vec \vec \)

Block cross-interactive residual smoothing for Lanczos-type solvers
for linear systems with multiple right-hand sides

Kensuke Aihara, Akira Imakura, Keiichi Morikuni

Abstract

Block Lanczos-type solvers, such as the block BiCGSTAB method [3], for large sparse linear systems

\begin{equation*} AX = B, \quad A \in \mathbb {R}^{n\times n}, \quad B \in \mathbb {R}^{n\times s}, \quad s \ll n \end{equation*}

often exhibit large oscillations in the residual norms. In finite precision arithmetic, the large oscillations lead to a large residual gap (the difference \(G_{R_k}\) between the recursively updated residual \(R_k\) and the explicitly computed residual \(B-AX_k\)) and a loss of attainable accuracy of the approximations, as observed in

\begin{equation*} \| G_{R_k} \| - \| R_k \| \leq \| B - A X_k \| \leq \| G_{R_k} \| + \| R_k \|, \quad G_{R_k} := (B-AX_k)-R_k. \end{equation*}

This problem is addressed by using cross-interactive residual smoothing (CIRS). Just as the standard Lanczos-type solvers for a single linear system have been extended to their global and block versions for solving systems with multiple right-hand sides, similar extensions of CIRS are naturally considered. While we have developed the global CIRS scheme (Gl-CIRS) in our previous study [1], we propose a block version (Bl-CIRS) in this study. Then, we demonstrate the effectiveness of Bl-CIRS from various perspectives, such as theoretical insights into the convergence behaviors of the residual and approximation norms, numerical experiments on model problems, and a detailed rounding error analysis for the residual gap. In particular, we show for the case of Bl-CIRS that orthonormalizing the columns of direction matrices plays an important role in reducing the residual gap. The presented analysis also complements our previous study above that includes an evaluation for the residual gap of the block Lanczos-type solvers.

Advances in residual smoothing Block Lanczos-type solvers typically update the \(k\)th approximation \(X_k\) and residual \(R_k\) by using the recursion formulas

\begin{equation*} X_{k+1} = X_k + P_k \alpha _k^\square , \quad R_{k+1} = R_k - (AP_k) \alpha _k^\square , \quad k = 0,1,\dots , \end{equation*}

respectively, where \(P_k \in \mathbb {R}^{n\times s}\) is a direction matrix and \(\alpha _k^\square \in \mathbb {R}^{s\times s}\) is determined under a certain condition. Residual smoothing was introduced by Schönauer [7] to Lanczos-type solvers for a single linear system to get a non-increasing sequence of residual norms [8]. A block version of the simple residual smoothing (Bl-SRS) was presented by Jbilou [5]. Let \(X_k\) and \(R_k\) be the \(k\)th primary approximation and residual, respectively. Then, new sequences of approximations \(Y_k\) and the corresponding smoothed residuals \(S_k (:= B - A Y_k)\) are generated by the recursion formulas

\begin{equation*} Y_{k+1} = Y_k + (X_{k+1} - Y_k) \eta _{k+1}^\square , \quad S_{k+1} = S_k + (R_{k+1} - S_k) \eta _{k+1}^\square , \quad k = 0,1,\dots , \end{equation*}

respectively, where \(Y_0 = X_0\), \(S_0 = R_0\), and \(\eta _k^\square \in \mathbb {R}^{s\times s}\) is a parameter matrix. With a local minimization of the smoothed residual norm in choosing \(\eta _k^\square \), a monotonically decreasing sequence of \(\| S_k \|\) is obtained.

Studies on the relationship between residual smoothing and the residual gap have an interesting history. For a single right-hand side case, Gutknecht and Rozložník [4] clarified that the conventional smoothing schemes (including the Zhou–Walker implementation [10]) do not help to improve the attainable accuracy. To be more specific, rounding errors accumulated in the primary sequences propagate to the smoothed sequences, and the smoothed true residual norms stagnate at the same order of magnitude as the primary ones. In order to remedy this phenomenon, Komeyama et al. [2, 6] modified the Zhou–Walker implementation so that the primary and smoothed sequences influence one another. This modification is referred to as cross-interactive residual smoothing (CIRS) and is indeed effective in reducing the residual gap and increasing attainable accuracy. As SRS has been extended to global and block versions [9, 5], CIRS is also extended. In this perspective, our previous study [1] presented a global version of CIRS (Gl-CIRS) for the global- and block-type solvers, and therefore, we propose a block version of CIRS (Bl-CIRS) in this study. Table 1 summarizes the recursion formulas for the aforementioned residual smoothing schemes. This table shows that this study fills a gap in the literature of the CIRS schemes.

Table 1: Difference in the recursion formulas for updating smoothed residuals.

.
Type SRS scheme CIRS scheme
.
Standard
\(A\boldsymbol {x}=\boldsymbol {b}\)
.
[7, 8]
\(\boldsymbol {s}_k = \boldsymbol {s}_{k-1} + \eta _k (\boldsymbol {r}_k-\boldsymbol {s}_{k-1})\),
\(\boldsymbol {r}_k, \boldsymbol {s}_{k-1} \in \mathbb {R}^n\),  \(\eta _k \in \mathbb {R}\)
.
[2, 6]
\(\boldsymbol {s}_k = \boldsymbol {s}_{k-1} - \eta _k A \boldsymbol {v}_k\),
\(\boldsymbol {v}_k, \boldsymbol {s}_{k-1} \in \mathbb {R}^n\),  \(\eta _k \in \mathbb {R}\)
.
Global
\(AX=B\)
.
[9]
\(S_k = S_{k-1} + \eta _k (R_k-S_{k-1})\),
\(R_k, S_{k-1} \in \mathbb {R}^{n\times s}\),  \(\eta _k \in \mathbb {R}\)
.
[1]
\(S_k = S_{k-1} - \eta _k AV_k\),
\(V_k, S_{k-1} \in \mathbb {R}^{n\times s}\),  \(\eta _k \in \mathbb {R}\)
.
Block
\(AX=B\)
.
[5]
\(S_k = S_{k-1} + (R_k-S_{k-1}) \eta _k^\square \),
\(R_k, S_{k-1} \in \mathbb {R}^{n\times s}\),  \(\eta _k^\square \in \mathbb {R}^{s\times s}\)
.
Present study
\(S_k = S_{k-1} - (A\tilde {Q}_k)\tilde {\eta }_k^\square \),
\(\tilde {Q}_k,S_{k-1}\in \mathbb {R}^{n\times s}\),  \(\tilde {\eta }_k^\square \in \mathbb {R}^{s\times s}\)

Block cross-interactive residual smoothing This study proposes updating approximations and the corresponding residuals by the recursion formulas

\begin{align*} & \text {smoothed} & Y_{k+1} & = Y_k + V_{k+1} \eta _{k+1}^\square , & S_{k+1} & = S_k - (AV_{k+1}) \eta _{k+1}^\square , \\ & \text {primary} & X_{k+1} & = Y_{k+1} + V_{k+1} (I_s - \eta _{k+1}^\square ), & R_{k+1} & = S_{k+1} - (AV_{k+1}) (I_s -\eta _{k+1}^\square ), \end{align*} respectively, for \(k = 0,1,\dots \) with \(Y_0 = X_0\) and \(S_0 = R_0\) so that the primary and smoothed sequences influence one another, where \(V_{k+1} = V_k (I_s - \eta _k^\square ) + \tilde {P}_k\) is an auxiliary matrix for \(\eta _0^\square = O \in \mathbb {R}^{s\times s}\) and \(V_0 = O \in \mathbb {R}^{n\times s}\). Here, \(\tilde {P}_k \in \mathbb {R}^{n\times s}\) is a direction matrix in the recursion formula \(X_{k+1} = X_k + \tilde {P}_k\). Again with a local minimization of the smoothed residual norm in choosing \(\eta _k^\square \), a monotonically decreasing sequence of \(\| S_k \|\) is obtained. Note that the essential difference of Bl-CIRS from Gl-CIRS [1, Algorithm 3.1] is that the smoothing parameter \(\eta _k^\square \) of Bl-CIRS is an \(s\)-by-\(s\) matrix instead of a scalar.

For numerical stability, Bl-CIRS needs to orthonormalize the columns of the auxiliary matrix \(V_k\) for each iteration. Let \(V_k = \tilde {Q}_k \tilde {\xi }_k\) be the QR decomposition of \(V_k\), where \(\tilde {Q}_k \in \mathbb {R}^{n\times s}\) and \(\tilde {\xi }_k \in \mathbb {R}^{s\times s}\) are the Q- and R-factors, respectively. With the auxiliary matrix \(\tilde {\eta }_k^\square := \tilde {\xi }_k \eta _k^\square \), the above formulas are equivalently rewritten as

\begin{align*} & \text {smoothed} & Y_{k+1} & = Y_k + \tilde {Q}_{k+1} \tilde {\eta }_{k+1}^\square , & S_{k+1} & = S_k - (A\tilde {Q}_{k+1}) \tilde {\eta }_{k+1}^\square , & \\ & \text {primary} & X_{k+1} & = Y_{k+1} + \tilde {Q}_{k+1} (\tilde {\xi }_{k+1} - \tilde {\eta }_{k+1}^\square ), & R_{k+1} & = S_{k+1} - (A \tilde {Q}_{k+1}) (\tilde {\xi }_{k+1} - \tilde {\eta }_{k+1}^\square ) & \end{align*} together with \(V_{k+1} = \tilde {Q}_k (\tilde {\xi }_k - \tilde {\eta }_k^\square ) + \tilde {P}_k\) for \(k = 0,1,\dots \).

Our main results via a rounding error analysis shows that Bl-CIRS with the orthonormalization strategy suppresses the residual gap.

  • Theorem 1. In finite precision arithmetic, let \(X_k \in \mathbb {F}^{n\times s}\) and \(R_k \in \mathbb {F}^{n\times s}\) be the \(k\)th approximation and residual generated by the recursion formulas

    \begin{equation} X_k = X_{k-1} + \widehat {Q}_{k-1} \alpha _{k-1}^\square , \quad R_k = R_{k-1} - (A\widehat {Q}_{k-1}) \alpha _{k-1}^\square , \quad k = 0,1,\dots , \label {eq:X_k,R_k} \end{equation}

    respectively, where \(\mathbb {F} \subset \mathbb {R}\) is a set of floating point numbers and \(\widehat {Q}_{k-1}\) is a Q-factor of the direction matrix \(P_{k-1}\) obtained from the QR decomposition with Givens rotations or Householder transformations. Then, with \(X_0 = O\), the residual gap \(G_{R_k}\) satisfies the bound

    \begin{equation*} \| G_{R_k} \| < \left ( 8 \sqrt {s} \gamma _{m+3s} + \gamma _1 \right ) k \| A \| \max _{0<i\leq k} \| X_i \| + k\gamma _1 \max _{0<i\leq k} \| R_i \|, \end{equation*}

    where \(\gamma _k := k \mathbf {u} / (1-k\mathbf {u})\) with a unit roundoff \(\mathbf {u}\) and \(m\) is the maximum number of nonzero entries per row of \(A\).

In the case of Bl-CIRS, replacing \(X_k\) and \(R_k\) by \(Y_k\) and \(S_k\), respectively, in (1), this theorem holds for the residual gap \(G_{S_k} = (B-AY_k)-S_k\). Therefore, even when using an inexact orthonormalization for the columns of iteration matrices, Bl-CIRS in which the residual and approximation norms converge smoothly is indeed useful to reduce the residual gap. This theoretical result is consistent with our numerical results. Numerical experiments demonstrate that Bl-CIRS is effective for suppressing the residual gap and improving the attainable accuracy of the approximations.

References