# A Flexible Generalized Conjugate Residual Method with Inner Orthogonalization and Deflated Restarting

Luiz Mariano Carvalho, Serge Gratton, Rafael Lago, Xavier Vasseur

► To cite this version:

Luiz Mariano Carvalho, Serge Gratton, Rafael Lago, Xavier Vasseur. A Flexible Generalized Conjugate Residual Method with Inner Orthogonalization and Deflated Restarting. SIAM Journal on Matrix Analysis and Applications, 2011, 32 (4), pp.1212 - 1235. 10.1137/100786253. hal-00650239v2

HAL Id: hal-00650239

https://hal.science/hal-00650239v2

Submitted on 3 Oct 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A FLEXIBLE GENERALIZED CONJUGATE RESIDUAL METHOD WITH INNER ORTHOGONALIZATION AND DEFLATED RESTARTING*

L. M. CARVALHO†, S. GRATTON‡, R. LAGO§, AND X. VASSEUR¶

**Abstract.** This work is concerned with the development and study of a minimum residual norm subspace method based on the generalized conjugate residual method with inner orthogonalization (GCRO) method that allows flexible preconditioning and deflated restarting for the solution of non-symmetric or non-Hermitian linear systems. First we recall the main features of flexible generalized minimum residual with deflated restarting (FGMRES-DR), a recently proposed algorithm of the same family but based on the GMRES method. Next we introduce the new inner-outer subspace method named FGCRO-DR. A theoretical comparison of both algorithms is then made in the case of flexible preconditioning. It is proved that FGCRO-DR and FGMRES-DR are algebraically equivalent if a collinearity condition is satisfied. While being nearly as expensive as FGMRES-DR in terms of computational operations per cycle, FGCRO-DR offers the additional advantage to be suitable for the solution of sequences of slowly changing linear systems (where both the matrix and right-hand side can change) through subspace recycling. Numerical experiments on the solution of multidimensional elliptic partial differential equations show the efficiency of FGCRO-DR when solving sequences of linear systems.

**Key words.** flexible or inner-outer Krylov subspace methods, variable preconditioning, deflation, iterative solver

**AMS subject classifications.** 65F10, 65N22, 15A06

**DOI.** 10.1137/100786253

**1. Introduction.** In recent years, several authors studied inner-outer Krylov subspace methods that allow variable preconditioning for the iterative solution of large sparse linear systems of equations. One of the first papers describing a subspace method with variable preconditioning is due to Axelsson and Vassilevski, who proposed the generalized conjugate gradient method [1]. See also [2, section 12.3] for additional references. Since then, numerous methods have been proposed to address the symmetric, nonsymmetric, or non-Hermitian cases; these include flexible conjugate gradient [20], flexible GMRES (FGMRES) [24], flexible QMR [31], and GMRESR [34], among others. This class of methods is required when preconditioning with a different (possibly nonlinear) operator at each iteration of a subspace method is considered. This notably occurs when adaptive preconditioners using information obtained from previous iterations [3, 14] are used or when inexact solutions of the preconditioning system using, e.g., adaptive cycling strategy in multigrid [19] or approximate interior solvers in domain decomposition methods [32, section 4.3] are considered. The latter

<footnote>†Department of Applied Mathematics, IME-UERJ, R. S. F. Xavier, 524, 629D, 20559-900, Rio de Janeiro, RJ, Brazil (luizmc@gmail.com). This author's work was supported by grant CNPq-473420/2007-4, coordinated by Professor Nelson Maculan.</footnote>

<footnote>‡INPT-IRIT, University of Toulouse and ENSEEIHT, 2 rue Camichel, BP 7122, F-31071 Toulouse Cedex 7, France, and CERFACS, 42, Avenue Gaspard Coriolis, F-31057 Toulouse Cedex 1, France (serge.gratton@enseeiht.fr).</footnote>

<footnote>§CERFACS, 42, Avenue Gaspard Coriolis, F-31057 Toulouse Cedex 1, France (lago@cerfacs.fr).</footnote>

<footnote>¶CERFACS and HiePACS project joint INRIA-CERFACS Laboratory, 42, Avenue Gaspard Coriolis, F-31057 Toulouse Cedex 1, France (vasseur@cerfacs.fr).</footnote>
situation is frequent when solving very large systems of linear equations resulting from
the discretization of partial differential equations in three dimensions. Thus flexible
Krylov subspace methods have gained a considerable interest in the recent years and
are the subject of both theoretical and numerical studies [27]. We refer the reader to
[29, section 10] for additional comments on flexible methods.

When nonvariable preconditioning is considered, the full GMRES method [23] is
often chosen for the solution of nonsymmetric or non-Hermitian linear systems be-
cause of its robustness and its minimum residual norm property [26]. Nevertheless to
control both the memory requirements and the computational cost of the orthogonal-
ization scheme, restarted GMRES is preferred; it corresponds to a scheme where the
maximal dimension of the approximation subspace is fixed. It means in practice that
the orthonormal basis built is thrown away at the end of the cycle. Since some infor-
mation is discarded at the restart, the convergence may stagnate and is expected to be
slower compared to full GMRES. Nevertheless to retain the convergence rate a num-
ber of techniques have been proposed; they fall in the class of augmented and deflated
methods; see, e.g., [4, 10, 11, 16, 25]. Deflated methods compute spectral informa-
tion at a restart and use this information to improve the convergence of the subspace
method. One of the most recent procedures based on a deflation approach is GMRES
with deflated restarting (GMRES-DR) [18]. This method reduces to restarted GM-
RES when no deflation is applied, but may provide a much faster convergence than
restarted GMRES for well-chosen deflation spaces as described in [18].

Quite recently a new minimum residual norm subspace method based on GMRES
allowing deflated restarting and variable preconditioning has been proposed in [15]. It
mainly attempted to combine the numerical features of GMRES-DR and the flexibil-
ity property of FGMRES. Numerical experiments in [15] have shown the efficiency of
FGMRES with deflated restarting (FGMRES-DR) on both academic and industrial
examples. In this paper we study a new minimum residual norm subspace method
based on the generalized conjugate method with inner orthogonalization (GCRO) [9]
allowing deflated restarting and variable preconditioning. It is named flexible gener-
alized conjugate residual method with inner orthogonalization and deflated restarting
(FGCRO-DR) and can be viewed as an extension of GCRO-DR [22] to the case of
variable preconditioning. A major advantage of FGCRO-DR over FGMRES-DR is
its ability to solve sequences of linear systems (where both the left- and right-hand
sides can change) through recycling [22]. In [22] Parks et al. mentioned that GCRO-
DR and GMRES-DR were algebraically equivalent, i.e., both methods produce the
same iterates in exact arithmetic when solving the same given linear system starting
from the same initial guess. When variable preconditioning is considered, it seems
therefore natural to ask whether FGCRO-DR and FGMRES-DR could also be alge-
braically equivalent. We address this question in this paper, and the main theoretical
developments that are proposed will help us to answer this question. The main
contributions of the paper are then twofold. First we prove that FGCRO-DR and
FGMRES-DR can be considered algebraically equivalent if a collinearity condition
between two certain vectors is satisfied at each cycle. When considering nonvari-
able preconditioning, these theoretical developments will also allow us to show the
algebraic equivalence between GCRO-DR and GMRES-DR that was stated without
proof in [22]. Second we carefully analyze the computational cost of FGCRO-DR and
show that the proposed method is nearly as expensive as FGMRES-DR in terms of
operations per cycle. Furthermore it is explained how to include subspace recycling
into FGCRO-DR, and numerical experiments are reported showing the efficiency of
FGCRO-DR.
This paper is organized as follows. In section 2 we introduce the general background of this study. We briefly recall the main properties of FGMRES-DR and then introduce the FGCRO-DR method both from a mathematical and an algorithmic point of view. Section 3 is mainly devoted to the analysis of both flexible methods. Therein we show that both methods can be algebraically equivalent in the flexible case if a certain collinearity condition is satisfied at each cycle. In section 4 we compare FGCRO-DR and FGMRES-DR in terms of computational operations per cycle and storage and discuss the solution of sequences of linear systems through subspace recycling. Finally we draw some conclusions and perspectives in section 5.

## 2. Flexible Krylov methods with restarting.

### 2.1. General setting.

*Notation.* Throughout this paper we denote by $||.||$ the Euclidean norm, by $I_k \in \mathbb{C}^{k \times k}$ the identity matrix of dimension $k$, and by $0_{i \times j} \in \mathbb{C}^{i \times j}$ the zero rectangular matrix with $i$ rows and $j$ columns. Given $N \in \mathbb{C}^{n \times m}$, $\Pi_{N^\perp} = I_n - N N^\dagger$ will represent the orthogonal projector onto range($N$)$^\perp$, where the superscript $\dagger$ refers to the Moore-Penrose inverse. Finally, given $Z_m = [z_1, \dots, z_m] \in \mathbb{C}^{n \times m}$, we will usually decompose $Z_m$ into two submatrices defined as $Z_k = [z_1, \dots, z_k] \in \mathbb{C}^{n \times k}$ and $Z_{m-k} = [z_{k+1}, \dots, z_m] \in \mathbb{C}^{n \times (m-k)}$.

*Setting.* We focus on minimum residual norm based subspace methods that allow flexible preconditioning for the iterative solution of

$$ (2.1) \qquad Ax = b, \quad A \in \mathbb{C}^{n \times n}, \quad x, b \in \mathbb{C}^n, $$

given an initial vector $x_0 \in \mathbb{C}^n$. In this paper $A$ is supposed to be nonsingular. Flexible methods refer to a class of methods where the preconditioner is allowed to vary at each iteration. We refer the reader to, e.g., [29] for a general introduction on Krylov subspace methods and to [29, section 10] and [26, section 9.4] for a review on flexible methods. The minimum residual norm GMRES method [23] has been extended by Saad [24] to allow variable preconditioning. The resulting algorithm known as FGMRES($m$) relies on the Arnoldi relation

$$ (2.2) \qquad AZ_m = V_{m+1}\bar{H}_m, $$

where $Z_m \in \mathbb{C}^{n \times m}$, $V_{m+1} \in \mathbb{C}^{n \times (m+1)}$ has orthonormal columns and $\bar{H}_m \in \mathbb{C}^{(m+1) \times m}$ is upper Hessenberg. We denote by $M_j$ the preconditioning operator at iteration $j$ and remark that $M_j$ may be a nonlinear preconditioning function. We will then denote by $M_j(v)$ the action of $M_j$ on a vector $v$. In (2.2), the columns of $V_{m+1}$ form an orthonormal basis of the subspace spanned by the vectors

$$ \{r_0, Az_1, \dots, Az_m\} \text{ with } r_0 = b - Ax_0, $$

whereas $Z_m = [z_1, \dots, z_m]$ and $V_m = [v_1, \dots, v_m]$ are related by

$$ Z_m = [M_1(v_1), \dots, M_m(v_m)] \quad \text{with} \quad v_1 = \frac{r_0}{\|r_0\|}. $$

At the end of the cycle an approximate solution $x_m \in \mathbb{C}^n$ is then found by minimizing the residual norm $||r_0 - AZ_m y||$ over the space $x_0 + \text{range}(Z_m)$. Thus we obtain that

$$ x_m = x_0 + Z_m y^*, $$

where $y^*$ is the solution of the following least-squares problem of size $(m+1) \times m$:

$$ y^* = \underset{y \in \mathbb{C}^m}{\operatorname{argmin}} ||r_0 - AZ_m y|| = \underset{y \in \mathbb{C}^m}{\operatorname{argmin}} ||r_0|| \quad e_1 - \bar{H}_m y^*, $$
where $e_1$ is the first canonical vector of $\mathbb{C}^{m+1}$. Flexible subspace methods with restarting are based on a procedure where the construction of the subspace is stopped after a certain number of steps (denoted by $m$ in this paper with $m < n$). The method is then restarted mainly to control both the memory requirements and the cost of the orthogonalization scheme. In FGMRES($m$) the restarting consists in taking as an initial guess the last iterate of the cycle ($x_m$).

The main focus of this paper is to present minimum residual norm subspace methods with deflated restarting that allow flexible preconditioning. Deflated restarting aims at determining an approximation subspace of dimension $m$ as a direct sum of two subspaces of smaller dimension, where one of these subspaces will contain relevant spectral information that will be kept for the next cycle. We refer the reader to, e.g., [25] and [29, section 9] for a review of augmented and deflated methods. Flexible methods with deflated restarting will notably satisfy the following flexible Arnoldi relation:

