Admissible and attainable convergence behavior of block Arnoldi and GMRES

TL;DR

This paper extends the understanding of convergence behaviors for block GMRES and block Arnoldi methods, establishing what is achievable and how to construct matrices for desired convergence patterns, with implications for eigenvalue approximation.

Contribution

It develops convergence results for block GMRES and block Arnoldi by framing the problem over a ring of matrices, providing admissible convergence behaviors and construction methods.

Findings

Characterizes admissible convergence behaviors for block GMRES.

Provides methods to construct matrices and right-hand sides for specific convergence patterns.

Shows near independence of convergence between block Arnoldi and block GMRES for the same matrix.

Abstract

It is well-established that any non-increasing convergence curve is possible for GMRES and a family of pairs $(A,b)$ can be constructed for which GMRES exhibits a given convergence curve with $A$ having arbitrary spectrum. No analog of this result has been established for block GMRES, wherein multiple right-hand sides are considered. By reframing the problem as a single linear system over a ring of square matrices, we develop convergence results for block Arnoldi and block GMRES. In particular, we show what convergence behavior is admissible for block GMRES and how the matrices and right-hand sides producing any admissible behavior can be constructed. Moreover, we show that the convergence of the block Arnoldi method for eigenvalue approximation can be almost fully independent of the convergence of block GMRES for the same coefficient matrix and the same starting vectors.

Tables

Table 1. Table 1: Entities and operations on ℂ ℂ \mathbb{C} and their counterparts on 𝕊 𝕊 \mathbb{S} . CholU stands for the upper triangular factor of the Cholesky decomposition.

standard	block
$ℂ$	$𝕊\simeq {ℂ}^{s\times s}$
(commutative ∗-ring)	(noncommutative ∗-algebra over $ℂ$)
${ℝ}^{+}$	${𝕊}^{+}$…upper-$\Delta$ with positive diag. entries
(multiplicative topological group)	(multiplicative topological group)
${ℝ}_0^{+}$	${𝕊}_0^{+}$…upper-$\Delta$ with nonnegative diag. entries
$0$	singular $s\times s$ matrix
$1$	$I\in 𝕊$
$a,b\in ℂ$	$A,B\in 𝕊$
$\left|a\right|\equiv \sqrt{a^{\ast }a}\in {ℝ}_0^{+}$	$\left|A\right|=\sqrt{A^{\ast }A}\equiv \text{cholU}\left(A^{\ast }A\right)\in {𝕊}_0^{+}$
$\left|a\right|\in {ℝ}^{+}⟺a\ne 0$	$\left|A\right|\in {𝕊}^{+}⟺A\text{ nonsingular}$
$\left|a\left|b\right|\right|=\left|a\right|\left|b\right|$	$\left|A\left|B\right|\right|=\left|A\right|\left|B\right|$
$𝒙,𝒚\in {ℂ}^n$	$𝑿,𝒀\in {𝕊}^n\simeq {ℂ}^{ns\times s}$
$⟨𝒙,𝒚⟩\equiv {𝒚}^{\ast }𝒙\in ℂ$	$⟨⟨𝑿,𝒀⟩⟩\equiv {𝒀}^{\ast }𝑿\in 𝕊$
$⟨𝒙,𝒚⟩={⟨𝒚,𝒙⟩}^{\ast }$	$⟨⟨𝑿,𝒀⟩⟩={⟨⟨𝒀,𝑿⟩⟩}^{\ast }$
$⟨𝒙a,𝒚⟩=⟨𝒙,𝒚⟩a$	$⟨⟨𝑿A,𝒀⟩⟩=⟨⟨𝑿,𝒀⟩⟩A$
$⟨𝒙,𝒚a⟩=a^{\ast }⟨𝒙,𝒚⟩$	$⟨⟨𝑿,𝒀A⟩⟩=A^{\ast }⟨⟨𝑿,𝒀⟩⟩$
$\Vert 𝒙\Vert \equiv \sqrt{⟨𝒙,𝒙⟩}\in {ℝ}_0^{+}$	$\Vert \left|𝑿\right|\Vert \equiv \sqrt{⟨⟨𝑿,𝑿⟩⟩}\in {𝕊}_0^{+}$
$\Vert 𝒙\left|a\right|\Vert =\Vert 𝒙\Vert \left|a\right|$	$\Vert \left|𝑿\left|A\right|\right|\Vert =\Vert \left|𝑿\right|\Vert \left|A\right|$
$I_n\in {ℂ}^{n\times n}$…identity matrix	${𝓘}_n\in {𝕊}^{n\times n}\simeq {ℂ}^{ns\times ns}$…block identity matrix
${𝒆}_k$…$k$th column of $I$	${𝑬}_k$…$k$th block column of $𝓘$
$\{{𝒗}_1,\dots ,{𝒗}_n\}$  linearly independent	$\{{𝑽}_1,\dots ,{𝑽}_n\}$  linearly independent
$\Vert {\sum }_{i=1}^n{𝒗}_i{c}_i\Vert =0\Rightarrow c_i=0\forall i$	$\Vert \left|{\sum }_{i=1}^n{𝑽}_i{C}_i\right|\Vert \text{singular}\Rightarrow C_i\text{singular}\forall i$
$\{{𝒘}_1,\dots ,{𝒘}_n\}$  orthonormal basis of ${ℂ}^n$	$\{{𝑾}_1,\dots ,{𝑾}_n\}$  orthonormal basis of ${𝕊}^n$
$⟨{𝒘}_i,{𝒘}_j⟩={\delta }_{ij}$	$⟨⟨{𝑾}_i,{𝑾}_j⟩⟩={\delta }_{ij}I$
$𝒙={\sum }_{k=1}^n{𝒘}_k⟨𝒙,{𝒘}_k⟩$,  $\forall 𝒙\in {ℂ}^n$	$𝑿={\sum }_{k=1}^n{𝑾}_k⟨⟨𝑿,{𝑾}_k⟩⟩$,  $\forall 𝑿\in {𝕊}^n$

Equations

$$A\mathbf{x}=\mathbf{b},A\in {\mathbb{C}}^{m\times m},\mathbf{b}\in {\mathbb{C}}^m,$$

$$AX=B,A\in {\mathbb{C}}^{m\times m},B\in {\mathbb{C}}^{m\times s},$$

$${\mathcal{K}}_k(A,R_0)\equiv \text{colspan}\{R_0,A{R}_0,\dots ,A^{k-1}{R}_0\}\text{mbox}with{R}_0=B-A{X}_0.$$

$$X_k=\sum _{i=1}^k{V}_i{D}_i,D_i\in {\mathbb{C}}^{s\times s},$$

$$\hat{A}⟶\left[\begin{array}{cc}A& 0\\ 0& I\end{array}\right]\in {\mathbb{C}}^{ns\times ns},\hat{B}⟶\left[\begin{array}{c}B\\ 0\end{array}\right]\in {\mathbb{C}}^{ns\times s},\hat{X}⟶\left[\begin{array}{c}X\\ 0\end{array}\right]\in {\mathbb{C}}^{ns\times s},$$

$$\underbrace{\left[{\footnotesize\begin{array}[]{c|c|c}\begin{matrix}a_{1,1}&a_{1,2}\\ a_{2,1}&a_{2,2}\\ \end{matrix}&\begin{matrix}a_{1,3}&a_{1,4}\\ a_{2,3}&a_{2,4}\\ \end{matrix}&\begin{matrix}a_{1,5}&a_{1,6}\\ a_{2,5}&a_{2,6}\\ \end{matrix}\\ \hline\cr\begin{matrix}a_{3,1}&a_{3,2}\\ a_{4,1}&a_{4,2}\\ \end{matrix}&\begin{matrix}a_{3,3}&a_{3,4}\\ a_{4,3}&a_{4,4}\\ \end{matrix}&\begin{matrix}a_{3,5}&a_{3,6}\\ a_{4,5}&a_{4,6}\\ \end{matrix}\\ \hline\cr\begin{matrix}a_{5,1}&a_{5,2}\\ a_{6,1}&a_{6,2}\\ \end{matrix}&\begin{matrix}a_{5,3}&a_{5,4}\\ a_{6,3}&a_{6,4}\\ \end{matrix}&\begin{matrix}a_{5,5}&a_{5,6}\\ a_{6,5}&a_{6,6}\\ \end{matrix}\end{array}}\right]}_{\left[\begin{array}[]{cccc}A_{1,1}&A_{1,2}&A_{1,3}\\ A_{2,1}&A_{2,2}&A_{2,3}\\ A_{3,1}&A_{3,2}&A_{3,3}\end{array}\right]}\underbrace{\left[{\footnotesize\begin{array}[]{cc}x_{1,1}&x_{1,2}\\ x_{2,1}&x_{2,2}\\ \hline\cr x_{3,1}&x_{3,2}\\ x_{4,1}&x_{4,2}\\ \hline\cr x_{5,1}&x_{5,2}\\ x_{6,1}&x_{6,2}\end{array}}\right]}_{\left[\begin{array}[]{c}X_{1}\\ X_{2}\\ X_{3}\end{array}\right]}=\underbrace{\left[{\footnotesize\begin{array}[]{cc}b_{1,1}&b_{1,2}\\ b_{2,1}&b_{2,2}\\ \hline\cr b_{3,1}&b_{3,2}\\ b_{4,1}&b_{4,2}\\ \hline\cr b_{5,1}&b_{5,2}\\ b_{6,1}&b_{6,2}\end{array}}\right]}_{\left[\begin{array}[]{c}B_{1}\\ B_{2}\\ B_{3}\end{array}\right]},$$