$$ (2.3) \qquad AZ_m = V_{m+1}\bar{H}_m \quad \text{with} \quad V_{m+1}^H V_{m+1} = I_{m+1}, $$

where $\bar{H}_m \in \mathbb{C}^{(m+1) \times m}$ is not necessarily of upper Hessenberg form. In this paper we call this relation a flexible Arnoldi-like relation due to its similarity to relation (2.2).

*Stagnation and breakdown.* We refer the reader to [27, section 6] for general comments and a detailed discussion on the possibility of both breakdown and stagnation in flexible inner-outer Krylov subspace methods. Although important, these issues are not addressed in this paper, and we assume that no breakdown occurs in the inner-outer subspace methods that will be proposed.

**2.2. Flexible GMRES with deflated restarting.** A number of techniques have been proposed to compute spectral information at a restart and use this information to improve the convergence rate of the Krylov subspace methods; see, e.g., [16, 17, 18, 25]. These techniques have been exclusively developed in the case of a fixed preconditioner. GMRES-DR is one of these methods. It focuses on removing (or deflating) the eigenvalues of smallest magnitude. A full subspace of dimension $k$, $k < m$ (and not only the approximate solution with minimum residual norm) is now retained at the restart, and the success of this approach has been demonstrated in many academic examples [16]. Approximations of eigenvalues of smallest magnitude are obtained by computing harmonic Ritz pairs of $A$ with respect to a certain subspace [18]. We present here a definition of a harmonic Ritz pair equivalent to the one introduced in [21, 30]; it will be of key importance when defining appropriate deflation strategies.

**DEFINITION 2.1 (harmonic Ritz pair).** Consider a subspace $U$ of $\mathbb{C}^n$. Given $B \in \mathbb{C}^{n \times n}$, $\theta \in \mathbb{C}$, and $y \in U$, $(\theta, y)$ is a harmonic Ritz pair of $B$ with respect to $U$ if and only if

$$ By - \theta y \perp BU $$

or, equivalently, for the canonical scalar product,

$$ \forall w \in \text{range}(BU) \quad w^H (By - \theta y) = 0. $$

We call $y$ a harmonic Ritz vector associated with the harmonic Ritz value $\theta$.

As in the case of fixed preconditioning, deflated restarting may also improve the convergence rate of flexible subspace methods. In [15] a deflated restarting procedure has been proposed for the FGMRES algorithm. The $i$th cycle of the resulting algorithm, called FGMRES-DR, is now briefly described, and we denote by
$r_0^{(i-1)} = b - Ax_0^{(i-1)}$, $V_{m+1}, \bar{H}_m$, and $Z_m$ the residual and matrices obtained at the end of the $(i-1)$th cycle.

Based on the Arnoldi-like relation (2.3), the deflation procedure proposed in [15, Proposition 1] relies on the use of $k$ harmonic Ritz vectors $Y_k = V_m P_k$ of $AZ_m V_m^H$ with respect to range($V_m$), where $Y_k \in \mathbb{C}^{n \times k}$ and $P_k \in \mathbb{C}^{m \times k}$. In Lemma 2.2 shown in [15, Lemma 3.1], we recall a useful relation satisfied by the harmonic Ritz vectors.

**LEMMA 2.2.** In FGMRES-DR, the harmonic Ritz vectors are given by $Y_k = V_m P_k$ with corresponding harmonic Ritz values $\lambda_k$. $P_k \in \mathbb{C}^{m \times k}$ satisfies the following relation:

$$ (2.4) \qquad AZ_m P_k = V_{m+1} \begin{bmatrix} P_k \\ 0_{1 \times k} \end{bmatrix}, c - \bar{H}_m y^* \begin{bmatrix} \operatorname{diag}(\lambda_1, \dots, \lambda_k) \\ \alpha_{1 \times k} \end{bmatrix}, $$

$$ (2.5) \qquad AZ_m P_k = [V_m P_k, r_0^{(i-1)}] \begin{bmatrix} \operatorname{diag}(\lambda_1, \dots, \lambda_k) \\ \alpha_{1 \times k} \end{bmatrix}, $$

where $r_0^{(i-1)} = V_{m+1}(c - \bar{H}_m y^*)$ and $\alpha_{1 \times k} = [\alpha_1, \dots, \alpha_k] \in \mathbb{C}^{1 \times k}$.

Next, the QR factorization of the $(m+1) \times (k+1)$ matrix appearing on the right-hand side of relation (2.4) is performed as

$$ (2.6) \qquad \begin{bmatrix} P_k \\ 0_{1 \times k} \end{bmatrix}, c - \bar{H}_m y^* = QR, $$

where $Q \in \mathbb{C}^{(m+1) \times (k+1)}$ has orthonormal columns and $R \in \mathbb{C}^{(k+1) \times (k+1)}$ is upper triangular, respectively. We write the matrix $Q$ obtained in relation (2.6) as

$$ (2.7) \qquad Q = \begin{bmatrix} Q_{m \times k} \\ 0_{1 \times k} \end{bmatrix}, \frac{\bar{\rho}}{\|\bar{\rho}\|}, $$

where $Q_{m \times k} \in \mathbb{C}^{m \times k}$ and $\bar{\rho} \in \mathbb{C}^{m+1}$ is defined as

$$ (2.8) \qquad \bar{\rho} = \left( I_{m+1} - \begin{bmatrix} Q_{m \times k} \\ 0_{1 \times k} \end{bmatrix} \begin{bmatrix} Q_{m \times k} \\ 0_{1 \times k} \end{bmatrix}^H \right) (c - \bar{H}_m y^*). $$

**PROPOSITION 1.** In FGMRES-DR, the flexible Arnoldi relation

$$ (2.9) \qquad A Z_k = V_{k+1} \bar{H}_k, $$

$$ (2.10) \qquad V_{k+1}^H V_{k+1} = I_{k+1}, $$

$$ (2.11) \qquad \operatorname{range} \left( \begin{bmatrix} Y_k, r_0^{(i-1)} \end{bmatrix} \right) = \operatorname{range} (V_{k+1}) $$

holds at the $i$th cycle with matrices $Z_k, V_k \in \mathbb{C}^{n \times k}$ and $\bar{H}_k \in \mathbb{C}^{(k+1) \times k}$ defined as

$$ (2.12) \qquad Z_k = Z_m Q_{m \times k}, $$

$$ (2.13) \qquad V_{k+1} = V_{m+1} Q, $$

$$ (2.14) \qquad \bar{H}_k = Q^H \bar{H}_m Q_{m \times k}, $$

where $V_{m+1}, Z_m$, and $\bar{H}_m$ refer to matrices obtained at the end of the $(i-1)$th cycle.

*Proof.* Relations (2.9), (2.10), (2.12), (2.13), and (2.14) have been shown in [15, Proposition 2]. From relations (2.13) and (2.6), respectively, we deduce

$$ (2.15) \qquad V_{k+1} R = V_{m+1} \begin{bmatrix} P_k \\ 0_{1 \times k} \end{bmatrix}, c - \bar{H}_m y^*, $$

$$ (2.15) \qquad V_{k+1} R = \begin{bmatrix} V_m P_k, r_0^{(i-1)} \end{bmatrix}, $$
which finally shows that range([Y_k, r_0^{(i-1)}]) = range(V_{k+1}) since R is supposed to be nonsingular. □

FGMRES-DR then carries out m-k Arnoldi steps with flexible preconditioning and starting vector $v_{k+1}$ while maintaining orthogonality to $V_k$, leading to

$$A [z_{k+1}, \dots, z_m] = [v_{k+1}, \dots, v_{m+1}] \bar{H}_{m-k} \quad \text{and} \quad V_{m+1}^H V_{m+1} = I_{m+1}.$$

We note that $\bar{H}_{m-k} \in \mathbb{C}^{(m-k+1) \times (m-k)}$ is upper Hessenberg. At the end of the ith cycle this gives the flexible Arnoldi-like relation

$$A [Z_k, Z_{m-k}] = [V_{m+1}] \left[ \begin{bmatrix} \bar{H}_k \\ 0_{m-k \times k} \end{bmatrix} \begin{bmatrix} B_{k \times m-k} \\ \bar{H}_{m-k} \end{bmatrix} \right],$$

where $V_{m+1} \in \mathbb{C}^{n \times (m+1)}$, $\bar{H}_m \in \mathbb{C}^{(m+1) \times m}$, and $B_{k \times m-k} \in \mathbb{C}^{k \times (m-k)}$ results from the orthogonalization of $[v_{k+2}, \dots, v_{m+1}]$ against $V_{k+1}$. We note that $\bar{H}_m$ is no longer upper Hessenberg due to the leading dense $(k+1) \times k$ submatrix $\bar{H}_k$. At the end of the ith cycle, an approximate solution $x_0^{(i)} \in \mathbb{C}^n$ is then found by minimizing the residual norm $\|b - A(x_0^{(i)} + Z_m y)\|$ over the space $x_0^{(i-1)} + \text{range}(Z_m)$, the corresponding residual being $r_0^{(i)} = b - Ax_0^{(i)}$, with $r_0^{(i)} \in \text{range}(V_{m+1})$. We refer the reader to [15] for the complete derivation of the method and numerical experiments showing the efficiency of FGMRES-DR in both academic and industrial examples.

**2.3. Flexible GCRO with deflated restarting.** GCRO-DR [22]—a combination of GMRES-DR and GCRO—is a Krylov subspace method that allows deflated restarting and subspace recycling simultaneously. This latter feature is particularly interesting when solving sequences of linear systems with possibly different left- or right-hand sides. As pointed out in [22], GCRO-DR is attractive because any subspace may be recycled. In this paper we restrict the presentation to the case of a single linear system as proposed in (2.1).

GCRO and GCRO-DR belong to the family of inner-outer methods [2, Chap. 12] where the outer iteration is based on GCR, a minimum residual norm method proposed by Eisenstat, Elman, and Schultz [13]. To this end GCR maintains a correction subspace spanned by $\text{range}(Z_m)$ and an approximation subspace spanned by $\text{range}(V_m)$, where $Z_m, V_m \in \mathbb{C}^{n \times m}$ satisfy

$$A Z_m = V_m, \\ V_m^H V_m = I_m.$$

The optimal solution of the minimization problem $\min \|b - Ax\|$ over the subspace $x_0 + \text{range}(Z_m)$ is then found as $x_m = x_0 + Z_m V_m^H r_0$. Consequently $r_m = b - A x_m$ satisfies

$$r_m = r_0 - V_m V_m^H r_0 = \Pi_{V_m^\perp} r_0, \quad r_m \perp \text{range}(V_m).$$

In [9] de Sturler proposed an improvement to GMRESR [34], an inner-outer method based on GCR in the outer part and GMRES in the inner part, respectively. He suggested that the inner iteration takes place in a subspace orthogonal to the outer Krylov subspace. In this inner iteration the projected residual equation

$$(I_n - V_m V_m^H) Az = r_m$$

is solved only approximately. If a minimum residual norm subspace method is used in the inner iteration to solve this projected residual linear system, the residuals over
both the inner and outer subspaces are minimized. This leads to the GCRO Krylov subspace method [9]. Numerical experiments [9] indicate that the resulting method may perform better than other inner-outer methods (without orthogonalization) in some cases.

The GCRO method with deflated restarting (named GCRO-DR) based on harmonic Ritz value information was proposed in [22]. An approximate invariant subspace is used for deflation following closely the GMRES-DR method. We refer the reader to [22] for a description of this method, algorithms, and implementation details. We present now a new variant of GCRO-DR that allows flexible preconditioning by explaining the different steps occurring during the ith cycle. Again we denote by $r_0^{(i-1)} = b - Ax_0^{(i-1)}$, $V_{m+1}, \bar{H}_m$, and $Z_m$ the residual and matrices obtained at the end of the (i-1)th cycle.

We suppose that a flexible Arnoldi-like relation of type (2.3) holds. As in section 2.2 an important point is to specify which harmonic Ritz information is selected. Given a certain matrix $W_m \in \mathbb{C}^{n \times m}$ to be specified later on, such as range($W_m$) = range($V_m$), the deflation procedure relies on the use of $k$ harmonic Ritz vectors $Y_k = W_m P_k$ of $AZ_m W_m^\dagger$ with respect to range($W_m$), where $Y_k \in \mathbb{C}^{n \times k}$ and $P_k \in \mathbb{C}^{m \times k}$. $W_m$ will notably satisfy a property detailed in Lemma 3.3 and we point out that the calculation of $W_m^\dagger$ is not needed in the practical implementation of the algorithm (see section 4.1.1). In Lemma 2.3 we detail a useful relation satisfied by the harmonic Ritz vectors.