$$|A|\prec |B|$$

$$\mathcal{A}\mathcal{V}=\mathcal{V}\mathcal{H},\mathcal{H}\in {\mathbb{S}}^{n\times n}\text{mbox}blockupperHessenberg,$$

$${\mathbf{X}}_k=\left[\begin{array}{ccc}{\mathbf{V}}_1& \cdots & {\mathbf{V}}_k\end{array}\right]{\mathbf{Y}}_k,$$

$${\mathbf{Y}}_k^F={\left({\mathcal{H}}^{(k)}\right)}^{-1}{\mathbf{E}}_1|||{\mathbf{R}}_0|||.$$

$${\mathbf{Y}}_k^G={\left({\underline{\mathcal{H}}}^{(k)}\right)}^{†}{\mathbf{E}}_1|||{\mathbf{R}}_0|||,$$

$${\mathbf{Y}}_k^F={\mathrm{argmin}}_{\mathbf{Y}}\left\{\Vert {{\mathbf{E}}_{\mathbf{k}}}^T\mathbf{Y}{\Vert }_F;\mathbf{Y}={\left({\mathcal{H}}^{(k)}\right)}^{-}{\mathbf{E}}_1|||{\mathbf{R}}_0|||\right\},$$

$$\mathbf{V}=\left[\begin{array}{c}{V}_1\\ V_2\end{array}\right]⟶\mathcal{Q}\mathbf{V}=\left[\begin{array}{c}\tilde{V}\\ 0\end{array}\right],$$

$$\mathcal{Q}=\left[\begin{array}{cc}\bar{C}& \bar{S}\\ -S& C\end{array}\right]:=\left[\begin{array}{cc}X{Z}^{\ast }& X\\ -Y& YZ\end{array}\right],$$

$$Z=V_1{V}_2^{-1},X^{\ast }X=(I+Z^{\ast }Z{)}^{-1},\text{mbox}and{Y}^{\ast }Y=(I+Z{Z}^{\ast }{)}^{-1},$$

$$X={\left(\sqrt{(I+Z^{\ast }Z)}\right)}^{-\ast }\text{mbox}andY=\sqrt{(I+Z{Z}^{\ast }{)}^{-1}},$$

$${\mathcal{G}}^{(k)}=\left[\begin{array}{cccc}{\mathcal{I}}_{k-1}& & & \\ & {\bar{C}}_k& {\bar{S}}_k& \\ & -S_k& C_k& \\ & & & {\mathcal{I}}_{n-k-1}\end{array}\right],$$

$$\left[\begin{array}{cc}\overline{C_k}& \overline{S_k}\\ -S_k& C_k\end{array}\right]=\left[\begin{array}{cc}{\sqrt{(I+Z_k^{\ast }{Z}_k)}}^{-\ast }{Z}_k^{\ast }& {\sqrt{(I+Z_k^{\ast }{Z}_k)}}^{-\ast }\\ -\sqrt{(I+Z_k{Z}_k^{\ast }{)}^{-1}}& \sqrt{(I+Z_k{Z}_k^{\ast }{)}^{-1}}{Z}_k\end{array}\right],Z_k=H_{k,k}^{(k-1)}{H}_{k+1,k}^{-1},$$

$$\mathcal{Q}={\mathcal{G}}^{(n-1)}{\mathcal{G}}^{(n-2)}\cdots {\mathcal{G}}^{(1)}.$$

$$\left[\begin{array}{cc}\overline{C_k}& \overline{S_k}\\ -S_k& C_k\end{array}\right]\left[\begin{array}{c}{H}_{k,k}^{(k-1)}\\ H_{k+1,j}\end{array}\right]=\left[\begin{array}{c}\Xi \\ 0\end{array}\right]$$

$$⟹H_{k+1,k}={\bar{S}}_k^{\ast }\Xi$$

$$⟹\Xi ={\bar{S}}_k^{-\ast }{H}_{k+1,k}.$$

$$\left[\begin{array}{cc}{\bar{C}}_k& {\bar{S}}_k\\ -S_k& C_k\end{array}\right]\left[\begin{array}{c}|||{\mathbf{R}}_{k-1}^G|||\\ 0\end{array}\right]=\left[\begin{array}{c}\ast \\ -S_k|||{\mathbf{R}}_{k-1}^G|||\end{array}\right]$$

$$|||{\mathbf{R}}_k^G|||=S_k|||{\mathbf{R}}_{k-1}^G|||,\text{mbox}whichimplies|||{\mathbf{R}}_k^G|||=S_k\cdots S_1|||{\mathbf{R}}_0|||.$$

$$|||{\mathbf{R}}_k^F|||$$

$$=\left|H_{k+1,k}{\left(H_{k,k}^{(k-1)}\right)}^{†}\right||||{\mathbf{R}}_{k-1}^G|||$$

$$=|C_k^{†}{S}_k||||{\mathbf{R}}_{k-1}^G|||=|C_k^{†}|S_k|||{\mathbf{R}}_{k-1}^G|||=|C_k^{†}||||{\mathbf{R}}_k^G|||,$$

$$⟨⟨{\mathbf{R}}_k^F,{\mathbf{R}}_k^F⟩{⟩}^{†}=⟨⟨{\mathbf{R}}_k^G,{\mathbf{R}}_k^G⟩{⟩}^{-1}-⟨⟨{\mathbf{R}}_{k-1}^G,{\mathbf{R}}_{k-1}^G⟩{⟩}^{-1}$$

$$⟨⟨{\mathbf{R}}_k^G,{\mathbf{R}}_k^G⟩{⟩}^{-1}=\sum _{i=0}^k⟨⟨{\mathbf{R}}_i^F,{\mathbf{R}}_i^F⟩{⟩}^{†}.$$

$$⟨⟨{\mathbf{R}}_k^F,{\mathbf{R}}_k^F⟩⟩$$

$$=|||{\mathbf{R}}_k^G||{|}^{\ast }{\left(I-|||{\mathbf{R}}_k^G||||||{\mathbf{R}}_{k-1}^G||{|}^{-1}|||{\mathbf{R}}_{k-1}^G||{|}^{-\ast }|||{\mathbf{R}}_k^G||{|}^{\ast }\right)}^{†}|||{\mathbf{R}}_k^G|||.$$

$$⟨⟨{\mathbf{R}}_k^F,{\mathbf{R}}_k^F⟩{⟩}^{†}$$

$$=|||{\mathbf{R}}_k^G||{|}^{-1}|||{\mathbf{R}}_k^G||{|}^{-\ast }-|||{\mathbf{R}}_{k-1}^G||{|}^{-1}|||{\mathbf{R}}_{k-1}^G||{|}^{-\ast }$$

Full text

{AMS}

65F10, 65F15

The second author would like to dedicate this work to the late Richard Timoney, who was very generous with his time in explaining the basics of ∗-algebras to an uninitiated new colleague.

1 Introduction

The celebrated GMRES algorithm [21] is an effective, widely-used iterative method for solving linear systems

[TABLE]

with non-Hermitian coefficient matrices. It has been shown that, in contrast to the case of Hermitian matrices, for non-Hermitian coefficient matrices, we cannot guarantee certain convergence speed based solely on our knowledge of the spectrum. In particular, for any given non-increasing sequence of positive numbers $f_{0}\geq f_{1}\geq\cdots\geq f_{n-1}$, one can construct a non-Hermitian matrix with arbitrary spectrum for which GMRES produces residuals whose norms correspond to this sequence [15]. In [19], it was observed this does not mean eigenvalues are meaningless for the convergence of GMRES applied to non-normal problems. This result simply establishes one extreme of what role eigenvalues can play in the residual convergence of GMRES.

If instead we solve multiple systems with the same $A$ (i.e., a system with multiple right-hand sides)

[TABLE]

a direct generalization of GMRES called block GMRES exists which produces approximate solutions simultaneously for all right-hand sides. We would like to explore if block GMRES residual convergence admits a similar characterization as in [15].