**LEMMA 2.3.** In flexible GCRO with deflated restarting (FGCRO-DR), the harmonic Ritz vectors are given by $Y_k = W_m P_k$ with corresponding harmonic Ritz values $\theta_k$. The matrix $P_k = [p_1, \dots, p_k] \in \mathbb{C}^{m \times k}$ satisfies the following relation:

$$ (2.16) \qquad AZ_m P_k = \begin{bmatrix} W_m P_k, r_0^{(i-1)} \end{bmatrix} \begin{bmatrix} \operatorname{diag}(\theta_1, \dots, \theta_k) \\ \beta_{1 \times k} \end{bmatrix}, $$

where $r_0^{(i-1)} = V_{m+1}(c - \bar{H}_m y^*)$ and $\beta_{1 \times k} = [\beta_1, \dots, \beta_k] \in \mathbb{C}^{1 \times k}$.

*Proof.* According to Definition 2.1, the harmonic residual vectors $AZ_m W_m^\dagger W_m p_j - \theta_j W_m p_j$ and the residual vector $r_0^{(i-1)} = V_{m+1}(c - \bar{H}_m y^*)$ all belong to a subspace of dimension $m+1$ (spanned by the columns of $V_{m+1}$) and are orthogonal to the same subspace of dimension $m$ (spanned by the columns of $AZ_m$ subspace of range($V_{m+1}$)), so they must be collinear. Consequently there exist $k$ coefficients noted $\beta_j \in \mathbb{C}$ with $1 \le j \le k$ such that

$$ (2.17) \qquad \forall j \in \{1, \dots, k\} \quad AZ_m p_j - \theta_j W_m p_j = \beta_j r_0^{(i-1)}. $$

Setting $\beta_{1 \times k} = [\beta_1, \dots, \beta_k] \in \mathbb{C}^{1 \times k}$, the collinearity expression (2.17) can be written in matrix form as

$$ AZ_m P_k = \begin{bmatrix} W_m P_k, r_0^{(i-1)} \end{bmatrix} \begin{bmatrix} \operatorname{diag}(\theta_1, \dots, \theta_k) \\ \beta_{1 \times k} \end{bmatrix}. $$

Due to the flexible Arnoldi-like relation (2.3), relation (2.16) can be also expressed as

$$ (2.18) \qquad V_{m+1} \bar{H}_m P_k = \begin{bmatrix} W_m P_k, r_0^{(i-1)} \end{bmatrix} \begin{bmatrix} \operatorname{diag}(\theta_1, \dots, \theta_k) \\ \beta_{1 \times k} \end{bmatrix}. $$

If required, $\beta_{1 \times k}$ can be deduced from (2.18) by

$$ (2.19) \qquad (c - \bar{H}_m y^*)^H (\bar{H}_m P_k - V_{m+1}^H W_m P_k \operatorname{diag}(\theta_1, \dots, \theta_k)) \\ = (c - \bar{H}_m y^*)^H (c - \bar{H}_m y^*) \beta_{1 \times k}. \quad \square $$
Next, the QR factorization of the $(m+1) \times k$ matrix $\bar{H}_m P_k$ appearing in relation (2.18) is performed as $\bar{H}_m P_k = QR$ with $Q \in \mathbb{C}^{(m+1) \times k}$ and $R \in \mathbb{C}^{k \times k}$.

**PROPOSITION 2.** In FGCRO-DR, the relation $AZ_k = V_k$ with $V_k^H V_k = I_k$ holds at the ith cycle with matrices $Z_k, V_k \in \mathbb{C}^{n \times k}$ defined as

$$ 
\begin{aligned}
Z_k &= Z_m P_k R^{-1}, \\
V_k &= V_{m+1} Q,
\end{aligned}
$$

where $V_{m+1}$ and $Z_m$ refer to matrices obtained at the end of the $(i-1)$th cycle. In addition $V_k^H r_0^{(i-1)} = 0$ holds during the ith cycle.

*Proof.* By using information related to the QR factorization of $\bar{H}_m P_k$ and the flexible Arnoldi relation (2.3) exclusively, we obtain

$$ 
\begin{aligned}
A Z_k &= AZ_m P_k R^{-1} \\
&= V_{m+1} \bar{H}_m P_k R^{-1} \\
&= V_{m+1} Q \\
&= V_k.
\end{aligned}
$$

Since both $V_{m+1}$ and $Q$ have orthonormal columns, $V_k$ satisfies $V_k^H V_k = I_k$. Finally since $r_0^{(i-1)}$ is the optimum residual at the $(i-1)$th cycle, i.e., $(AZ_m)^H r_0^{(i-1)} = 0$, we obtain

$$ 
\begin{aligned}
P_k^H (AZ_m)^H r_0^{(i-1)} &= 0, \\
(V_{m+1} \bar{H}_m P_k)^H r_0^{(i-1)} &= 0, \\
R^H V_k^H r_0^{(i-1)} &= 0.
\end{aligned}
$$

This finally shows that $V_k^H r_0^{(i-1)} = 0$ since $R$ is supposed to be nonsingular. □

To complement the subspaces, the inner iteration is based on the approximate solution of

$$ (I_n - V_k V_k^H) A z = (I_n - V_k V_k^H) r_0^{(i-1)} = r_0^{(i-1)}, $$

where the last equality is due to Proposition 2. For that purpose FGCRO-DR then carries out $m-k$ steps of the Arnoldi method with flexible preconditioning, leading to

$$ 
\begin{aligned}
(I_n - V_k V_k^H) A [z_{k+1}, \dots, z_m] &= [v_{k+1}, \dots, v_{m+1}] \bar{H}_{m-k}, \\
(I_n - V_k V_k^H) A Z_{m-k} &= V_{m-k+1} \bar{H}_{m-k}
\end{aligned}
$$

with $v_{k+1} = r_0^{(i-1)}/\|r_0^{(i-1)}\|$. At the end of the cycle this gives the flexible Arnoldi-like relation

$$ 
\begin{aligned}
A [Z_k, Z_{m-k}] &= [V_k, V_{m-k+1}] \begin{bmatrix} I_k & V_k^H A Z_{m-k} \\ 0_{m-k+1 \times k} & \bar{H}_{m-k} \end{bmatrix}, \\
A Z_m &= V_{m+1} \bar{H}_m,
\end{aligned}
$$

where $Z_m \in \mathbb{C}^{n \times m}$, $V_{m+1} \in \mathbb{C}^{n \times (m+1)}$, and $\bar{H}_m \in \mathbb{C}^{(m+1) \times m}$. At the end of the ith cycle, an approximate solution $x_0^{(i)} \in \mathbb{C}^n$ is then found by minimizing the residual norm $\|b - A(x_0^{(i-1)} + Z_m y)\|$ over the space $x_0^{(i-1)} + \text{range}(Z_m)$, the corresponding residual being $r_0^{(i)} = b - Ax_0^{(i)}$, with $r_0^{(i)} \in \text{range}(V_{m+1})$.
**2.4. Algorithms.** Details of the FGCRO-DR method are given in Algorithm 1, where MATLAB-like notations are adopted (for instance, in step 7b, $Q(1:m, 1:k)$ denotes the submatrix made of the first $m$ rows and first $k$ columns of matrix $Q$ noted $Q_{m \times k}$ in (2.7)). For the sake of completeness the FGMRES-DR algorithm has also been described with notation chosen as closely as possible to FGCRO-DR to make code comparison easier. Concerning Algorithm 1 we make the following comments:

* As will be discussed later, the computation of $W_m^\dagger$ in step 5a is not required thanks to the definition of the harmonic Ritz pair (see Definition 2.1).

* As pointed out by Morgan [18] and Parks et al. [22] we might have to adjust $k$ during the algorithm to include both the real and imaginary parts of complex eigenvectors.

* In steps 10a and 10b $M_j^{(i)}$ denotes the possibly nonlinear preconditioning operator at iteration $j$ during the $i$th cycle.

**Algorithm 1.** FGCRO-DR($m, k$) and FGMRES-DR($m, k$).

1: Choose $m, k, \text{tol}$, and $x_0$

2: $r_0 = b - Ax_0$, $\beta = \|r_0\|$, $v_1 = r_0/\beta$, $c = \beta e_1$, $i \leftarrow 0$

3: Apply FGMRES(m) to obtain $\bar{H}_m$, $Z_m$, $V_{m+1}$ such that $AZ_m = V_{m+1}\bar{H}_m$, $y^* = \arg \min_{y \in \mathbb{C}^m} \|c - \bar{H}_m y\|$, $x_0^{(0)} = x_0 + Z_m y^*$, $r_0^{(0)} = b - Ax_0^{(0)} = V_{m+1}(c - \bar{H}_m y^*)$, $W_m = V_m$

4: **while** $\|r_0^{(i)}\| > \|b\| \times \text{tol}$ **do** $i \leftarrow i+1$

**FGCRO-DR**

5a: Compute $k$ harmonic Ritz vectors of $AZ_m W_m^\dagger$ with respect to range($W_m$) and store them in $Y_k$. Define $P_k$ such that $Y_k = W_m P_k$

6a: $Q R = \bar{H}_m P_k$

7a: $W_k = W_m P_k R^{-1}$

8a: $V_k = V_{m+1} Q$

9a: $Z_k = Z_m P_k R^{-1}$

10a: Apply $m-k$ flexible preconditioned Arnoldi steps with $(I_n - V_k V_k^H) A$ and $v_{k+1} = r_0^{(i-1)}/\|r_0^{(i-1)}\|$ such that $(I_n - V_k V_k^H) A [z_{k+1}, \dots, z_m] = [v_{k+1}, \dots, v_{m+1}] \bar{H}_{m-k}$ with $z_j = M_j^{(i)}(v_j)$

11a: Set $\bar{H}_m = \begin{bmatrix} I_k & V_k^H A Z_{m-k} \\ 0_{m-k+1 \times k} & \bar{H}_{m-k} \end{bmatrix}$ yielding $A[z_1, \dots, z_m] = [v_1, \dots, v_{m+1}] \bar{H}_m$ and define $W_m = [W_k \ V_m(1:n,k+1:m)]$

12: $y^* = \arg \min_{y \in \mathbb{C}^m} \|c - \bar{H}_m y\|$ with $c = V_{m+1}^H r_0^{(i-1)}$

13: $x_0^{(i)} = x_0^{(i-1)} + Z_m y^*$

14: $r_0^{(i)} = b - Ax_0^{(i)}$

15: **end while**

**FGMRES-DR**

5b: Compute $k$ harmonic Ritz vectors of $AZ_m V_m^H$ with respect to range($V_m$) and store them in $Y_k$. Define $P_k$ such that $Y_k = V_m P_k$

6b: $QR = \begin{bmatrix} P_k \\ 0_{1 \times k} \end{bmatrix} c - \bar{H}_m y^*$

7b: $\bar{H}_k = Q^H \bar{H}_m Q(1:m, 1:k)$

8b: $V_{k+1} = V_{m+1} Q$

9b: $Z_k = Z_m Q(1:m, 1:k)$

10b: Apply $m-k$ flexible preconditioned Arnoldi steps with $A$ and $v_{k+1}$ while maintaining orthogonality to $V_k$ such that $A[z_{k+1}, \dots, z_m] = [v_{k+1}, \dots, v_{m+1}] \bar{H}_{m-k}$

with $z_j = M_j^{(i)}(v_j)$ and $V_{m+1}^H V_{m+1} = I_{m+1}$

11b: Set $\bar{H}_m = \begin{bmatrix} \bar{H}_k \\ 0_{m-k \times k} \end{bmatrix} \begin{bmatrix} B_{k \times m-k} \\ \bar{H}_{m-k} \end{bmatrix}$ yielding $A[z_1, \dots, z_m] = [v_1, \dots, v_{m+1}] \bar{H}_m$
**3. Analysis of FGMRES-DR and FGCRO-DR.** We compare now the flexible variants of GMRES-DR and GCRO-DR introduced in sections 2.2 and 2.3, respectively. In the following we use ~ to denote quantities related to the FGMRES-DR algorithm, e.g., $\tilde{Y}_k$ denotes the set of harmonic Ritz vectors computed in the FGMRES-DR algorithm. When analyzing both algorithms we will suppose that identical preconditioning operators are used in steps 10a and 10b, respectively, i.e.,