Most often, analysis of block GMRES takes the view of the method as a minimization of each individual residual over a sum of spaces. This enables some basic convergence analysis but such analysis fails to capture the full picture of block GMRES behavior, which is influenced by the interaction between the different right-hand sides. The nature of this interaction is quite difficult to describe when considering block GMRES as a method treating a collection of scalar linear systems.

Indeed, some authors have taken a different approach, discussing these methods in terms of vector blocks in $\mathbb{C}^{m\times s}$, namely [11] and [24]. We demonstrate here that embracing a totally block view of this iteration greatly simplifies analysis of block GMRES and allows us to obtain clean convergence results.

The aim of the paper is to extend well-known GMRES/Arnoldi convergence results developed in [1, 7, 16, 15] (along with many excellent follow-up papers by subsets of the same authors) to the block case using the framework of [11, 12].

The paper is organized as follows. Section 2 introduces the block methods that will be analyzed and the framework used to perform the analysis. In Section 3, we generalize the notion of non-increasing convergence curve to the block setting. Section 4 provides a characterization of matrices and starting vectors exhibiting prescribed convergence behavior, including discussion of the spectral properties of the obtained coefficient matrices. We summarize our results and formulate open questions in Section 5. Throughout the paper, we assume exact arithmetic. We also introduce some specialized notation for this setting in Section 2.1.

2 Preliminaries

In each step, block Krylov subspace methods look for an approximation of each individual solution in the space

[TABLE]

For simplicity of presentation, we assume $X_{0}=0$ leading to $R_{0}=B$ throughout the paper.

If the columns of the block vectors $V_{1},\ldots,V_{k}$, $V_{j}\in\mathbb{C}^{m\times s}$, form a basis of the subspace $\mathcal{K}_{k}(A,B)$, we say that $V_{1},\ldots,V_{k}$ is a block basis of $\mathcal{K}_{k}(A,B)$. The block solution $X_{k}$ can be represented as a block linear combination of these block vectors:

[TABLE]

and we say in this setting that $X_{k}\in\text{blockspan}\{B,AB,\ldots,A^{k-1}B\}$. The particular choice of $\{D_{i}\}_{i=1}^{k}$ is defined by the conditions which the method imposes on the residual $R_{k}=B-AX_{k}$; cf. Section 2.2. For an overview of block Krylov subspace methods, see, e.g., [20, sec. 6.12].

Comparing Eq. 4 with the standard definition of a linear combination of a set of vectors, we see that the $s\times s$ matrices here play the role of complex numbers $\mathbb{C}$. Let $n=\lceil\frac{m}{s}\rceil$. Padding the matrix and the right-hand side by zeros as follows

[TABLE]

the unknown $X$ and block GMRES approximations thereof do not change and we can view the new matrix also as an $n\times n$ array of $s\times s$ matrices, and similarly the new right-hand side and the solution as a vector of length $n$ of $s\times s$ matrices. Thus, to simplify the presentation and without loss of (except, perhaps, at the final iteration), we assume hereafter that $m=ns$, and decompose $A$, $B$, and $X$ into blocks, which we demonstrate on the following example.

Let $m=6$ and $s=2$, then $n=3$, and we can write

[TABLE]

where $A_{i,j}$, $X_{j}$, and $B_{i}$ are all $2\times 2$ matrices.

In this paper, we restrict our analysis to the case when $\dim(\mathcal{K}_{n}(A,B))=ns$.

2.1 Framework

To generalize the results from standard Krylov subspace methods to block ones, we follow [11, sec. 2] and replace $\mathbb{C}$ by the non-commutative ∗-algebra $\mathbb{S}$ of complex $s\times s$ matrices.111The authors in [11] considered various other subrings of $\mathbb{S}$ to classify different block-type Krylov subspace methods in a common framework. We observe that objects and operations over $\mathbb{C}$ not relying on commutativity have counterparts in $\mathbb{S}$ that are relevant for the analysis of block Krylov subspace methods; some of them are shown in Table 1. We use the notation $\mathbb{S}\simeq\mathbb{C}^{s\times s}$ to say that the elements of $\mathbb{S}$ are from $\mathbb{C}^{s\times s}$ but equipped with additional operations.

We emphasize that to comply with the standard block Krylov subspace method notation, the generalization of positive real numbers denoted by $\mathbb{S}^{+}$ corresponds to the upper triangular matrices with positive diagonal entries.222The standard subalgebra of positive entries is represented by Hermitian positive-definite matrices, whose Cholesky factors have a one-to-one correspondence with elements of $\mathbb{S}^{+}$ used here. This choice becomes particularly important for the ∗-algebra formulation of the block Arnoldi algorithm, see Section 2.2. Note also that in $\mathbb{S}^{n}$, similarly to $\mathbb{C}^{n}$, we distinguish between the block absolute value $|\cdot|$ and the block norm $|||\cdot|||$. Both these operators map to $\mathbb{S}_{0}^{+}$, however the first is defined only for the elements of $\mathbb{S}$, whereas the second one is defined for an arbitrary vector with elements in $\mathbb{S}$, see Table 1.

{notation}

In order to keep our notation clear, we denote matrices over $\mathbb{S}$ with caligraphic letters, vectors over $\mathbb{S}$ by bold capital letters, and elements of $\mathbb{S}$ by unbolded capital letters. Vectors over $\mathbb{C}$ are denoted by lower-case bold letters and scalar elements of $\mathbb{C}$ are denoted by lower-case unbolded letters.

In order to be able to investigate convergence behavior, we generalize the ordering of nonnegative real numbers, denoted by $\mathbb{R}_{0}^{+}$, to that of upper triangular complex matrices with real nonnegative diagonal entries, denoted by $\mathbb{S}_{0}^{+}$, as follows:

[TABLE]

where $\overset{\text{{\tiny Loewner}}}{\prec}$ and $\overset{\text{{\tiny Loewner}}}{\preceq}$ is the Loewner (partial) ordering of Hermitian matrices.333Defined by the relation $A\overset{\text{{\tiny Loewner}}}{\prec}B\iff B-A$ is Hermitian positive definite and $A\overset{\text{{\tiny Loewner}}}{\preceq}B\iff B-A$ is Hermitian positive semi-definite; see [18]. Most of the results presented in the upcoming sections can be formulated also in terms of the decay of the “squares”, i.e., the traditional Loewner order. To make the results formally as similar as possible to the classical results for standard GMRES, we have chosen to use the non-increasing order of the norms rather than their squares.

2.2 Arnoldi-based methods in the new framework

Every Krylov subspace method relies on construction of a well-conditioned (ideally orthonormal) basis of the Krylov subspace. For standard GMRES, this orthonormal basis is computed using the Arnoldi algorithm. For block GMRES, where our aim is to compute an orthonormal block basis, see the last part of Table 1, we use the block Arnoldi algorithm, iteratively producing orthonormal block Arnoldi vectors $\bm{V}_{1},\ldots,\bm{V}_{k}$, $k=1,\ldots,n,$ spanning the Krylov subspace Eq. 3. If we assume no premature breakdown, it can be carried out for $n$ iterations to produce a basis for the space $\mathbb{S}^{n}$. In Algorithm 1, following [11], we provide a pseudocode using the notation introduced in Table 1.

There are various strategies for how to handle rank-deficient Arnoldi vectors; see, e.g., [2, 10, 25]. Hereafter, we will however assume that there is no breakdown in Algorithm 1.

Taken to the $n$th iteration, the block Arnoldi algorithm yields the block Arnoldi relation

[TABLE]

with $\bm{\mathcal{V}}$ having as columns an orthonormal basis of $\mathbb{S}^{n}$, and the relation Eq. 8 represents a full orthogonal Hessenberg factorization of the matrix $\bm{\mathcal{A}}$. Since the blocks on the subdiagonal of the block Hessenberg matrix $\bm{\mathcal{H}}$ are from $\mathbb{S}^{+}$, our definition of $\mathbb{S}^{+}$ assures that the matrix $\bm{\mathcal{H}}$ has exactly $s$ subdiagonals and the entries of the last one are positive. Therefore, the block Arnoldi algorithm formulated through the new framework computes the same Arnoldi decomposition as the block Arnoldi algorithm performed in a conventional way; see, e.g., [20, Section 6.12]. Terminating Algorithm 1 after iteration $k$ produces a basis of the Krylov subspace Eq. 3.

Using the block Arnoldi method (blArnoldi) for eigenvalue problems, the eigenvalues of $\bm{\mathcal{A}}$ are in each step $k$ approximated by the Ritz values, i.e., the $sk$ eigenvalues of the $k$th principal submatrix $\bm{\mathcal{H}}^{(k)}\in\mathbb{S}^{k\times k}$ of $\bm{\mathcal{H}}$.444By the eigenvalue of a matrix $\bm{\mathcal{A}}$ we mean the standard eigenvalue, i.e., $\lambda\in\mathbb{C}$ satisfying $\bm{\mathcal{A}}\bm{v}=\lambda\bm{v}$ for some nonzero $\bm{v}\in\mathbb{C}^{ns}$. We note that considering the linear system over $\mathbb{S}$ induces a notion of a block eigenvalue decomposition, i.e., there exist $\bm{\Lambda}\in\mathbb{S}$ such that $\bm{\mathcal{A}}\bm{V}=\bm{V}\bm{\Lambda}$ for some $\bm{V}\in\mathbb{S}^{n}$ with linearly independent columns. We do not treat such a notion in this paper as it is highly technical, due in part to the non-commutitivity of multiplication in $\mathbb{S}$. Furthermore, as we show in Section 4.5, the correct notion of block eigenvalues in this setting is that of block solvents of specified $\lambda$-matrices. For systems of linear algebraic equations, the approximate block solutions $\bm{X}_{k}$ are constructed as a linear combination of the computed orthonormal basis of $\mathcal{K}_{k}(\bm{\mathcal{A}},\bm{B})$, i.e.,

[TABLE]

where the block vector $\bm{Y}_{k}\in\mathbb{S}^{k}$ is obtained from small projected problems. We analyze two different methods: block FOM (blFOM) and block GMRES (blGMRES). The blFOM method is a Galerkin method, keeping the individual residuals in each step orthogonal to all previous, i.e.,

[TABLE]

The blGMRES method minimizes the Euclidean norm of each of the individual residuals555The application of the Moore-Penrose pseudoinverse to each column individually gives the least-squares solution for each right-hand side, and applying it to each column individually is equivalent to applying it to the block right-hand side. Further, minimizing the Euclidean norm of individual residuals is also equivalent to minimizing the Frobenius norm of the block residual. This was utilized, for example, in the analysis of blGMRES in [24]., i.e.,

[TABLE]

where $\underline{\bm{\mathcal{H}}}^{(k)}\in\mathbb{S}^{(k+1)\times k}$ is the upper-left $(k+1)\times k$ block of $\bm{\mathcal{H}}$ and $\cdot^{\dagger}$ denotes here and hereafter the Moore-Penrose pseudoinverse.

While the blGMRES solution always exists, the blFOM solution is not defined when the matrix $\bm{\mathcal{H}}^{(k)}$ is singular. In this case, we will consider the generalized blFOM solution and define

[TABLE]

where $\cdot^{-}$ denotes a generalized inverse; see also [26]. Since the true and the generalized FOM solution coincide for nonsingular $\bm{\mathcal{H}}^{(k)}$, we denote both by $\bm{X}_{k}^{F}$ to simplify the notation.

3 Admissible convergence behavior of blGMRES

In this section, we provide a block definition of admissible convergence behavior of the norm-minimizing method by generalizing some of the well-known relations between the residuals of the norm-minimizing and the Galerkin (residual orthogonalizing) method.

3.1 Block Givens transformation

Each step of standard GMRES requires computation of one new elementary rotation to eliminate the $k$th subdiagonal entry of $\bm{\mathcal{H}}$. These rotations facilitate analysis of GMRES convergence and its relationship to FOM [20, sec. 6.5.7]. For blGMRES, a similar analysis is available using the product of Householder transformations [26], but a block analog to Givens rotations provides additional and clearer results.

For the elimination in blGMRES, $s^{2}$ standard elementary Givens rotations are needed [20, Section 6.12]. A product of (elementary) Givens rotations is an orthogonal transformation but loses the properties of rotation, except for the very special case when all the individual rotations commute; see [3, 14, 6]. Therefore, the block Givens transformation will not be a true generalization of a Givens rotation, but rather a block representation of the product of the $s^{2}$ standard elementary Givens rotations.666We emphasize that here we are concerned with mathematical properties of such block Givens transformations. In practical computations, the individual Givens sines and cosines are stored, and the product is rarely computed explicitly. 777When $s$ is large, the subdiagonal entries of the block Hessenberg matrix are eliminated using Householder reflections. There are two main possible generalizations of Householder reflections to the block case. The first, see [22], preserves the properties of a reflection, but is only able to eliminate the subdiagonal block. The second one, see [23], is able to eliminate all subdiagonal entries and is a block representation of a product of $s$ standard elementary Householder reflections, but does not have properties of a reflection. In any case, they are of little use when generalizing relations from the standard case, where Givens rotations always are performed.

We follow the idea of [17]. Assume first a general orthogonal transformation $\bm{\mathcal{Q}}$ that eliminates an entry of a block vector, i.e.,

[TABLE]

where $V_{1},V_{2}\in\mathbb{S}$ and $|\widetilde{V}|=|||\bm{V}|||$. If $V_{2}$ is invertible, then the unitary matrix $\bm{\mathcal{Q}}$ eliminating the block $V_{2}$ can be constructed as

[TABLE]

where

[TABLE]

which can be verified by simple computation.

Note that there is freedom in the choice of $X$ and $Y$, since they can be arbitrary right factors of the matrices $(I+Z^{*}Z)^{-1}$ and $(I+ZZ^{*})^{-1}$, respectively. This allows us also to control the non-zero pattern of the matrix $\bm{\mathcal{Q}}$. If we take

[TABLE]

then, due to our definition of $\sqrt{\,\cdot\,}$ , the matrix $\bm{\mathcal{Q}}$ in Eq. 14 becomes $(2s+1)$-diagonal. In this case, we can see $\bm{\mathcal{Q}}$ as a product of individual standard Givens rotations applied to $\bm{V}$ to eliminate the subdiagonal entries. If $V_{1},V_{2}\in\mathbb{C}$, the choice Eq. 16 will lead to the standard elementary Givens rotation with $\bar{S}$ and $\bar{C}$ being the complex conjugates of $C$ and $S$, respectively. This is however generally not the case for $V_{1},V_{2}\in\mathbb{S}$, i.e., $\bar{C}\neq C^{*}$ and $\bar{S}\neq S^{*}$. The following proposition explains why the choice Eq. 16 is important for the solution of Eq. 11 by factorizing the block upper Hessenberg matrix $\bm{\mathcal{H}}$ from Eq. 8.

Proposition 3.1**.**

Let $\bm{\mathcal{H}}\in\mathbb{S}^{n\times n}$ be an upper block Hessenberg matrix with blocks $H_{i,j}$, and let further the subdiagonal blocks be upper triangular with positive diagonal entries, i.e., $H_{k+1,k}\in\mathbb{S}^{+}$, $k=1,\ldots,n-1$. Define, for $k=1,\ldots,n-1$, the block Givens transformations $\bm{\mathcal{G}}^{(k)}$ as

[TABLE]

where

[TABLE]

with $\bm{\mathcal{H}}^{(0)}=\bm{\mathcal{H}}$ and $\bm{\mathcal{H}}^{(k)}=\bm{\mathcal{G}}^{(k)}\bm{\mathcal{H}}^{(k-1)}$. Define further the unitary matrix $\bm{\mathcal{Q}}$ as

[TABLE]

*Then $\bm{\mathcal{QH}}$ is a block upper triangular matrix, with the first $n-1$ block diagonal entries in $\mathbb{S}^{+}$. *

Proof 3.2**.**

The unitarity of $\bm{\mathcal{Q}}$ and the upper block triangular form of $\bm{\mathcal{QH}}$ follows directly from the construction and discussion above. Further, observe that

[TABLE]

*Since $\mathbb{S}^{+}$ is a multiplicative group and both $\bar{S}_{k}^{-*}\in\mathbb{S}^{+}$ and $H_{k+1,k}\in\mathbb{S}^{+}$, also $\Xi\in\mathbb{S}^{+}$, which yields the desired statement. *

Proposition 3.1 shows that the $n-1$ block Givens transformations $\bm{\mathcal{G}}_{1},\ldots,\bm{\mathcal{G}}_{n-1}$ defined in Eq. 17, with additional normalization of the last block entry of $\bm{\mathcal{H}}^{(n-1)}$, enable reduction of the Hessenberg matrix $\bm{\mathcal{H}}$ to the upper triangular form.

3.2 Residual norms and the peak-plateau relation

In this section, we use the block Givens transformations to generalize some of the well-know results about GMRES and FOM residuals. Note that for the generalized blFOM solution Eq. 12, the analysis only addresses the component of the true residual that satisfies the desired Galerkin condition, i.e., $\bm{R}_{k}^{F}:=\bm{V}_{k+1}\bm{V}_{k+1}^{*}\left(\bm{B}-\bm{\mathcal{A}}\bm{X}_{k}^{F}\right).$ Since the blGMRES solution always exists, it is defined in the conventional way as $\bm{R}_{k}^{G}:=\bm{B}-\bm{\mathcal{A}}\bm{X}_{k}^{G}$.