$$ (3.1) \qquad \forall i, \forall j \in \{k+1, \dots, m\}, \quad \mathcal{M}_j^{(i)}(\cdot) = \tilde{\mathcal{M}}_j^{(i)}(\cdot). $$

**3.1. Equivalent preconditioning matrix.**

**DEFINITION 3.1** (equivalent preconditioning matrix). Suppose that $V_p = [v_1, \dots, v_p] \in \mathbb{C}^{n \times p}$ and $Z_p = [\mathcal{M}_1(v_1), \dots, \mathcal{M}_p(v_p)] \in \mathbb{C}^{n \times p}$ obtained during a cycle of a flexible method with (standard or deflated) restarting (with $1 \le p \le m < n$) are both of full rank, i.e., $\text{rank} V_p = \text{rank} Z_p = p$. We will then denote by $M_{V_p} \in \mathbb{C}^{n \times n}$ a nonsingular equivalent preconditioning matrix defined as

$$ (3.2) \qquad Z_p \stackrel{\text{def}}{=} M_{V_p} V_p. $$

Such a matrix represents the action of the nonlinear operators $\mathcal{M}_j$ on the set of vectors $v_j$ (with $j = 1, \dots, p$). It can be chosen, e.g., as $M_{V_p} = [Z_p Z_p^{-1}] [V_p V_p^{-1}]^{-1}$, where $Z_p$ (respectively, $V_p$) denotes an orthogonal complement of $Z_p$ (respectively, $V_p$) in $\mathbb{C}^n$.

**3.2. Relations between $Z_m$ and $W_m$ and $\tilde{Z}_m$ and $\tilde{V}_m$.** We denote by $M_{W_m}^{(0)}$ and $\tilde{M}_{\tilde{V}_m}^{(0)}$ the equivalent preconditioning matrices used in the initialization phase of both algorithms (step 3 in Algorithm 1). With this notation we remark that the following relations hold:

$$ (3.3) \qquad Z_m = M_{W_m}^{(0)} W_m = \tilde{Z}_m = \tilde{M}_{\tilde{V}_m}^{(0)} \tilde{V}_m. $$

We first analyze the relation between $\tilde{Z}_m$ and $\tilde{V}_m$.

**LEMMA 3.2.** At the end of the ith cycle of the FGMRES-DR method, $\tilde{Z}_m$ and $\tilde{V}_m$ satisfy

$$ (3.4) \qquad \tilde{Z}_m = \tilde{M}_{\tilde{V}_m}^{(i)} \tilde{V}_m = [\tilde{M}_{\tilde{V}_m}^{(i-1)} \tilde{V}_k, \tilde{M}_{\tilde{V}_{m-k}}^{(i)} \tilde{V}_{m-k}]. $$

*Proof.* The initialization phase leads to the relation $\tilde{Z}_m = \tilde{M}_{\tilde{V}_m}^{(0)} \tilde{V}_m$. We suppose that at the end of the $(i-1)$th cycle the following relation holds: $\tilde{Z}_m = \tilde{M}_{\tilde{V}_m}^{(i-1)} \tilde{V}_m$. At step 9b of the ith cycle, $\tilde{Z}_k$ is defined as

$$ \tilde{Z}_k = \tilde{Z}_m \tilde{Q}_{m \times k} = \tilde{M}_{\tilde{V}_m}^{(i-1)} \tilde{V}_m \tilde{Q}_{m \times k} = \tilde{M}_{\tilde{V}_m}^{(i-1)} \tilde{V}_k. $$

The proof is then complete since $\tilde{Z}_{m-k} = [\tilde{M}_{k+1}^{(i)}(\tilde{v}_{k+1}), \dots, \tilde{M}_m^{(i)}(\tilde{v}_m)] = \tilde{M}_{\tilde{V}_{m-k}}^{(i)} \tilde{V}_{m-k}$ at the end of step 10b. □

The next lemma details a relation between $Z_m$ and $W_m$ that is satisfied in the FGCRO-DR method.

**LEMMA 3.3.** At the end of the ith cycle of the FGCRO-DR method, $Z_m$ and $W_m$ satisfy

$$ (3.5) \qquad Z_m = M_{W_m}^{(i)} W_m = [M_{W_m}^{(i-1)} W_k, M_{W_{m-k}}^{(i)} W_{m-k}]. $$
Proof. The initialization phase leads to the relation $Z_m = M_{W_m}^{(0)} W_m$. We suppose that at the end of the (i-1)th cycle the following relation holds: $Z_m = M_{W_m}^{(i-1)} W_m$. At step 9a of the ith cycle, $Z_k$ is defined as

$$ 
\begin{aligned}
Z_k &= Z_m P_k R^{-1} \\
&= M_{W_m}^{(i-1)} W_m P_k R^{-1} \\
&= M_{W_m}^{(i-1)} W_k.
\end{aligned}
$$

The proof is then complete since $Z_{m-k} = [M_{k+1}^{(i)}(w_{k+1}), \dots, M_m^{(i)}(w_m)] = M_{W_{m-k}}^{(i)} W_{m-k}$ at the end of step 11a. □

Lemmas 3.2 and 3.3 show that $\tilde{Z}_m, \tilde{V}_m, Z_m$, and $W_m$ satisfy relations that will play a central role in section 3.3. We investigate next the relation between $Z_m$ and $V_m$.

LEMMA 3.4. At the end of the ith cycle of the FGCRO-DR method, $Z_m$ and $V_m$ satisfy

$$ (3.6) \qquad [AZ_k, Z_{m-k}] = [V_k, M_{V_{m-k}}^{(i)} V_{m-k}]. $$

Proof. We use the relation $AZ_k = V_k$ satisfied in the FGCRO-DR method shown in Proposition 2. The proof is then complete since $Z_{m-k} = [M_{k+1}^{(i)}(v_{k+1}), \dots, M_m^{(i)}(v_m)] = M_{V_{m-k}}^{(i)} V_{m-k}$ at the end of step 11a. □

We conclude this section by presenting a technical lemma related to the FGMRES-DR method.

LEMMA 3.5. During the ith cycle of the FGMRES-DR method, $\tilde{v}_{k+1}$ satisfies the relation

$$ (3.7) \qquad \tilde{v}_{k+1} = \tilde{v}_{k+1} / \|\tilde{v}_{k+1}\| \quad \text{with} \quad \tilde{v}_{k+1} = \Pi_{[\tilde{Y}_k]^{\perp}} \tilde{r}_0^{(i-1)}, $$

where $\tilde{r}_0^{(i-1)} = b - A\tilde{x}_0^{(i-1)}$ denotes the residual obtained at the end of the (i-1)th cycle.

Proof. Using Proposition 1 and relation (2.8) we obtain

$$ 
\begin{aligned}
\tilde{v}_{k+1} &= \tilde{V}_{m+1}\bar{\rho} = \tilde{r}_0^{(i-1)} - \tilde{V}_{m+1} \begin{bmatrix} \tilde{Q}_{m \times k} \\ 0_{1 \times k} \end{bmatrix} \begin{bmatrix} \tilde{Q}_{m \times k} \\ 0_{1 \times k} \end{bmatrix}^H \tilde{V}_{m+1}^H \tilde{r}_0^{(i-1)}, \\
\tilde{v}_{k+1} &= \tilde{V}_{m+1}\bar{\rho} = \tilde{r}_0^{(i-1)} - \tilde{V}_m \tilde{Q}_{m \times k} (\tilde{V}_m \tilde{Q}_{m \times k})^H \tilde{r}_0^{(i-1)}.
\end{aligned}
$$

Since $\tilde{V}_m \tilde{Q}_{m \times k}$ has orthonormal columns, this last expression now becomes

$$ \bar{v}_{k+1} = \Pi_{[\tilde{V}_m \tilde{Q}_{m \times k}]^{\perp}} \tilde{r}_0^{(i-1)}. $$

Because $\tilde{Q}_{m \times k}$ is the orthogonal factor of the QR decomposition of $\tilde{P}_k$, we obtain

$$ \mathrm{range}(\tilde{V}_m \tilde{P}_k) = \mathrm{range}(\tilde{V}_m \tilde{Q}_{m \times k}). $$

Since from Lemma 2.3 $\tilde{Y}_k = \tilde{V}_m \tilde{P}_k$, the proof is then complete. □
### 3.3. Analysis of the FGMRES-DR and FGCRO-DR methods. Lemma

3.3 has already described an important property satisfied by $W_m$ in the FGCRO-DR method proposed in Algorithm 1. We will analyze further the relation between the FGMRES-DR and FGCRO-DR methods. The next theorem states that the two flexible methods generate the same iterates in exact arithmetic under some conditions involving notably two vectors.

**THEOREM 3.6.** We denote by $r_0^{(i)} = b - Ax_0^{(i)}$ the residual obtained at the end of the ith cycle of the FGCRO-DR method (see step 14 of Algorithm 1). We suppose that Definition 3.1 holds and that the same equivalent preconditioning matrix is obtained at the end of the ith cycle of both the FGCRO-DR and FGMRES-DR algorithms, i.e., $M_{W_m}^{(i)} = \tilde{M}_{\tilde{V}_m}^{(i)}$. Under this assumption the harmonic Ritz vectors $\tilde{Y}_k$ and $Y_k$ can be chosen equal during the (i+1)th cycle. If in addition there exists a real-valued positive coefficient $\eta_{i+1}$ such that

$$ (3.8) \quad \Pi_{[Y_k, r_0^{(i)}/\|r_0^{(i)}\|]^{\perp}} A \mathcal{M}_{k+1}^{(i+1)} \left( \Pi_{Y_k^{\perp}} r_0^{(i)}/\|\Pi_{Y_k^{\perp}} r_0^{(i)}\| \right) \\ = \eta_{i+1} \Pi_{[Y_k, r_0^{(i)}/\|r_0^{(i)}\|]^{\perp}} A \mathcal{M}_{k+1}^{(i+1)} \left( r_0^{(i)}/\|r_0^{(i)}\| \right) $$

in the FGCRO-DR algorithm, then both algorithms generate the same iterates in exact arithmetic and

$$ (3.9) \qquad \text{range}(V_{m+1}) = \text{range}(\tilde{V}_{m+1}), $$

$$ (3.10) \qquad \text{range}(Z_m) = \text{range}(\tilde{Z}_m), $$

with

$$ (3.11) \qquad V_{m+1} = [\tilde{V}_{k+1}\hat{Q}, v_{k+2}, \dots, v_{m+1}], \qquad \tilde{V}_{m+1} = [\tilde{V}_{k+1}, v_{k+2}, \dots, v_{m+1}], $$

$$ (3.12) \qquad Z_m = [\tilde{Z}_{k+1}\hat{X}, z_{k+2}, \dots, z_m], \qquad \tilde{Z}_m = [\tilde{Z}_{k+1}, z_{k+2}, \dots, z_m], $$

where $\hat{Q} \in \mathbb{C}^{(k+1) \times (k+1)}$ is a unitary matrix and $\hat{X} \in \mathbb{C}^{(k+1) \times (k+1)}$ is a nonsingular triangular matrix.

*Proof.* The whole proof is performed in three parts assuming that we analyze the (i + 1)th cycle of each algorithm. Suppose that at the beginning of the (i + 1)th cycle (step 4) there exist a unitary matrix $\hat{Q} \in \mathbb{C}^{(k+1) \times (k+1)}$ and a nonsingular matrix $\hat{X} \in \mathbb{C}^{(k+1) \times (k+1)}$ such that the following relations hold:

$$ (3.13) \qquad V_{k+1} = \tilde{V}_{k+1}\hat{Q}, $$

$$ (3.14) \qquad Z_{k+1} = \tilde{Z}_{k+1}\hat{X}, $$

$$ (3.15) \qquad [v_{k+2}, \dots, v_{m+1}] = [\tilde{v}_{k+2}, \dots, \tilde{v}_{m+1}], $$

$$ (3.16) \qquad [z_{k+2}, \dots, z_m] = [\tilde{z}_{k+2}, \dots, \tilde{z}_m]. $$

We will then prove the existence of a unitary matrix $\hat{Q}' \in \mathbb{C}^{(k+1) \times (k+1)}$ and of a nonsingular matrix $\hat{X}' \in \mathbb{C}^{(k+1) \times (k+1)}$ such that at the end of the (i + 1)th cycle