Let $\bm{R}_{0}=\bm{V}_{1}|||\bm{R}_{0}|||$ be the initial residual. Note that since we assume no premature breakdown in Algorithm 1, it holds that $|||\bm{R}^{G}_{k}|||$ is invertible for $k=0,\ldots,n-1$. We follow [20, sec. 6.5.7]. Since the $k$th Givens rotation modifies the right-hand side of the projected problem as

[TABLE]

and since $S_{k}\in\mathbb{S}^{+}$, we get that the blGMRES residual satisfies

[TABLE]

Using Proposition 3.1, we obtain that the (generalized) blFOM residual satisfies

[TABLE]

where $C_{k}^{\dagger}$ becomes $C_{k}^{-1}$, whenever the $k$th principal submatrix of $\bm{\mathcal{H}}$, and therefore also $H^{(k-1)}_{k,k}$, is invertible; cf. Eq. 12.

Since blGMRES and blFOM form a norm-minimizing/Galerkin pair, one expects that they satisfy some form of peak-plateau relation, see, e.g., [20, sec. 6.5.7] or [4]. The following proposition shows that this is indeed the case.

Proposition 3.3**.**

The residuals of blGMRES and blFOM satisfy

[TABLE]

and

[TABLE]

Proof 3.4**.**

Using Eq. 25, we have

[TABLE]

By taking the pseudoinverse, we obtain

[TABLE]

*which gives Eq. 28. Applying relation Eq. 28 recursively gives the second statement. *

Relation Eq. 28 can be viewed as a generalization of the peak-plateau relation, to which it reduces when $s=1$. For other relations between blGMRES and blFOM, see also [26].

Also note that $\langle\langle\bm{R}_{k}^{F},\bm{R}_{k}^{F}\rangle\rangle$ is singular if and only if $|||\bm{R}_{k}^{F}|||$ is singular, i.e., if either the standard FOM iterate does not exist for some of the right-hand sides, or if the individual FOM residuals are linearly dependent. This corresponds to the situation when the residual update in blGMRES is of rank smaller than $s$, see [26, p. 173]. Such situation will be of interest also in Section 4.3.

3.3 Admissible convergence behavior of blGMRES

The results of the previous section have some nontrivial consequences for the convergence behavior of blGMRES.

Theorem 3.5**.**

The blGMRES residuals satisfy

[TABLE]

Proof 3.6**.**

Since $\langle\langle\bm{R}_{k}^{F},\bm{R}_{k}^{F}\rangle\rangle$ is a positive semidefinite matrix, relation Eq. 28 directly implies that for $k=1,\ldots,n$

[TABLE]

Using the generalization of the Loewner ordering from Section 2 we can write that

[TABLE]

*Applying this relation recursively, we obtain Eq. 35. *

We call the sequence $|||\bm{R}_{0}|||,|||\bm{R}_{1}|||,\ldots,|||\bm{R}_{n-1}|||$ of elements in $\mathbb{S}^{+}$ an admissible convergence behavior of blGMRES, if it satisfies Eq. 35.

Remark 3.7**.**

Using $\|\bm{R}_{k}^{G}e_{j}\|^{2}=\langle\bm{R}_{k}^{G}e_{j},\bm{R}_{k}^{G}e_{j}\rangle=e_{j}^{T}\langle\langle\bm{R}_{k}^{G},\bm{R}_{k}^{G}\rangle\rangle e_{j}$, $j=1,\ldots,s$, relation Eq. 35 trivially implies monotonic convergence of the size of the individual residuals. However, relation Eq. 35 is generally stronger, because it takes into account the inter-residual relationships, expressed by the off-diagonal entries of the matrices $\langle\langle\bm{R}_{k}^{G},\bm{R}_{k}^{G}\rangle\rangle$, $k=0,\ldots,n-1$. If two individual initial residuals (of the same norm) are almost linearly dependent, one cannot expect radically different convergence behaviors for each right-hand side. We will demonstrate this with the following example. Let the initial residuals be almost linearly dependent:

[TABLE]

and let the size of the first residual be decreased to $\sqrt{\varepsilon}$ and of the second one to $\sqrt{1-\varepsilon}$:

[TABLE]

*Then there is no $p$ such that $\langle\langle\bm{R}_{0}^{G},\bm{R}_{0}^{G}\rangle\rangle\overset{\text{{\tiny Loewner}}}{\succeq}\langle\langle\bm{R}_{1}^{G},\bm{R}_{1}^{G}\rangle\rangle\overset{\text{{\tiny Loewner}}}{\succ}0$, and therefore such convergence behavior cannot be exhibited by blGMRES. Conversely, if two initial residuals are orthogonal, Eq. 37 does not give any further restriction on the convergence curve of the individual residuals. *

Remark 3.8**.**

*Note that the result of Theorem 3.5 is also very intuitive in the following sense: Any non-increasing sequence (in the Loewner sense) can be generated through the norms of the orthogonal projections of a (block) vector onto a sequence of embedded subspaces. And vice versa, norms of the orthogonal projections of a (block) vector to any sequence of embedded subspaces will generate a non-increasing sequence (in the Loewner sense). Since orthogonal projections of the residuals to the residual Krylov subspace are the very essence of the GMRES method, the Loewner ordering of its residual sizes is its inherent property. *

The Loewner ordering is also natural to consider when block Krylov methods are formulated via polynomial approximation of the inverse of $\bm{\mathcal{A}}$. The polynomial with matrix-valued coefficients applied to the right-hand side $\bm{B}$ to obtain the residual is optimal in the sense of the following theorem. We use the concept of matrix-valued polynomials, which we define in Eq. 111 and elaborate on in Appendix A, to give the proper generalization of the GMRES residual polynomial minimization.

Theorem 3.9**.**

Denote by $\mathcal{P}_{k}$ the set of all $\lambda$-matrices $P$ of degree at most $k$, with coefficients in $\mathbb{S}$, and satisfying $P(0)=I$, see Appendix A. The blGMRES residual norms satisfy

[TABLE]

*where $\circ$ is as in Eq. 120 and the minimum is considered in the generalized Loewner sense Eq. 7. *

Proof 3.10**.**

From the definition of the blGMRES iterate it follows that $\bm{R}_{k}^{G}=\widehat{P}_{k}(\bm{\mathcal{A}})\circ\bm{B}$ for some $\widehat{P}_{k}\in\mathcal{P}_{k}$. Therefore, it suffices to show that by changing the polynomial, one can only make the norm of the resulting block vector larger. Any $\lambda$-matrix in $\mathcal{P}_{k}$ can be written as $\widehat{P}_{k}(\lambda)+\lambda Q_{k-1}(\lambda)$, where $Q_{k-1}$ is a $\lambda$-matrix of degree at most $k-1$. Therefore, we can write

[TABLE]

Further, we observe that $\bm{\mathcal{A}}Q_{k-1}(\bm{\mathcal{A}})\circ\bm{B}\in\bm{\mathcal{A}}\mathcal{K}_{k}(\bm{\mathcal{A}},\bm{B})$. Since the residual $\bm{R}^{G}_{k}$ is orthogonal to this space, we obtain

[TABLE]

With this, Eq. 41 becomes

[TABLE]

*which finishes the proof. *

In the following section we show that any admissible convergence behavior Eq. 35 is actually attainable by blGMRES. For better readability, we drop the superscript $\cdot^{G}$ in the remainder of this paper.

4 Prescribing convergence behavior

In this section, we utilize the new framework for block Krylov subspace methods and generalize some of the results of [1, 7, 16, 15] for standard Arnoldi and GMRES to the block case. We show that, under moderate conditions, we simultaneously can prescribe the residual convergence of blGMRES($\bm{\mathcal{A}}$,$\bm{B}$) as well as the spectral properties of $\bm{\mathcal{A}}$ and the principal submatrices $\bm{\mathcal{H}}^{(k)}$ of the ultimate upper Hessenberg matrix $\bm{\mathcal{H}}$ produced by blArnoldi($\bm{\mathcal{A}}$,$\bm{B}$). These entities will be prescribed based on the framework introduced in Section 2.1 and have a slightly different form than in the standard case.

To prescribe the blGMRES residual norms, we require

[TABLE]

for some given sequence $\{F_{k}\}_{k=0}^{n-1}$, $F_{k}\in\mathbb{S}^{+}$, satisfying the admissibility condition, with

[TABLE]

The spectral properties of $\bm{\mathcal{A}}$ are in the standard case prescribed through the similarity to the companion matrix corresponding to the characteristic polynomial defined by these eigenvalues; see [15]. We will show that by considering block companion matrices of the form