$$ (3.17) \qquad V_{k+1} = \tilde{V}_{k+1}\hat{Q}', $$

$$ (3.18) \qquad Z_{k+1} = \tilde{Z}_{k+1}\hat{X}', $$

$$ (3.19) \qquad [v_{k+2}, \dots, v_{m+1}] = [\tilde{v}_{k+2}, \dots, \tilde{v}_{m+1}], $$

$$ (3.20) \qquad [z_{k+2}, \dots, z_m] = [\tilde{z}_{k+2}, \dots, \tilde{z}_m]. $$
Regarding FGCRO-DR we assume that at the beginning of the (i + 1)th cycle (step 4)

$$
(3.21) \qquad \text{range}(W_m) = \text{range}(V_m).
$$

We will also prove that relation (3.21) holds at the end of the (i+1)th cycle. Note that relations (3.9), (3.10), and (3.21) are obviously satisfied before the first cycle, because steps 1 to 3 are identical in both algorithms, yielding $V_{m+1} = \tilde{V}_{m+1}$, $Z_m = \tilde{Z}_m$, and $W_m = V_m$. Finally a consequence of (3.13), (3.15), (3.14), and (3.16) is that the residuals of the linear system $Ax = b$ in both algorithms are equal at the beginning of the (i + 1)th cycle, i.e., $r_0^{(i)} = \tilde{r}_0^{(i)}$. We will denote by $r_0$ this residual for ease of notation.

Part I: Steps 5a and 5b. In this part, we prove that we can choose $\tilde{Y}_k = Y_k$ with $Y_k = W_m P_k = \tilde{V}_m \tilde{P}_k$.

FGCRO-DR. Let $y_j = W_m p_j$ be the $j$th column of $Y_k$. Since $y_j$ is a harmonic Ritz vector of $AZ_m W_m^\dagger$ with respect to $\text{range}(W_m)$, the following relation holds (see Definition 2.1):

$$
(3.22) \qquad Z_m^H A^H (AZ_m p_j - \theta_j W_m p_j) = 0.
$$

Due to (3.14) and (3.16) there exists a nonsingular matrix $X \in \mathbb{C}^{m \times m}$ that relates $Z_m$ and $\tilde{Z}_m$:

$$
(3.23) \qquad Z_m = \tilde{Z}_m X.
$$

Using (3.23), the harmonic Ritz relation (3.22) now becomes

$$
X^H \tilde{Z}_m^H A^H (A \tilde{Z}_m X p_j - \theta_j W_m p_j) = 0.
$$

From Lemma 3.3 and relation (3.23) we deduce

$$
\begin{align*}
& X^H \tilde{Z}_m^H A^H (A \tilde{Z}_m X p_j - \theta_j M_{W_m}^{(i)^{-1}} Z_m p_j) = 0, \\
& X^H \tilde{Z}_m^H A^H (A \tilde{Z}_m X p_j - \theta_j \tilde{M}_{\tilde{V}_m}^{(i)^{-1}} \tilde{Z}_m X p_j) = 0,
\end{align*}
$$

where we have used explicitly the assumption on the equivalent preconditioning matrix
obtained at the end of the *i*th cycle, i.e., $M_{W_m}^{(i)} = \tilde{M}_{\tilde{V}_m}^{(i)}$. Next, the application of
Lemma 3.2 leads to

$$
(3.24) \qquad X^H \tilde{Z}_m^H A^H (A \tilde{Z}_m \tilde{V}_m^H \tilde{V}_m X p_j - \theta_j \tilde{V}_m X p_j) = 0.
$$

Since $X$ is nonsingular the last equality proves that $\tilde{V}_m X p_j$ is a harmonic Ritz vector
of $A \tilde{Z}_m \tilde{V}_m^H$ with respect to range($\tilde{V}_m$) associated to the Ritz value $\theta_j$. From relations
(3.22) and (3.24) we deduce that the harmonic Ritz vectors can be chosen to be equal
and correspond to the same harmonic Ritz values. In this case they notably satisfy
the following equality:

$$
(3.25) \qquad \forall j \in \{1, \dots, k\}, \quad \tilde{V}_m X p_j = W_m p_j, \quad \text{i.e.,} \quad \tilde{p}_j = X p_j.
$$

We will then denote by $Y = \tilde{Y}_k = Y_k$ the $k$ harmonic Ritz vectors computed in either FGCRO-DR or FGMRES-DR. We assume that the harmonic Ritz values $\theta_j$ ($1 \le j \le k$) are nonzero.
Part IIa: Steps 6a–10a, 6b–10b. We show that at the end of steps 10a and 10b the following relations hold: range($V_{k+1}$) = range($\tilde{V}_{k+1}$) = range([Y, $r_0^{(i)}/\|r_0^{(i)}\|$]). This result will help us to prove the existence of the matrix $\hat{Q}'$ introduced in relation (3.17).

FGCRO-DR. Since $AZ_mP_k = V_kR$ (Proposition 2), we deduce from Lemma 2.3

$$ (3.26) \quad \left[ V_k, r_0^{(i)} / \|r_0^{(i)}\| \right] = \left[ Y, r_0^{(i)} / \|r_0^{(i)}\| \right] \left[ \begin{array}{cc} \operatorname{diag}(\theta_1, \dots, \theta_k) R^{-1} & 0_{k \times 1} \\ \|r_0^{(i)}\| \beta_{1 \times k} R^{-1} & 1 \end{array} \right]. $$

This relation leads to the following result:

$$ (3.27) \quad \operatorname{range}(V_{k+1}) = \operatorname{range} \left( \left[ Y, r_0^{(i)} / \|r_0^{(i)}\| \right] \right). $$

Similarly $W_{k+1} = [W_k, \frac{r_0^{(i)}}{\|r_0^{(i)}\|}]$, using $Y = W_mP_k$, can be written as

$$ \begin{align} [W_k, r_0^{(i)} / \|r_0^{(i)}\|] &= \left[ W_m P_k R^{-1}, \frac{r_0^{(i)}}{\|r_0^{(i)}\|} \right] \nonumber \\ &= [Y R^{-1}, r_0^{(i)} / \|r_0^{(i)}\|] \nonumber \\ (3.28) \qquad &= [Y, r_0^{(i)} / \|r_0^{(i)}\|] \left[ \begin{array}{cc} R^{-1} & 0_{k \times 1} \\ 0_{1 \times k} & 1 \end{array} \right]. \end{align} $$

From relations (3.28) and (3.27) we deduce that

$$ (3.29) \quad \operatorname{range}(W_{k+1}) = \operatorname{range}(V_{k+1}). $$

This last result also proves that $\operatorname{range}(W_m) = \operatorname{range}(V_m)$ at the end of the cycle.

FGMRES-DR. In Proposition 1 we have shown that

$$ (3.30) \quad \operatorname{range}(\tilde{V}_{k+1}) = \operatorname{range} \left( \left[ Y, r_0^{(i)} / \|r_0^{(i)}\| \right] \right). $$

Since both $V_{k+1}$ and $\tilde{V}_{k+1}$ have orthonormal columns, we deduce from (3.27) and (3.30) that there exists a unitary matrix $\hat{Q}'$ such that

$$ (3.31) \quad V_{k+1} = \tilde{V}_{k+1} \hat{Q}', $$

which proves the relation proposed in (3.17).

Part IIb: Steps 6a–10a, 6b–10b. We show that at the end of steps 10a and 10b the following relation holds: range($Z_{k+1}$) = range($\tilde{Z}_{k+1}$). This result will help us to prove the existence of the matrix $\hat{X}'$ introduced in relation (3.18).

FGCRO-DR. Concerning $Z_{k+1} = [Z_k, z_{k+1}]$, there exists a nonsingular matrix $M^{(i+1)}_{[W_k, r_0^{(i)}/\|r_0^{(i)}\|]} \in \mathbb{C}^{n \times n}$ (see Definition 3.1) such that

$$ Z_{k+1} = M^{(i+1)}_{[W_k, r_0^{(i)}/\|r_0^{(i)}\|]} \left[ W_k, r_0^{(i)}/\|r_0^{(i)}\| \right]. $$

If $T \in \mathbb{C}^{(k+1) \times (k+1)}$ denotes the triangular matrix

$$ T = \begin{bmatrix} R & 0_{k \times 1} \\ 0_{1 \times k} & 1 \end{bmatrix} $$

due to relation (3.28), then $Z_{k+1}T$ can be written as

$$ (3.32) \quad Z_{k+1}T = M^{(i+1)}_{[W_k, r_0^{(i)}/\|r_0^{(i)}\|]} \left[ Y, r_0^{(i)}/\|r_0^{(i)}\| \right]. $$
FGMRES-DR. Similarly, from Lemma 3.2, $\tilde{Z}_{k+1}$ can be expressed as

$$ \tilde{Z}_{k+1} = \tilde{M}_{\tilde{V}_{k+1}}^{(i+1)} \tilde{V}_{k+1}, $$

where $\tilde{M}_{\tilde{V}_{k+1}}^{(i+1)} \in \mathbb{C}^{n \times n}$ is nonsingular (see Definition 3.1). If $\tilde{T} \in \mathbb{C}^{(k+1) \times (k+1)}$ denotes the triangular matrix

$$ \tilde{T} = \tilde{R} \begin{bmatrix} I_k & 0_{k \times 1} \\ 0_{1 \times k} & 1 / \|r_0^{(i)}\| \end{bmatrix}, $$

$\tilde{Z}_{k+1} \tilde{T}$ can be expressed as

$$ (3.33) \qquad \tilde{Z}_{k+1} \tilde{T} = \tilde{M}_{\tilde{V}_{k+1}}^{(i+1)} [Y, r_0^{(i)} / \|r_0^{(i)}\|] $$

thanks to relation (2.15). Relations (3.32) and (3.33) characterize $Z_{k+1}T$ and $\tilde{Z}_{k+1} \tilde{T}$ with respect to $[Y, r_0^{(i)} / \|r_0^{(i)}\|]$. We can further improve this result by showing the following equality:

$$ (3.34) \qquad M_{[W_k, r_0^{(i)} / \|r_0^{(i)}\|]}^{(i+1)} [Y, r_0^{(i)} / \|r_0^{(i)}\|] = \tilde{M}_{\tilde{V}_{k+1}}^{(i+1)} [Y, r_0^{(i)} / \|r_0^{(i)}\|]. $$

Lemma 3.3 and Lemma 3.2, respectively, give us two useful relations for $M_{[W_k, r_0^{(i)} / \|r_0^{(i)}\|]}^{(i+1)}$ [Y, $r_0^{(i)} / \|r_0^{(i)}\|$] and $\tilde{M}_{\tilde{V}_{k+1}}^{(i+1)}$ [Y, $r_0^{(i)} / \|r_0^{(i)}\|$], i.e.,

$$ (3.35) \qquad M_{[W_k, r_0^{(i)} / \|r_0^{(i)}\|]}^{(i+1)} [Y, r_0^{(i)} / \|r_0^{(i)}\|] = [M_{W_m}^{(i)} Y, \mathcal{M}_{k+1}^{(i+1)} (r_0^{(i)} / \|r_0^{(i)}\|)], $$

$$ (3.36) \qquad \tilde{M}_{\tilde{V}_{k+1}}^{(i+1)} [Y, r_0^{(i)} / \|r_0^{(i)}\|] = [\tilde{M}_{\tilde{V}_m}^{(i)} Y, \tilde{M}_{k+1}^{(i+1)} (r_0^{(i)} / \|r_0^{(i)}\|)]. $$

Using the assumption on the equivalent preconditioning matrix obtained at the end of the $i$th cycle, i.e., $M_{W_m}^{(i)} = \tilde{M}_{\tilde{V}_m}^{(i)}$, we have

$$ (3.37) \qquad M_{W_m}^{(i)} Y = \tilde{M}_{\tilde{V}_m}^{(i)} Y. $$

The fact that identical (possibly nonlinear) preconditioning operators are used in steps 10a and 10b of Algorithm 1 (see relation (3.1)) allows us to write

$$ (3.38) \qquad \mathcal{M}_{k+1}^{(i+1)} (r_0^{(i)} / \|r_0^{(i)}\|) = \tilde{\mathcal{M}}_{k+1}^{(i+1)} (r_0^{(i)} / \|r_0^{(i)}\|). $$

Relations (3.37) and (3.38) finally show the relation (3.34). Consequently from relations (3.32), (3.33), and (3.34) we deduce that there exists a nonsingular matrix $\hat{X}' \in \mathbb{C}^{(k+1) \times (k+1)}$ such that