[TABLE]

we can proceed analogously in the block case. Note that the matrix $\bm{\mathcal{C}}$ defined in Eq. 49 is the block companion matrix corresponding to the $\lambda$-matrix

[TABLE]

We will be looking for the matrices $\bm{\mathcal{A}}$ and the right-hand sides $\bm{B}$ annihilating the polynomial $M$, i.e.,

[TABLE]

where the operation $\circ$ is as in Eq. 120. Using $\bm{\mathcal{A}}=\bm{\mathcal{V}}\bm{\mathcal{H}}\bm{\mathcal{V}}^{*}$ and $\bm{B}=\bm{V}_{1}|||\bm{R}_{0}|||$, Eq. 51 can be rewritten as

[TABLE]

The relation between the spectral properties of $\bm{\mathcal{A}}$ and condition Eq. 51 is discussed in more detail in Section 4.5.

Due to [7], it is known that in the standard case, not only the eigenvalues of the matrix $\bm{\mathcal{A}}$ and therefore also the ultimate upper Hessenberg matrix $\bm{\mathcal{H}}$ can be prescribed, but also the eigenvalues of all the principal block submatrices $\bm{\mathcal{H}}^{(k)}$ of this Hessenberg matrix, i.e., the Ritz values. In the block case, spectral properties of the submatrices $\bm{\mathcal{H}}^{(k)}$ will be enforced analogously to the spectral properties of $\bm{\mathcal{A}}$ through

[TABLE]

where $M^{(k)}$ is defined as

[TABLE]

4.1 Attaining prescribed spectral properties of $\bm{\mathcal{A}}$

Defining the Krylov matrix

[TABLE]

condition Eq. 51 is equivalent to

[TABLE]

yielding

[TABLE]

Consider the unique $QR$-decomposition (with separate diagonal scaling matrix $\bm{\mathcal{D}}$) of the matrix $\bm{\mathcal{K}}$,

[TABLE]

with $\bm{\mathcal{V}}$ unitary, $\bm{\mathcal{D}}$ non-singular block diagonal with diagonal entries in $\mathbb{S}^{+}$, and $\bm{\mathcal{U}}$ non-singular upper triangular with unit diagonal entries. With this decomposition, equations Eq. 57 become

[TABLE]

Similar to the standard case, see, e.g., [19, sec. 3], any product of the form $\bm{\mathcal{D}}\bm{\mathcal{U}}\bm{\mathcal{C}}\bm{\mathcal{U}}^{-1}\bm{\mathcal{D}}^{-1}$ is a block upper Hessenberg matrix with subdiagonal entries in $\mathbb{S}^{+}$ and vice versa. Moreover, $\bm{\mathcal{DU}}$ satisfies

[TABLE]

where $D_{1}$ is the first block entry of $\bm{\mathcal{D}}$. 888Note that contrary to the standard case, for a given upper Hessenberg matrix $\bm{\mathcal{H}}$, the matrices $\bm{\mathcal{U}}$ and $\bm{\mathcal{C}}$ are not defined uniquely, since they depend on the choice of $D_{1}$. For the same reason, the scaling by $|||\bm{B}|||$ cannot be easily omitted in Eq. 53.

Using the decomposition Eq. 59, we provide a complete characterization of matrices and right-hand sides providing prescribed convergence behavior. The convergence behavior of blArnoldi and blGMRES is unitarily invariant; therefore the choice of $\bm{\mathcal{V}}$ plays no role in the analysis. Note that from Eq. 59, we already have that $D_{1}=|||\bm{B}|||$. In the following sections, we determines the $(n-1)(n-2)/2$ block entries of $\bm{\mathcal{U}}$ and the last $n-1$ diagonal block entries of $\bm{\mathcal{D}}$ so that the prescribed blGMRES and blArnoldi behavior is met.999In [15, 1], the characterization of matrices and right-hand sides providing prescribed behavior is obtained through the factorization of the matrix corresponding to the Krylov residual subspaces $\bm{\mathcal{A}}\mathcal{K}_{n}(\bm{\mathcal{A}},\bm{B})$. This factorization provides alternative formulation of the results presented in the subsequent sections.

4.2 Attaining prescribed blArnoldi convergence

We show that the blArnoldi convergence Eq. 53 is encoded solely in the upper triangular matrix $\bm{\mathcal{U}}$.

Using the definition of $M^{(k)}(\lambda)$, condition Eq. 53 can be rewritten using Eq. 54 as

[TABLE]

Since $\bm{\mathcal{H}}=\bm{\mathcal{D}}\bm{\mathcal{U}}\bm{\mathcal{C}}\bm{\mathcal{U}}^{-1}\bm{\mathcal{D}}^{-1}$ and $\bm{\mathcal{U}}^{-1}\bm{\mathcal{D}}^{-1}\bm{E}_{1}|||\bm{B}|||=\bm{E}_{1}$, we observe that

[TABLE]

where $\bm{\mathcal{D}}^{(k)}$ and $\bm{\mathcal{U}}^{(k)}$ are the $k$th principal submatrices of $\bm{\mathcal{D}}$ and $\bm{\mathcal{U}}$, respectively. Substituting into Eq. 61, we have, for $k=1,\ldots,n-1$,

[TABLE]

Using the factorization of $\bm{\mathcal{H}}$ and the structure of $\bm{\mathcal{D}}$, $\bm{\mathcal{U}}$, and $\bm{\mathcal{C}}$, the left-hand side of Eq. 63 can be rewritten as

[TABLE]

Furthermore, using the block inversion formula and the fact that that $\bm{\mathcal{U}}$ is upper triangular with identity matrices on the main diagonal, we have

[TABLE]

Substituting Eqs. 64 and 70 back into Eq. 63 and premultiplying both sides by the inverse of $\bm{\mathcal{D}}^{(k)}\bm{\mathcal{U}}^{(k)}$, we finally obtain

[TABLE]

Conditions Eq. 71 uniquely determine the upper triangular matrix $\bm{\mathcal{U}}$ as

[TABLE]

4.3 Attaining prescribed blGMRES convergence

We show that for a given $\bm{\mathcal{U}}$, the prescribed blGMRES convergence Eq. 47 can be achieved by a proper choice of the block diagonal matrix $\bm{\mathcal{D}}$, which also defines the subdiagonal entries of $\bm{\mathcal{H}}$.

Let $\bm{\mathcal{W}}$ be any matrix such that its columns $\bm{W}_{1},\ldots,\bm{W}_{k}$ form an orthonormal basis of $\bm{\mathcal{A}}\mathcal{K}_{k}(\bm{\mathcal{A}},\bm{B})$, $k=1,\ldots,n$. To satisfy Eq. 47, it has to hold that

[TABLE]

Now we relate the columns of $\bm{\mathcal{W}}$ to those of $\bm{\mathcal{V}}$. First, we see that

[TABLE]

where $\bm{\mathcal{Q}}\in\mathbb{S}^{n\times n}$ is unitary block diagonal and ${\bm{\mathcal{R}}},\widehat{\bm{\mathcal{R}}}\in\mathbb{S}^{(n-1)\times(n-1)}$ are nonsingular upper block triangular matrices. Combining with Eq. 58, we now have two factorizations of $\bm{\mathcal{K}}$, i.e.,

[TABLE]

The right-hand side of Eq. 76 is in the form of the $QR$-decomposition; we study that of the left-hand side by looking at the structure of the Cholesky factorization of $\bm{\mathcal{G}}^{*}\bm{\mathcal{G}}$. We observe that

[TABLE]

where $\bm{\mathcal{R}}_{G}$ is the unique upper triangular Cholesky factor of $(\bm{\mathcal{I}}-\widehat{\bm{G}}F_{0}^{-1}F_{0}^{-*}\widehat{\bm{G}}^{*})$, which is positive definite since $|||\widehat{\bm{G}}F_{0}^{-1}|||\preceq I$. There must be equality between the $R$-factors of the unique $QR$-decomposition of both sides of Eq. 76. Note that the matrix $\widehat{\bm{\mathcal{R}}}$ has arbitrary nonsingular block entries on the diagonal. To obtain the unique $R$-factor with entries in $\mathbb{S}^{+}$, further transformation by a block diagonal unitary matrix (here denoted by $\widehat{\Gamma}$) is needed. Using convenient decompositions of the involved matrices, we obtain from Eq. 76 that

[TABLE]

see also [7, sec. 3] for the analog in the standard case.

We proceed by comparing individual block entries. Equality of the first diagonal block entries gives

[TABLE]

which is satisfied trivially. Equality of the second diagonal block entries gives an expression for $\widehat{\bm{\mathcal{R}}}$

[TABLE]

Substituting from Eq. 81 and Eq. 82 to the equation given by the off-diagonal entry, we obtain

[TABLE]

We now investigate the objects in equation Eq. 83. First, using the block inversion formula on Eq. 72, we have

[TABLE]

Second, we apply Lemma B.1 to $\widehat{\bm{G}}F_{0}^{-1}$, and together with the definition of $\widehat{\bm{G}}$, we obtain

[TABLE]

This substituting to Eq. 83 and defining $\widehat{\bm{\mathcal{Q}}}:=\text{diag}(Q_{1}^{*},\ldots,Q_{n-1}^{*})$ gives

[TABLE]

We now have $n-1$ equations for the $n-1$ block entries of the block diagonal matrix $\widehat{\bm{\mathcal{Q}}}\widehat{\Gamma}\widehat{\bm{\mathcal{D}}}$. To ensure that a non-singular $\widehat{\bm{\mathcal{Q}}}\widehat{\Gamma}\widehat{\bm{\mathcal{D}}}$ satisfying Eq. 88 exists, there must be some consistency between $\langle\langle F_{k},F_{k}\rangle\rangle^{-1}-\langle\langle F_{k-1},F_{k-1}\rangle\rangle^{-1}$ and $C^{(k)}_{0}$. More precisely, it has to hold that, for $k=1,\ldots,n-1$,

[TABLE]

or alternatively101010Using the fact that $\text{Range}(R^{\ast})=\text{Range}(R^{\ast}R)$.

[TABLE]

If $\text{Range}\left(C^{(k)}_{0}\right)=\text{Range}\left(\langle\langle F_{k},F_{k}\rangle\rangle^{-1}-\langle\langle F_{k-1},F_{k-1}\rangle\rangle^{-1}\right)=\mathbb{C}^{s}$, then the entry $\bm{E}_{k}^{T}\widehat{\bm{\mathcal{D}}}\bm{E}_{k}$ is defined uniquely as

[TABLE]

In other cases satisfying (90), there is certain freedom in the components corresponding to the null space of $\sqrt{\langle\langle F_{k},F_{k}\rangle\rangle^{-1}-\langle\langle F_{k-1},F_{k-1}\rangle\rangle^{-1}}\,^{*}$.

In the sense of Section 3.2, Eq. 90 implies that the FOM residual norm must satisfy

[TABLE]

4.4 Final result

In the standard case, Eq. 90 reduces to the well-known condition that GMRES stagnates if and only if we obtain at least one zero Ritz value. In the block case, (partial) stagnation also appears if and only if $\bm{\mathcal{H}}^{(k)}$ is singular. But in addition, the rank deficiency of $\langle\langle F_{k},F_{k}\rangle\rangle^{-1}-\langle\langle F_{k-1},F_{k-1}\rangle\rangle^{-1}$, encoding the stagnation, and the rank deficiency of $C^{(k)}_{0}$, encoding the singularity of $\bm{\mathcal{H}}^{(k)}$, must have the same structure, i.e., the corresponding matrices must share the same range. If this is the case, we can prescribe the convergence of the block Arnoldi method and the block GMRES method at the same time, as summarized in the following theorem, which generalizes [7, Th. 3.6].

Theorem 4.1**.**

Let $\left\{M^{(k)}\right\}_{k=1}^{n}$ be any sequence of $\lambda$-matrices,

[TABLE]

$C_{0}^{(n)}$* nonsingular, and let $\{F_{k}\}_{k=0}^{n-1}$, $F_{k}\in\mathbb{S}^{+}$, be any sequence satisfying*

[TABLE]

Under the assumption that the two sequences satisfy the consistency condition

[TABLE]

the following two assertions are equivalent:

The residuals of blGMRES($\bm{\mathcal{A}}$,$\bm{B}$) satisfy

[TABLE]

and the principal submatrices of the Hessenberg matrix generated by blArnoldi($\bm{\mathcal{A}}$,$\bm{B}$) satisfy

[TABLE] 2. 2.

The matrix $\bm{\mathcal{A}}$ and the starting vector/right-hand side $\bm{B}$ are of the form

[TABLE]

where $\bm{\mathcal{V}}$ is a unitary matrix, $\bm{\mathcal{C}}$ is the block companion matrix corresponding to $M^{(n)}$,

[TABLE]

and $\bm{\mathcal{D}}$ is a block diagonal matrix with entries in $\mathbb{S}^{+}$ satisfying

[TABLE]

for some $Q_{k}\in\mathbb{S}$, $|Q_{k}|=I$.

Proof 4.2**.**

*The proof follows from the construction in Sections 4.1, 4.2 and 4.3. *

Theorem 4.1 shows that, similar to the standard case, we can prescribe the convergence behavior for blArnoldi independently of the residual convergence of blGMRES, as long as the (partially) stagnating iterations are reflected in the corresponding $\lambda$-matrix and vice versa.

Remark 4.3**.**

If we do not prescribe the blGMRES convergence, we only have a condition on $\bm{\mathcal{U}}$, and $\bm{\mathcal{D}}$ can be an arbitrary block diagonal matrix with entries in $\mathbb{S}^{+}$. Similarly, if the blArnoldi convergence is not prescribed, Eq. 101 only gives a condition on the first row of $(\bm{\mathcal{D}}\bm{\mathcal{U}})^{-1}$, which has to satisfy

[TABLE]

*cf. [8, Th. 1]. *

Remark 4.4**.**

The Ritz value companion transform $\bm{\mathcal{U}}$ makes the intermediate $\lambda$-matrices $M^{(k)}$ completely independent of the last $\lambda$-matrix $M^{(n)}$ and these are all independent of the matrix $\bm{\mathcal{D}}$ defining the subdiagonal entries of $\bm{\mathcal{H}}$. In the standard case, this implies that the (often used) residual measure of the Ritz value convergence

[TABLE]

*may provide little or no useful information about the convergence of the Ritz values to the eigenvalues of $A$; see [7, pp. 964–965]. With the block generalization of the Jordan form, see, e.g., [13], a similar conclusion is possible for the block version. Further analysis of this topic is however beyond the scope of this paper. *

In the next section, we discuss some of the fundamental differences between the results presented in the preceding papers on standard Arnoldi and GMRES, and the results regarding their block counterparts presented here.

4.5 The role of polynomials with matrix coefficients

In Section 4, the spectral properties of $\bm{\mathcal{A}}$ and the submatrices $\bm{\mathcal{H}}^{(k)}$ are prescribed through the $\lambda$-matrices $M^{(k)}(\lambda)$. The relation Eq. 51, equivalent to Eq. 56, together with the assumption that $\mathcal{K}_{n}(\bm{\mathcal{A}},\bm{B})$ is of full rank, means that $M(\bm{\mathcal{A}})$ is zero when evaluated on $n$ linearly independent block vectors, as it satisfies

[TABLE]

In the standard case, Eq. 106 implies $M(A)=0$, i.e., $M$ is the characteristic polynomial of $A$, and $M(A)\circ\bm{v}=0$, $\forall\bm{v}\in\mathbb{C}^{n}$. This can be also seen from

[TABLE]

In other words, requiring $M(A)\circ\bm{b}=0$ is equivalent to prescribing the eigenvalues of $A$.

In the block case, the situation is different. First, $M(\bm{\mathcal{A}})$ itself cannot be defined, because of the clash of dimensions. Further, despite the fact that any vector $\bm{V}\in\mathbb{S}^{n}$ can be written as a block linear combination of $\bm{V}_{1},\ldots,\bm{V}_{n}$, i.e., $\bm{V}=\sum_{i=1}^{n}\bm{V}_{i}D_{i}$, Eq. 106 does not imply $M(\bm{\mathcal{A}})\circ\bm{V}=0$, $\forall\bm{V}\in\mathbb{S}^{n}$. This is because

[TABLE]

The transition between Eq. 109 and Eq. 110 is only possible when $D_{i}$ and $C_{k}$ commute for all $i$ and $k$. This makes the standard case different from the block case.

From Eq. 57, it is clear that the eigenvalues of $\bm{\mathcal{A}}$ are defined by the eigenvalues of $\bm{\mathcal{C}}$ and coincide with the latent roots of $M$; see also Appendix A. For a given block companion matrix ${\bm{\mathcal{C}}\in\mathbb{S}^{n\times n}}$, the manifold of block companion matrices similar to $\bm{\mathcal{C}}$ has dimension $ns^{2}-ns$; see, e.g., [9]. Therefore the eigenvalues do not define $M$ uniquely. This can be seen also from the fact that the coefficients $C_{0},\ldots,C_{n-1}$ have $ns^{2}$ free parameters, while there only are $ns$ eigenvalues.