$$ (3.39) \qquad Z_{k+1} = \tilde{Z}_{k+1} \hat{X}'. $$

This proves the relation proposed in (3.18). Since $T$ and $\tilde{T}$ are both triangular, we note that $\hat{X}' = \tilde{T} T^{-1}$ is also triangular.
Part IIIa: Steps 10a and 10b. We first show that $\tilde{v}_{k+2} = v_{k+2}$ by expressing these two quantities as a function of $r_0^{(i)}$ and $Y$.

**FGCRO-DR.** The Arnoldi relation (step 10a) yields $v_{k+2} = \bar{v}_{k+2}/||\bar{v}_{k+2}||$, where $\bar{v}_{k+2} = (I_n - v_{k+1}v_{k+1}^H)(I_n - V_kV_k^H)A\mathcal{M}_{k+1}^{(i+1)}(r_0^{(i)}/||r_0^{(i)}||)$. Since from Proposition 2 $V_k^H r_0^{(i)} = 0$ in the $(i+1)$th cycle, $(I_n - v_{k+1}v_{k+1}^H)$ and $(I_n - V_kV_k^H)$ commute, and from Part IIa of the proof, the following expression can be derived:

$$ (3.40) \quad \bar{v}_{k+2} = \Pi_{V_{k+1}^\perp} A\mathcal{M}_{k+1}^{(i+1)}(r_0^{(i)}/||r_0^{(i)}||) = \Pi_{[Y, r_0^{(i)}/||r_0^{(i)}||]^\perp} A\mathcal{M}_{k+1}^{(i+1)}(r_0^{(i)}/||r_0^{(i)}||). $$

**FGMRES-DR.** The following expression for $\tilde{v}_{k+2} = \bar{v}_{k+2}/||\bar{v}_{k+2}||$ is obtained using Lemma 3.5:

$$ (3.41) \quad \tilde{v}_{k+2} = (I_n - \tilde{V}_{k+1} \tilde{V}_{k+1}^H) A\mathcal{M}_{k+1}^{(i+1)}(\tilde{v}_{k+1}) = \Pi_{[Y, r_0^{(i)}/||r_0^{(i)}||]^\perp} A\mathcal{M}_{k+1}^{(i+1)}(\Pi_{Y^\perp} r_0^{(i)}/||\Pi_{Y^\perp} r_0^{(i)}||). $$

Due to the assumption (3.8) of Theorem 3.6 we deduce from (3.40) and (3.41) that $\bar{v}_{k+2} = \eta_{i+1} \tilde{v}_{k+2}$ with $\eta_{i+1}$ positive, and therefore $v_{k+2} = \tilde{v}_{k+2}$.

**Part IIIb:** Steps 10a and 10b. In this part we continue the analysis of the Arnoldi procedure with flexible preconditioning and show that $v_{k+2+j} = \tilde{v}_{k+2+j}$ for $j = 1, \dots, m-k-1$.

For the case $j=1$, we introduce $\bar{v}_{k+3}$ and $\tilde{v}_{k+3}$ such that $v_{k+3} = \bar{v}_{k+3}/||\bar{v}_{k+3}||$ and $\tilde{v}_{k+3} = \bar{v}_{k+3}/||\bar{v}_{k+3}||$. The application of the Arnoldi procedure in both algorithms leads to

$$ \bar{v}_{k+3} = (I_n - v_{k+2}v_{k+2}^H)(I_n - V_{k+1}V_{k+1}^H) A\mathcal{M}_{k+2}^{(i+1)}(\bar{v}_{k+2}), $$

$$ \tilde{v}_{k+3} = (I_n - \tilde{v}_{k+2}\tilde{v}_{k+2}^H)(I_n - \tilde{V}_{k+1}\tilde{V}_{k+1}^H) A\mathcal{M}_{k+2}^{(i+1)}(\tilde{v}_{k+2}). $$

Thus from Parts II and IIIa of the proof we obtain that $v_{k+3}$ and $\tilde{v}_{k+3}$ are equal. The proof can then be completed by induction.

Results from Parts II and III justify relation (3.19), i.e., $[v_{k+2}, \dots, v_{m+1}] = [\tilde{v}_{k+2}, \dots, \tilde{v}_{m+1}]$. Consequently from Lemma 3.2, Lemma 3.4, and relation (3.1) we deduce relation (3.20). This finally shows the main relations (3.9) and (3.10) of Theorem 3.6 that are satisfied at the end of the $(i+1)$th cycle. □

### 3.3.1. First consequence of Theorem 3.6.

**COROLLARY 3.7.** If the same flexible preconditioning operators are used in both Arnoldi procedures (steps 10a and 10b of Algorithm 1) and if at each cycle $i$ there exists a real-valued positive coefficient $\eta_i$ such that

$$ \Pi_{[Y, r_0^{(i-1)}/||r_0^{(i-1)}||]^\perp} A\mathcal{M}_{k+1}^{(i)}(\Pi_{Y^\perp} r_0^{(i-1)}/||\Pi_{Y^\perp} r_0^{(i-1)}||) \\ = \eta_i \Pi_{[Y, r_0^{(i-1)}/||r_0^{(i-1)}||]^\perp} A\mathcal{M}_{k+1}^{(i)}(r_0^{(i-1)}/||r_0^{(i-1)}||), $$

or, equivalently (from relations (3.40) and (3.41)), such that $\tilde{v}_{k+2} = \eta_i \bar{v}_{k+2}$, FGCRO-DR and FGMRES-DR are algebraically equivalent.

**Proof.** We have already emphasized that $M_{W_m}^{(0)} = \tilde{M}_{\tilde{V}_m}^{(0)}$ in relation (3.3). In Theorem 3.6 we have analyzed the $(i+1)$th cycle of both algorithms assuming that $M_{W_m}^{(i)} = \tilde{M}_{\tilde{V}_m}^{(i)}$. First we have proved in Part IIb the relation (3.34), and second in Parts IIIa and IIIb that $[v_{k+2}, \dots, v_m] = [\tilde{v}_{k+2}, \dots, \tilde{v}_m]$ and $[z_{k+2}, \dots, z_m] = [\tilde{z}_{k+2}, \dots, \tilde{z}_m]$, respectively. Consequently the same equivalent preconditioner matrix is obtained at the end of the $(i+1)$th cycle, i.e., $M_{W_m}^{(i+1)}$ and $\tilde{M}_{\tilde{V}_m}^{(i+1)}$ can be chosen equal. We deduce that FGCRO-DR and FGMRES-DR are algebraically equivalent. □
**3.3.2. About GCRO-DR and GMRES-DR.** We propose a second consequence of Theorem 3.6 analyzed now with a fixed preconditioning matrix M.

**COROLLARY 3.8.** *When a fixed right preconditioner is used, the GCRO-DR and GMRES-DR methods sketched in Algorithm 1 are algebraically equivalent.*

*Proof.* We denote by M the fixed right preconditioning operator. A straightforward reformulation of Lemma 3.3 in this context leads to the relation $Z_m = MW_m$ in GCRO-DR. Exploiting now Lemma 2.3 allows us to derive the following relation, which holds during the $(i+1)$th cycle:

$$AMW_m P_k = AMY = [Y, r_0^{(i)}] \begin{bmatrix} \operatorname{diag}(\theta_1, \dots, \theta_k) \\ \beta_{1 \times k} \end{bmatrix}.$$

Thus

$$(3.42) \qquad \Pi_{[Y, r_0^{(i)}] \perp AMY} = 0.$$

Due to (3.42) and Part IIIa of the proof of Theorem 3.6 we deduce the following development:

$$\begin{align*}
\bar{v}_{k+2} &= \Pi_{[Y, r_0^{(i)}] \perp AM} (r_0^{(i)} - Y Y^{\dagger} r_0^{(i)}), \\
\bar{v}_{k+2} &= \Pi_{[Y, r_0^{(i)}] \perp AM} \Pi_{Y \perp r_0^{(i)}}, \\
\bar{v}_{k+2} &= \tilde{\bar{v}}_{k+2}.
\end{align*}$$

By induction it is possible to deduce the rest of the proof regarding $\bar{v}_{k+j}, j > 2$. Using $\mathrm{range}(\tilde{V}_{k+1}) = \mathrm{range}(V_{k+1})$ obtained in Part IIa we deduce that

$$(3.43) \qquad \mathrm{range}(\tilde{V}_m) = \mathrm{range}(V_m) = \mathrm{range}(W_m).$$

A straightforward reformulation of Lemma 3.2 leads to the relation $\tilde{Z}_m = M\tilde{V}_m$ in GMRES-DR. From relation (3.43) we finally deduce that

$$\mathrm{range}(\tilde{Z}_m) = \mathrm{range}(Z_m).$$

Consequently the minimization problem $\min \|r_0^{(i)} - AZ_m y\|$ leads to the same solution for both algorithms at each cycle: GCRO-DR and GMRES-DR sketched in Algorithm 1 are thus algebraically equivalent. □

**3.3.3. A numerical illustration.** In this section we intend to illustrate the results shown in sections 3.3.1 and 3.3.2 on a simple numerical example. We consider a real symmetric positive definite matrix $A = Q D Q^T$ of size 200 with $Q$ orthonormal and $D$ diagonal with entries ranging from $10^{-4}$ to 1. The spectrum of $A$ contains eigenvalues of small magnitude,<fnref n="1" /> and consequently the use of deflation techniques should improve the convergence rate of Krylov subspace methods if the harmonic Ritz values of smallest modulus are taken into account. In this experiment we consider a polynomial preconditioner represented by two iterations of unpreconditioned GMRES for the solution of $Ax = b$ with $b$ given by $b = \frac{Ae}{\|Ae\|_2}$ ($e \in \mathbb{R}^{200}$ denoting the vector with all components equal to one) starting from a zero initial guess. Figure 3.1 shows the histories of convergence of various flexible methods minimizing

<footnote n="1">¹The eigenvalues of $A$ are logarithmically spaced $(10^{-4}, 10^{-3}, 10^{-2})$ and linearly distributed between 0.02 and 1 with step 1/200.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 3.1. Convergence histories of different flexible methods applied to $Ax = b$, where $A \in \mathbb{R}^{200 \times 200}$ is symmetric positive definite with some eigenvalues of small magnitude.</figcaption></figure>

<caption>TABLE 3.1</caption>

Scalar product $v_{k+2}^T \tilde{v}_{k+2}$ during the first five cycles of FGCRO-DR(10,6) when solving the linear system considered in section 3.3.3 using a variable preconditioner (two iterations of GMRES) and a fixed right preconditioner (diagonal preconditioning).

<table><thead><tr><th>Cycle index</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th></tr></thead><tbody><tr><td>Variable preconditioner</td><td>0.92</td><td>0.89</td><td>0.45</td><td>0.90</td><td>0.90</td></tr><tr><td>Fixed right preconditioner</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td><td>1.00</td></tr></tbody></table>

over a subspace of same dimension, i.e., FGMRES(10), FGMRES-DR(10,6), FGCRO-DR(10,6), respectively, and full flexible GMRES with such a variable preconditioner. Flexible methods with deflated restarting are found to be efficient since they are close to the full flexible GMRES method in terms of performances. We also remark that the convergence histories of FGCRO-DR(10,6) and FGMRES-DR(10,6) are different. According to Corollary 3.7 we compute the scalar product of $v_{k+2}$ and $\tilde{v}_{k+2}$ (which are both vectors of unit norm) to determine the cosine of the angle between these two vectors. The values are reported in Table 3.1 for the first five cycles. With such a variable preconditioner it is found that the methods are not equivalent in the first cycle already since the collinearity condition between $v_{k+2}$ and $\tilde{v}_{k+2}$ is not fulfilled. The situation is similar during the following cycles, which explains why different convergence histories for FGMRES-DR(10,6) and FGCRO-DR(10,6) observed in Figure 3.1 are obtained in such a case. As expected from section 3.3.2, if a fixed right preconditioner is used, the convergence histories of GMRES-DR(10,6) and FGCRO-DR(10,6) are found to be exactly the same (results not shown here). In such a case $v_{k+2}$ and $\tilde{v}_{k+2}$ fulfill the collinearity condition; this is confirmed in Table 3.1 when a diagonal preconditioning is used.
4. Further features of FGCRO-DR(m, k). In this section we first compare FGCRO-DR(m, k) with FGMRES-DR(m, k) presented in Algorithm 1 from both a computational and a storage point of view. Then we detail how subspace recycling can be used in FGCRO-DR(m, k) when solving a sequence of linear systems.

4.1. Computational cost. We first analyze the computational cost related to the generalized eigenvalue problem to deduce harmonic Ritz information and then detail the total cost of the proposed method.

4.1.1. Harmonic Ritz information. The generalized eigenvalue problem (3.22) can also be written as

$$ (4.1) \qquad \bar{H}_m^H \bar{H}_m y = \theta \bar{H}_m^H V_{m+1}^H W_m y. $$

Since $W_m = [W_{k+1}, v_{k+2}, \dots, v_m]$, $V_{m+1}^H W_m$ can be decomposed at the end of the cycle as

$$ (4.2) \qquad V_{m+1}^H W_m = \begin{bmatrix} V_{k+1}^H & W_{k+1} & 0_{(k+1)\times(m-k-1)} \\ 0_{(m-k-1)\times(k+1)} & I_{m-k-1} & 0_{1\times(m-k-1)} \\ 0_{1\times(k+1)} & 0_{1\times(m-k-1)} & 0_{(k+1)\times(m-k-1)} \end{bmatrix}, $$

where the structure of the $(k+1) \times (k+1)$ block $V_{k+1}^H W_{k+1}$ is as follows:

$$ V_{k+1}^H W_{k+1} = \begin{bmatrix} V_k^H & W_k & V_k^H & w_{k+1} \\ v_{k+1}^H & W_k & v_{k+1}^H & w_{k+1} \end{bmatrix} = \begin{bmatrix} V_k^H & W_k & 0_{k \times 1} \\ v_{k+1}^H & W_k & 1 \end{bmatrix}. $$

$V_k^H W_k$ is a $k \times k$ matrix that satisfies the following relation at the end of the $i$th cycle:

$$ (V_k^H W_k)^{(i)} = Q^H (V_{m+1}^H W_m)^{(i-1)} P_k R^{-1}, $$

where the superscript is related to the cycle index. Thus storing the $(m+1) \times m$ matrix $(V_{m+1}^H W_m)^{(i-1)}$ can be used to obtain $(V_k^H W_k)^{(i)}$ at a cost that is independent of $n$. Next we analyze how to compute efficiently $v_{k+1}^H W_k$ during the $i$th cycle. From relation (2.18) shown in Lemma 2.3 and Proposition 2, respectively, we deduce the relation

$$ (4.3) \qquad v_{k+1}^H V_k R = v_{k+1}^H W_k R \operatorname{diag}(\theta_1, \dots, \theta_k) + v_{k+1}^H r_0^{(i-1)} \beta_{1 \times k}. $$

Due to Proposition 2 and the definition of $v_{k+1}$, we have $v_{k+1}^H V_k = 0$. Thus we finally obtain that

$$ (4.4) \qquad v_{k+1}^H W_k = - \|(c - \bar{H}_m y^*)^{(i-1)}\|_2 \beta_{1 \times k} (R \operatorname{diag}(\theta_1, \dots, \theta_k))^{-1}, $$

where $\beta_{1 \times k}$ is obtained from relation (2.19), which does only involve projected quantities. This allows us to deduce $v_{k+1}^H W_k$ at a cost independent of $n$. From this development we draw two important consequences from a computational point of view. First, $(V_{m+1}^H W_m)^{(i)}$ can be obtained recursively at a cost that is independent of the problem size $n$. Second, storing $W_m$ (which would represent $m$ additional vectors of size $n$) is not mandatory; only $V_{m+1}^H W_m$—matrix of size $(m+1) \times m$—is required.
<caption>TABLE 4.1</caption>

Computational cost of a generic cycle of FGMRES-DR(m, k) and FGCRO-DR(m, k). C represents the total cost of FGCRO-DR(m, k) and corresponds to $C = (m - k)(op_M + op_A) + n(2(m + 1)k + 1 + 2mk + (m - k)(2m + 2k + 6))$.

<table><thead><tr><th>Computation of</th><th>FGMRES-DR(m, k)</th><th>FGCRO-DR(m, k)</th></tr></thead><tbody><tr><td>$V_m(:, 1:k+1)$</td><td>$2n(m+1)(k+1)$</td><td>$2n(m+1)k+n$</td></tr><tr><td>$Z_m(:, 1:k)$</td><td>$2nmk$</td><td>$2nmk$</td></tr><tr><td>$V_m(:, k+2:m+1)$</td><td>$(m-k)op_A + n(m-k)(2m+2k+5)$</td><td>$(m-k)op_A + n(m-k)(2m+2k+6)$</td></tr><tr><td>$Z_m(:, k+1:m)$</td><td>$(m-k)op_M$</td><td>$(m-k)op_M$</td></tr><tr><td>Total cost</td><td>$C+n(m+k+1)$</td><td>$C$</td></tr></tbody></table>

### 4.1.2. Cost of a cycle.
We summarize in Table 4.1 the main computational costs associated with each generic cycle of FGMRES-DR(m, k) and FGCRO-DR(m, k). In FGCRO-DR(m, k), an Arnoldi method based on the modified Gram–Schmidt procedure has been assumed.<fnref n="2" /> We have included only the costs proportional to the size of the original problem n which is supposed to be much greater than m and k. We denote by $op_A$ and $op_M$ the floating point operation counts for the matrix-vector product and the preconditioner application, respectively.

The generalized eigenvalue problem in FGCRO-DR(m, k) has been ignored in Table 4.1 since it can be performed at a cost independent of n as outlined in section 4.1.1. Furthermore the computation of c (required at step 12 of Algorithm 1) has not been included in Table 4.1 since in both methods it can be obtained at a cost independent of n (see Proposition 3 in [15] for FGMRES-DR). From Table 4.1 we deduce that FGCRO-DR(m, k) requires slightly fewer operations per cycle than FGMRES-DR(m, k).

### 4.2. Storage requirements.
We consider only the storage proportional to the size of the original problem n. Similarly, as in FGMRES-DR(m, k) (see [15, section 3.2.2]), if the matrix multiplications $V_{m+1}Q$ and $Z_m P_k R^{-1}$ at steps 8a and 9a of Algorithm 1 are performed in place (i.e., overwriting $V_k$ and $Z_k$, respectively), FGCRO-DR(m, k) requires only the storage of $Z_m$ and $V_{m+1}$, which corresponds to $(2m+1)$ vectors of length n. The same storage cost is needed in FGMRES-DR(m, k) as detailed in [15].

### 4.3. Solution of sequence of linear systems.
As advocated in [22], GCRO-DR(m, k) is suited for the solution of a sequence of slowly changing linear systems defined as $A^l x^l = b^l$ where both the matrix $A^l \in \mathbb{C}^{n \times n}$ and the right-hand side $b^l \in \mathbb{C}^n$ change from one system to the next, and the linear systems may typically not be available simultaneously. Next, we analyze how subspace recycling can be used in FGCRO-DR(m, k). We suppose that FGCRO-DR(m, k) has been applied for the solution of a given linear system (indexed by $s-1$) in this sequence and that appropriate subspaces to be recycled, $Z_k^{s-1}$ and $W_k^{s-1}$, have been selected during a given cycle.

<footnote n="2">²In FGCRO-DR(m, k) (step 10a of Algorithm 1) the action of $(I_n - V_k V_k^H)$ requires $\sum_{j=k+1}^m (4nk + n)$ operations, the Arnoldi method based on modified Gram–Schmidt requires $\sum_{j=k+1}^m \sum_{i=k+1}^j (4n)$ operations, whereas norm computation and normalization cost $\sum_{j=k+1}^m (3n)$ operations. In FGMRES-DR(m, k) (step 10b of Algorithm 1) the Arnoldi method based on modified Gram–Schmidt requires $\sum_{j=k+1}^m \sum_{i=1}^j (4n)$ operations due to maintaining orthogonality to $V_k$, whereas norm computation and normalization cost $\sum_{j=k+1}^m (3n)$ operations.</footnote>
As explained in Proposition 2, the relations $A^{s-1}Z_k^{s-1} = V_k^{s-1}$ with $V_k^{s-1}V_k^{s-1} = I_k$ and range($W_k^{s-1}$) = range($V_k^{s-1}$) are then supposed to hold. Proposition 3 details how to consider subspace recycling in the initial phase of FGCRO-DR($m, k$), when solving the new linear system $A^s x^s = b^s$ with $x_0$ as an initial guess.

**PROPOSITION 3.** Suppose that $Z_k^{s-1}$ and $W_k^{s-1}$ are defined from solving a previous linear system $A^{s-1}x^{s-1} = b^{s-1}$ with FGCRO-DR($m, k$) and that $A^s x^s = b^s$ is the new linear system to be solved. In the initial phase of FGCRO-DR with subspace recycling, the relation $A^s Z_k^s = V_k^s$ with $V_k^{sH}V_k^s = I_k$ holds with matrices $V_k^s, Z_k^s \in \mathbb{C}^{n \times k}$ defined as

$$V_k^s = Q, \\ Z_k^s = Z_k^{s-1} R^{-1}$$

with $QR = A^s Z_k^{s-1}$, where $Q \in \mathbb{C}^{n \times k}$ has orthonormal columns and $R \in \mathbb{C}^{k \times k}$ is upper triangular. In addition we define $W_k^s \in \mathbb{C}^{n \times k}$ as $W_k^s = W_k^{s-1} R^{-1}$.

*Proof.* By using information related to the reduced QR factorization of $A^s Z_k^{s-1}$ and the relation $A^{s-1} Z_k^{s-1} = V_k^{s-1}$, respectively, we easily obtain

$$A^s Z_k^s = A^s Z_k^{s-1} R^{-1} = Q \\ = V_k^s.$$

Since $Q$ has orthonormal columns, $V_k^s$ satisfies $V_k^{sH}V_k^s = I_k$. Finally $W_k^s = W_k^{s-1} R^{-1}$ is imposed to make sure that the relation shown in Lemma 3.3 will hold at the end of the initial phase of FGCRO-DR($m, k$) with subspace recycling. □

In the case of a sequence where only the right-hand sides are changing, we note that the reduced QR factorization (step 3 in Algorithm 2) is not required. The complete construction of the initial generation of subspaces $V_{m+1}^s, Z_m^s, W_m^s$ is sketched in Algorithm 2. Once $V_{m+1}^s, Z_m^s$, and $W_m^s$ have been obtained, the main cycle of FGCRO-DR($m, k$) (lines 4 to 15 of Algorithm 1) can be applied straightforwardly.

**Algorithm 2.** Initial generation of $V_{m+1}^s, Z_m^s$, and $W_m^s$ when subspace recycling is used to solve $A^s x^s = b^s$.

1: Suppose that $V_k^{s-1}, Z_k^{s-1}$ and $W_k^{s-1}$ are defined from solving a previous linear system $A^{s-1}x^{s-1} = b^{s-1}$ and satisfy $A^{s-1}Z_k^{s-1} = V_k^{s-1}$ with $V_k^{s-1}V_k^{s-1} = I_k$ and range($W_k^{s-1}$) = range($V_k^{s-1}$)

2: $r_0 = b^s - A^s x_0$

3: $Q R = A^s Z_k^{s-1}$

4: $V_k^s = Q$

5: $Z_k^s = Z_k^{s-1} R^{-1}$

6: $W_k^s = W_k^{s-1} R^{-1}$

7: $x_0^{(0)} = x_0 + Z_k^s V_k^{sH} r_0$

8: $r_0^{(0)} = r_0 - V_k^s V_k^{sH} r_0$

9: Apply $m-k$ flexible preconditioned Arnoldi steps with $(I_n - V_k^s V_k^{sH})A^s$ and $v_{k+1}^s = r_0^{(0)} / \|r_0^{(0)}\|$ such that $(I_n - V_k^s V_k^{sH})A^s [z_{k+1}^s, \dots, z_m^s] = [v_{k+1}^s, \dots, v_{m-k}^s]$ $\bar{H}_{m-k}$ with $z_j^s = M_j^{(i)}(v_j^s)$

10: $d^* = \arg \min_{d \in Z_m^s} \|r_0^{(0)} - A^s d\|$, $x_0^{(1)} = x_0^{(0)} + d^*$, $r_0^{(1)} = b^s - A^s x_0^{(1)}$

11: $W_m^s = [W_k^s \quad V_m^s(1:n, k+1:m)]$
<caption>TABLE 4.2</caption>