Instead of focusing on the eigenvalues of $\bm{\mathcal{A}}$, it is advantageous to remain in the block setting and look at the solvents of the $\lambda$-matrix instead; see Appendix A. Prescribing the (chain of) solvents $S_{1},\ldots,S_{n}$ of the $\lambda$-matrix Eq. 50, the coefficients of the $\lambda$-matrix $M$ are defined uniquely through Eq. 117. The eigenvalues of $\bm{\mathcal{A}}$ are the eigenvalues of the solvents, but changing the eigenvectors of the individual solvents will have impact on the eigenvectors of $\bm{\mathcal{A}}$. To achieve $C_{0}$ nonsingular, we only require that each of the solvents $S_{1},\ldots,S_{n}$ of $M$ is nonsingular. A similar result holds for the submatrices $\bm{\mathcal{H}}^{(k)}$, with the exception that singular solvents are allowed.

Concluding, in the standard case, if the eigenvalues of $\bm{\mathcal{A}}$ and the convergence behavior of both GMRES and Arnoldi are prescribed, then, in case of no stagnation of GMRES, all matrices and right-hand sides satisfying these conditions will be identical up to a unitary transform. Each step in which GMRES stagnates provides one extra free parameter, represented by the entry of the diagonal matrix $\bm{\mathcal{D}}$ in Eq. 101. In the block case, prescribing the eigenvalues of $\bm{\mathcal{A}}$ and the convergence behavior of both blGMRES and blArnoldi will still give us certain freedom in the choice of the eigenvectors of $\bm{\mathcal{A}}$, plus again some extra free parameters in the stagnating iterations. From this point of view, the eigenvalues of $\bm{\mathcal{A}}$ are even less indicative regarding the residual convergence behavior of blArnoldi and blGMRES than they are in the standard case.

5 Conclusions and open questions

The analysis of block Krylov subspace methods has always presented a challenge beyond those encountered with classical non-block methods. This is due to the interaction between the right-hand sides. We have demonstrated here that the ∗-algebra approach introduced in [11] enables us to cleanly obtain the same sort of results one sees for non-block Krylov subspace methods, which have previously been unavailable. These results fill a certain gap in the understanding of the convergence behavior of block Krylov methods for non-symmetric matrices.

We have thus obtained block versions of the fundamental results regarding the admissible and attainable convergence of standard Arnoldi and GMRES presented in a series of papers published over the last 25 years. Extending the framework introduced in [11], we were able to keep the formal notation as close as possible to the original results, with the block generalization of the Givens transformation allowing for extension of well-known textbook relations for the residuals of GMRES and FOM to the block case. This framework allows us to see the interdependence of the residual convergence behaviors of the individual systems. We explicitly formulated conditions on the admissible convergence behavior of the residuals of blGMRES using an appropriate block generalization of the norm.

Under the assumption that blGMRES does not converge prematurely, we were then able to completely characterize matrices and right-hand sides producing any prescribed admissible convergence behavior. Furthermore, spectral properties of the matrix can be enforced through the similarity to a block companion matrix. We showed that arbitrary convergence of blArnoldi for the eigenvalue problem is possible, and that arbitrary convergence of blGMRES and blArnoldi can be, under moderate assumptions, achieved simultaneously. Combining these results with the theory of block companion matrices, we showed that in a certain sense, increasing the number of right-hand sides reduces the predictive value of the eigenvalues of $\bm{\mathcal{A}}$.

It should be noted that, as Meurant pointed out in [19], these results all concern residual convergence behavior. It is observed that for a class of matrices (each with different spectral properties) constructed to exhibit a specific admissible residual convergence behavior, the actual error convergence behavior may very well exhibit dependence on spectral properties of the matrix. Thus, one interpretation of the results presented in this paper and the work of the last 25 years on this topic is that spectral properties of a non-Hermitian matrix may not a priori tell us much about the behavior of GMRES with respect to our chosen method of measuring convergence. Thus one may connect the results in this and related works to the notion that one should measure (residual) error using an appropriate norm. However, this is beyond the scope of the current paper.

Certain important aspects of blGMRES and blArnoldi convergence also are beyond the scope of this paper. Since the Arnoldi algorithm is based on full orthogonalization, the amount of data that needs to be stored, as well as the computational complexity, grows with each iteration. In practical computations, restarting the orthogonalization process is therefore often unavoidable. For this reason, analysis of the influence of the restarts on the admissible convergence represents another important research direction. Furthermore, the blArnoldi process breaks down when fewer than $s$ linearly independent Arnoldi vectors are generated in the $k$th step. The situation when no single system has converged but rather a linear combination of the columns of $\bm{X}$ lies in $\mathcal{K}_{k}(\bm{\mathcal{A}},\bm{B})$ is particularly unpleasant. In the considered framework, resolving this situation by reducing the block size is not directly possible, since by this, we change $\mathbb{S}$. To avoid change in the block size, the linearly dependent vectors can be replaced by some auxiliary (random) vectors. This direction will be further explored elsewhere.

Appendix A Polynomials with matrix coefficients

In this section, we recall some terminology and fundamental results from the theory of polynomials with matrix coefficients. Statements presented in this section are adopted from [5]. For a more comprehensive overview of the topic, we recommend [13].

Let $C_{k}\in\mathbb{C}^{s\times s}$, $k=0,\ldots,n-1$. We call

[TABLE]

a $\lambda$-matrix and

[TABLE]

a matrix polynomial. We call $\lambda\in\mathbb{C}$ a latent root of the $\lambda$-matrix $M$ if $M(\lambda)$ is singular, and $S\in\mathbb{C}^{s\times s}$ a right solvent of the matrix polynomial $M$ if $M(S)=0$. The latent roots and solvents are related as follows.

Theorem A.1**.**

If $S$ is a solvent of the matrix polynomial $M$, then the $\lambda$-matrix $M$ can be factorized as

[TABLE]

Corollary A.2**.**

*The $s$ eigenvalues of the solvent $S$ of the matrix polynomial $M$ are all latent roots of the $\lambda$-matrix $M$. *

A sequence of matrices $S_{1},\ldots,S_{n}$ forms a chain of solvents if

[TABLE]

It follows directly from Eq. 114 that

[TABLE]

The matrix $C$, which has the form

[TABLE]

is called the block companion matrix associated with the $\lambda$-matrix Eq. 111 or equivalently the matrix polynomial Eq. 112. Eigenvalues of the block companion matrix $C$ and the latent roots are related as follows.

Theorem A.3**.**

$\det(C-\lambda I)=(-1)^{ns}\det(I\lambda^{n}-C_{n-1}\lambda^{n-1}-\cdots-C_{0})$*. *

Corollary A.4**.**

The eigenvalues of the block companion matrix are the latent roots of the associated $\lambda$-matrix, therefore $M$ has exactly $ns$ latent roots.

We define the action of a matrix polynomial on a block vector as

[TABLE]

see [24, p. 108] or [11, p. 107].

Appendix B Auxiliary lemma

We use MATLAB notation in this lemma. In particular, $Z_{j}$ denotes the $j$th block entry and $\bm{Z}_{1:j}$ denotes the first $j$ block entries of the block vector $\bm{Z}$.

Lemma B.1**.**

Let $\bm{Z}\in\mathbb{S}^{k}$ be such that $\bm{\mathcal{I}}-\bm{Z}\bm{Z}^{*}$ is positive definite, and let $\bm{\mathcal{R}}_{Z}$ be the unique upper triangular Cholesky factor of $\bm{\mathcal{I}}-\bm{Z}\bm{Z}^{*}$. Then

[TABLE]

Proof B.2**.**

Since $\bm{\mathcal{R}}_{Z}^{-*}$ is lower triangular, it holds that

[TABLE]

Therefore, it suffices to investigate $\bm{E}_{k}^{T}\bm{\mathcal{R}}_{Z}^{-*}\bm{Z}$ and apply the lemma recursively.

We observe that

[TABLE]

where we used the push-through identity (Woodbury matrix identity) and the fact that $\bm{\mathcal{R}}_{Z}$ is upper triangular.

The block entry $\bm{E}_{k}^{T}\bm{\mathcal{R}}_{Z}\bm{E}_{k}$, i.e., the $(k,k)$ block entry of the Cholesky factor, can be obtained as

[TABLE]

Thus

[TABLE]

*which gives the desired statement. *

Note that $\bm{\mathcal{I}}-\bm{Z}\bm{Z}^{*}$ is positive definite if and only if $\bm{\mathcal{I}}-\bm{Z}^{*}\bm{Z}$ is positive definite.

Acknowledgments

The authors would like to thank Gérard Meurant for his valuable suggestions regarding presentation of the result of Sections 4.1, 4.2 and 4.3. We are also grateful to two anonymous referees and the handling editor for their comments and suggestions.