Solution of a d-dimensional elliptic partial differential equation problem on a 16<sup>d</sup> grid with homogeneous Dirichlet boundary conditions (d = 2, ..., 5). Shown are the total number of matrix-vector products (#Mvp) required to solve a sequence of twelve linear systems with different flexible methods. The variable preconditioner is based on four iterations of unpreconditioned GMRES. The stopping criterion corresponds to a reduction of six orders of magnitude of the normalized residual in the Euclidean norm. Harmonic Ritz values of smallest modulus have been considered when deflating.

<table><thead><tr><th>Grid</th><th>16<sup>2</sup></th><th>16<sup>3</sup></th><th>16<sup>4</sup></th><th>16<sup>5</sup></th></tr></thead><tbody><tr><td>Problem size (n)</td><td>(225)</td><td>(3375)</td><td>(50625)</td><td>(759375)</td></tr><tr><td>Method</td><td>#Mvp</td><td>#Mvp</td><td>#Mvp</td><td>#Mvp</td></tr><tr><td>FGMRES(20)</td><td>972</td><td>1176</td><td>1272</td><td>1128</td></tr><tr><td>FGMRES-DR(20,10)</td><td>732</td><td>948</td><td>1020</td><td>876</td></tr><tr><td>FGCRO-DR(20,10) (no recycling)</td><td>732</td><td>948</td><td>1020</td><td>876</td></tr><tr><td>FGCRO-DR(20,10) (with recycling)</td><td>457</td><td>541</td><td>547</td><td>529</td></tr></tbody></table>

Subspace recycling can thus be easily used in FGCRO-DR(*m*, *k*) to solve sequences of linear systems.

**4.3.1. A numerical illustration.** As a numerical illustration we consider sequences of linear systems arising from the finite difference discretization of multi-dimensional elliptic partial differential equations (isotropic Laplace operator) posed on the [0, 1]<sup>d</sup> hypercube with homogeneous Dirichlet boundary conditions. These sequences correspond to situations where only the right-hand sides are changing for a given dimension *d*. An efficient solution method is of primary interest in certain applications related to, e.g., financial engineering, molecular biology, or quantum dynamics [5, 6]. In the numerical experiments reported here (performed in MATLAB) we have used second order finite difference discretization schemes leading to sparse matrices with at most 2*d*+1 nonzero elements per row. We analyze the performances of various flexible methods used with four iterations of unpreconditioned GMRES as a preconditioner. This polynomial preconditioner is a variable nonlinear function which thus requires a flexible Krylov subspace method as an outer method [28]. Table 4.2 collects the number of matrix-vector products of some flexible methods minimizing over a subspace of the same dimension, i.e., FGMRES(20), FGMRES-DR(20,10), FGCRO-DR(20,10), and FGCRO-DR(20,10) with subspace recycling, respectively. Using deflation helps to improve the convergence rate of FGMRES in this application since a reduction of approximately 20% to 25% in terms of matrix-vector products is obtained for FGMRES-DR(20,10) independently of the dimension *d*. FGCRO-DR(20,10) leads to numbers of matrix-vector products which are similar to FGMRES-DR(20,10) although the convergence histories are found to be different. Finally, using both deflation and recycling in FGCRO-DR leads to a significant decrease in terms of matrix-vector products. A reduction in the range of 40% to 45% is indeed obtained versus another flexible Krylov subspace method with deflated restarting (FGMRES-DR(*m*, *k*)). This can be considered as a primary advantage over FGMRES-DR(*m*, *k*) since FGMRES-DR(*m*, *k*) does not allow subspace recycling. It nicely extends to the flexible setting the advantage of GCRO-DR versus GMRES-DR previously illustrated in [22]. We note that the resulting method is factorization free and mostly relies on matrix-vector products, a nice feature if distributed memory platforms are targeted to address numerical problems of larger size in higher dimension.

**5. Conclusion and perspectives.** In this paper we have studied a new minimum residual norm subspace method with deflated restarting that allows flexible
preconditioning based on the GCRO subspace method. The resulting method, named FGCRO-DR, has been presented together with FGMRES-DR, a recently proposed algorithm of the same family but based on the GMRES subspace method. A theoretical comparison analysis of both algorithms has been performed in section 3, where Theorem 3.6—the main result of this paper—proves the algebraic equivalence of FGMRES-DR and FGCRO-DR if a certain collinearity condition holds at each cycle. Corollary 3.8 has also proved that GMRES-DR and GCRO-DR are algebraically equivalent when a fixed right preconditioner is used. Furthermore we have carefully analyzed the computational cost of a given cycle of FGCRO-DR and have shown that FGCRO-DR is nearly as expensive as FGMRES-DR in terms of operations. FGCRO-DR offers the additional advantage of being suitable for the solution of sequences of slowly changing linear systems (where both the matrix and right-hand side can change) through subspace recycling.

In [8] variants of FGCRO-DR have been proposed which only differ in the formulation of the projected generalized eigenvalue problem. In future work we plan to investigate the numerical properties of these variants on realistic problems of large size for both single and multiple left- or right-hand side situations. Of interest are applications related to, e.g., steady or unsteady simulations of nonlinear equations [7] or stochastic finite element methods [12, 33] in three dimensions where variable preconditioning using approximate solvers has to be usually considered. We also note that when all right-hand sides are available simultaneously and when the matrix is fixed, block subspace methods may be also suitable. Thus a perspective could be to propose a block variant of FGCRO-DR.

**Acknowledgments.** We would like to thank the referees and the associate editor for their careful readings and valuable suggestions that helped us to improve our manuscript significantly. We also thank Iain S. Duff and Xavier Pinel for fruitful discussions and comments. The first author would like to acknowledge the warm welcome he received at CERFACS, in the Parallel Algorithms Team, during his sabbatical leave where the initial work was completed.

## REFERENCES

[1] O. AXELSSON AND P. S. VASSILEVSKI, *A black box generalized conjugate gradient solver with inner iterations and variable-step preconditioning*, SIAM J. Matrix Anal. Appl., 12 (1991), pp. 625–644.

[2] O. AXELSSON, *Iterative Solution Methods*, Cambridge University Press, Cambridge, UK, 1994.

[3] J. BAGLAMA, D. CALVETTI, G. H. GOLUB, AND L. REICHEL, *Adaptively preconditioned GMRES algorithms*, SIAM J. Sci. Comput., 20 (1998), pp. 243–269.

[4] A. H. BAKER, E. R. JESSUP, AND T. MANTEUFFEL, *A technique for accelerating the convergence of restarted GMRES*, SIAM J. Matrix Anal. Appl., 26 (2005), pp. 962–984.

[5] G. BEYLKIN AND M. J. MOHLENKAMP, *Algorithms for numerical analysis in high dimensions*, SIAM J. Sci. Comput., 26 (2005), pp. 2133–2159.

[6] H. BIN ZUBAIR, C. W. OOSTERLEE, AND R. WIENANDS, *Multigrid for high-dimensional elliptic partial differential equations on non-equidistant grids*, SIAM J. Sci. Comput., 29 (2007), pp. 1613–1636.

[7] M. H. CARPENTER, C. VUIK, P. LUCAS, M. B. VAN GIJZEN, AND H. BIJL, *A General Algorithm for Reusing Krylov Subspace Information*. I. *Unsteady Navier-Stokes*, NASA/Technical Memorandum 2010-216190, NASA, Langley Research Center, 2010.

[8] L. M. CARVALHO, S. GRATTON, R. LAGO, AND X. VASSEUR, *A Flexible Generalized Conjugate Residual Method with Inner Orthogonalization and Deflated Restarting*, Technical Report TR/PA/10/10, CERFACS, Toulouse, France, 2010.

[9] E. DE STURLER, *Nested Krylov methods based on GCR*, J. Comput. Appl. Math., 67 (1996), pp. 15–41.
[10] E. DE STURLER, *Truncation strategies for optimal Krylov subspace methods*, SIAM J. Numer. Anal., 36 (1999), pp. 864-889.

[11] M. EIERMANN, O. G. ERNST, AND O. SCHNEIDER, *Analysis of acceleration strategies for restarted minimal residual methods*, J. Comput. Appl. Math., 123 (2000), pp. 261-292.

[12] M. EIERMANN, O. G. ERNST, AND E. ULLMANN, *Computational aspects of the stochastic finite element method*, Comput. Visual. Sci., 10 (2007), pp. 3-15.

[13] S. C. EISENSTAT, H. C. ELMAN, AND M. H. SCHULTZ, *Variational iterative methods for non-symmetric systems of linear equations*, SIAM J. Numer. Anal., 20 (1983), pp. 345-357.

[14] J. ERHEL, K. BURRAGE, AND B. POHL, *Restarted GMRES preconditioned by deflation*, J. Comput. Appl. Math., 69 (1996), pp. 303-318.

[15] L. GIRAUD, S. GRATTON, X. PINEL, AND X. VASSEUR, *Flexible GMRES with deflated restarting*, SIAM J. Sci. Comput., 32 (2010), pp. 1858-1878.

[16] R. B. MORGAN, *A restarted GMRES method augmented with eigenvectors*, SIAM J. Matrix Anal. Appl., 16 (1995), pp. 1154-1171.

[17] R. B. MORGAN, *Implicitly restarted GMRES and Arnoldi methods for nonsymmetric systems of equations*, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1112-1135.

[18] R. B. MORGAN, *GMRES with deflated restarting*, SIAM J. Sci. Comput., 24 (2002), pp. 20-37.

[19] Y. NOTAY AND P. S. VASSILEVSKI, *Recursive Krylov-based multigrid cycles*, Numer. Linear Algebra Appl., 15 (2008), pp. 473-487.

[20] Y. NOTAY, *Flexible conjugate gradients*, SIAM J. Sci. Comput., 22 (2000), pp. 1444-1460.

[21] C. C. PAIGE, B. N. PARLETT, AND H. A. VAN DER VORST, *Approximate solutions and eigenvalue bounds from Krylov subspaces*, Numer. Linear Algebra Appl., 2 (1995), pp. 115-134.

[22] M. L. PARKS, E. DE STURLER, G. MACKEY, D. D. JOHNSON, AND S. MAITI, *Recycling Krylov subspaces for sequences of linear systems*, SIAM J. Sci. Comput., 28 (2006), pp. 1651-1674.

[23] Y. SAAD AND M. H. SCHULTZ, *GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems*, SIAM J. Sci. Statist. Comput., 7 (1986), pp. 856-869.

[24] Y. SAAD, *A flexible inner-outer preconditioned GMRES algorithm*, SIAM J. Sci. Statist. Comput., 14 (1993), pp. 461-469.

[25] Y. SAAD, *Analysis of augmented Krylov subspace methods*, SIAM J. Matrix Anal. Appl., 18 (1997), pp. 435-449.

[26] Y. SAAD, *Iterative Methods for Sparse Linear Systems*, 2nd ed., SIAM, Philadelphia, 2003.

[27] V. SIMONCINI AND D. B. SZYLD, *Flexible inner-outer Krylov subspace methods*, SIAM J. Numer. Anal., 40 (2003), pp. 2219-2239.

[28] V. SIMONCINI AND D. B. SZYLD, *Theory of inexact Krylov subspace methods and applications to scientific computing*, SIAM J. Sci. Comput., 25 (2003), pp. 454-477.

[29] V. SIMONCINI AND D. B. SZYLD, *Recent computational developments in Krylov subspace methods for linear systems*, Numer. Linear Algebra Appl., 14 (2007), pp. 1-59.

[30] G. L. G. SLEIJPEN AND H. A. VAN DER VORST, *A Jacobi-Davidson iteration method for linear eigenvalue problems*, SIAM J. Matrix Anal. Appl., 17 (1996), pp. 401-425.

[31] D. B. SZYLD AND J. A. VOGEL, *FQMR: A flexible quasi-minimal residual method with inexact preconditioning*, SIAM J. Sci. Comput., 23 (2001), pp. 363-380.

[32] A. TOSELLI AND O. WIDLUND, *Domain Decomposition Methods: Algorithms and Theory*, Springer Ser. Comput. Math. 34, Springer, New York, 2004.

[33] E. ULLMANN, *Solution Strategies for Stochastic Finite Element Discretizations*, Ph.D. thesis, Technische Universität Bergakademie, Freiberg, Germany, 2008.

[34] H. A. VAN DER VORST AND C. VUIK, *GMRESR: A family of nested GMRES methods*, Numer. Linear Algebra Appl., 1 (1994), pp. 369-386.