A Survey of Low-Rank Updates of Preconditioners for Sequences of Symmetric Linear Systems

Abstract

The aim of this survey is to review some recent developments in devising efficient preconditioners for sequences of symmetric positive definite (SPD) linear systems $A_k{\mathbf{x}}_k={\mathbf{b}}_k,k=1,\dots$ arising in many scientific applications, such as discretization of transient Partial Differential Equations (PDEs), solution of eigenvalue problems, (Inexact) Newton methods applied to nonlinear systems, rational Krylov methods for computing a function of a matrix. In this paper, we will analyze a number of techniques of updating a given initial preconditioner by a low-rank matrix with the aim of improving the clustering of eigenvalues around 1, in order to speed-up the convergence of the Preconditioned Conjugate Gradient (PCG) method. We will also review some techniques to efficiently approximate the linearly independent vectors which constitute the low-rank corrections and whose choice is crucial for the effectiveness of the approach. Numerical results on real-life applications show that the performance of a given iterative solver can be very much enhanced by the use of low-rank updates.

1. Introduction

The solution of sequences of large and sparse linear systems $A_k{\mathbf{x}}_k={\mathbf{b}}_k,k=1,\dots$, is a problem arising in many realistic applications including time-dependent solution of Partial Differential Equations (PDEs), computation of a part of the spectrum of large matrices, solution to nonlinear systems by the Newton methods and its variants, evaluation of a matrix function on a vector, by rational Krylov methods. The large size and sparsity of the matrices involved make iterative methods as the preeminent solution methods for these systems and therefore call for robust preconditioners. Standard—full purpose—preconditioners for SPD matrices, such as the Incomplete Cholesky (IC) factorization or approximate inverses [1,2], as well as Multigrid Methods [3] are aimed at clustering eigenvalues of the preconditioned matrices around one.

In this paper, we will analyze a number of techniques of updating a given IC preconditioner (which we denote as $P_0$ in the sequel) by a low-rank matrix with the aim of further improving this clustering. The most popular low-rank strategies are aimed at removing the smallest eigenvalues (deflation, see e.g., References [4,5,6,7] for the nonsymmetric case) or at shifting them towards the middle of the spectrum. The low-rank correction is based on a (small) number of linearly independent vectors whose choice is crucial for the effectiveness of the approach. In many cases these vectors are approximations of eigenvectors corresponding to the smallest eigenvalues of the preconditioned matrix $P_{0}A$ [8,9].

We will also review some techniques to efficiently approximate these vectors when incorporated within a sequence of linear systems all possibly having constant (or slightly changing) coefficient matrices. Numerical results concerning sequences arising from discretization of linear/nonlinear PDEs and iterative solution of eigenvalue problems show that the performance of a given iterative solver can be very much enhanced by the use of low-rank updates.

The issue of updating (not necessarily by a low-rank matrix) a given preconditioner for solving sequences of linear systems has been considered by many authors, among which I would mention Reference [10,11], in which the authors modify a given LU factorization to precondition many systems in the sequence, and the work related to the solution of sequences of shifted linear systems, namely systems with matrices of the form $A_k=A_0+{\alpha }_{k}I$. Among the others we mention the works in Reference [12,13] and also Reference [7] when $A_k$ is nonsymmetric.

The synopsis of the paper is as follows: In Section 2 we will discuss various low-rank updates, namely, the deflation technique, the tuning strategy, described in an original fashion on the basis of a revisited secant condition; and the spectral preconditioners. In Section 3 we will be concerned with the computational complexity of the low-rank updates. Section 4 discusses the choice of the vectors which form the low-rank matrices while Section 5 will address the issue to cheaply assessing some of the leftmost eigenpairs of the coefficient matrix. In Section 6 we will give some hint on the possibility of extending these techniques to sequences of nonsymmetric linear systems. Section 7 is devoted to reporting some, partly new, numerical experiments in the framework of discretization of transient PDES and solution to eigenvalue problems. The conclusions are drawn in Section 8.

Notation. Throughout the paper we use the Euclidean norm for vectors and the induced spectral norm for matrices.

2. Low-Rank Updates

2.1. Deflation

The idea to use deflation to accelerate the Conjugate Gradient (CG) method goes back to the middle of the 1980s. The most cited article on this subject is a paper by Nicolaides [4] who developed a deflated Conjugate Gradient method based on a three-term recurrence. Other similar formulations are presented in Reference [14,15]. Some years later in Reference [6] a deflated Generalized Minimal Residual (GMRES) method was proposed by using an augmented basis. Rigorous bounds on the eigenvalues of the deflated-preconditioned matrices are given in Reference [16]. We finally mention a recent interesting survey on deflation methods in Reference [17].

We will follow here the formulation of the Deflated-CG algorithm developed in Reference [5]. Given a full column rank $n\times p$ matrix $W=\left[\begin{array}{ccc}{\mathbf{w}}_1& \dots & {\mathbf{w}}_p\end{array}\right]$ the A-orthogonal projector H is defined as

$$H=I-W{\left(W^{T}AW\right)}^{-1}{W}^{T}A,$$

(1)

which maintains the CG residuals orthogonal to the subspace spanned by $\{{\mathbf{w}}_1,\dots ,{\mathbf{w}}_p\}$, provided that also ${\mathbf{r}}_0$ satisfies $W^T{\mathbf{r}}_0=0$. The Deflated-CG algorithm is described in Algorithm 1.

Algorithm 1 Deflated CG.

Choose an $n\times p$ full column rank matrix W. Compute $\Pi =W^{T}AW$.

Choose ${\tilde{\mathbf{x}}}_0$; Set ${\mathbf{r}}_0=\mathbf{b}-A{\mathbf{x}}_0$ and ${\mathbf{x}}_0={\tilde{\mathbf{x}}}_0+W{\Pi }^{-1}{W}^T{\mathbf{r}}_0$.

Compute ${\mathbf{z}}_0=P_0{\mathbf{r}}_0$;

Set $k=0$, ${\mathbf{p}}_0={\mathbf{z}}_0-W{\Pi }^{-1}{W}^{T}A{\mathbf{z}}_0$

${\rho }_0={\mathbf{r}}_0^T{\mathbf{z}}_0$

repeat until convergence

${\alpha }_k={\displaystyle \frac{{\rho }_k}{{\mathbf{p}}_k^{T}A{\mathbf{p}}_k}}$

${\mathbf{x}}_{k+1}={\mathbf{x}}_k+{\alpha }_k{\mathbf{p}}_k$

${\mathbf{r}}_{k+1}={\mathbf{r}}_k-{\alpha }_{k}A{\mathbf{p}}_k$

$k=k+1$

${\mathbf{z}}_k=P_0{\mathbf{r}}_k$

${\rho }_k={\mathbf{r}}_k^T{\mathbf{z}}_k$

${\beta }_k={\displaystyle \frac{{\rho }_k}{{\rho }_{k-1}}}$

${\mathbf{p}}_{k+1}={\beta }_k{\mathbf{p}}_k+{\mathbf{z}}_k-W{\Pi }^{-1}{W}^{T}A{\mathbf{z}}_k$

end repeat.

Note that this algorithm is mathematically equivalent to the PCG method with preconditioner $P=H{P}_0{H}^T$ since

$${\mathbf{r}}_k^T{\mathbf{z}}_k={\mathbf{r}}_k^T{P}_0{\mathbf{r}}_k={\mathbf{r}}_k^{T}H{P}_0{H}^T{\mathbf{r}}_k,$$

being $W^T{\mathbf{r}}_k=0$ and therefore $H^T{\mathbf{r}}_k={\mathbf{r}}_k$. P is only symmetric positive semidefinite, however PCG could not breakdown as stated in the following result ([5], Theorem 4.2).

Theorem 1.

Let A be an SPD matrix and W an $n\times p$ matrix with full column rank. The Deflated-CG algorithm applied to the linear system $A\mathbf{x}=\mathbf{b}$ will not break down at any step. The approximate solution ${\mathbf{x}}^{\left(k\right)}$ is the unique minimizer of $\parallel {\mathbf{x}}^{\left(k\right)}-{\mathbf{x}}^{(\ast )}{\parallel }_A$ over the affine space ${\mathbf{x}}_0+span(W,{\mathcal{K}}_k\left({\mathbf{r}}_0\right))$.

The action of this deflated preconditioner is to put some of the eigenvalues of the preconditioner matrix to zero in fact:

$$PAW=H{P}_0{H}^{T}AW=0,$$

that is, all p columns of W are eigenvectors of $PA$ corresponding to the zero eigenvalue. Deflated CG guarantees that (at least in infinite precision arithmetic) all residuals satisfy $W^T{\mathbf{r}}_k=0$. If the columns of W are (approximate) eigenvectors of $P_{0}A$ corresponding to the smallest eigenvalues, then these eigenvalues are removed from the spectrum of $PA$ with obvious decrease of the condition number.

2.2. Tuning

The adjective tuned associated with a preconditioner was introduced in Reference [18] in the framework of iterative eigensolvers. In our context we redefine it in the following way:

Definition 1.

Given a preconditioner $P_0$ and a vector $\mathbf{w}$, a tuned preconditioner for matrix A is a matrix $P\equiv M^{-1}$ obtained by correcting $P_0$ by a rank one or rank two matrix depending on A and $\mathbf{w}$ and satisfying

$$PA\mathbf{w}=\mathbf{w}(M\mathbf{w}=A\mathbf{w}).$$

(2)

In this way matrix M agrees with A in the direction $\mathbf{w}$ and also the preconditioned matrix has the eigenvalue 1 corresponding to the eigenvector $\mathbf{w}$. This definition can be easily extended to multiple vectors:

Definition 2.

Given a preconditioner $P_0$ and an $n\times p$ matrix W with full column rank, a block tuned preconditioner for matrix A is a matrix P obtained by correcting $P_0$ by a low-rank matrix depending on A and W and satisfying

$$PAW=W.$$

(3)

In  Reference [19], the following single-vector tuned preconditioned for SPD linear systems was proposed (here $M_0=L{L}^T\approx A$ and $P=M^{-1}$ is computed by the Sherman-Morrison formula).

$$\begin{array}{ccc}\hfill \mathbf{u}=(A-M_0)\mathbf{w}& & \mathbf{z}=P_0\mathbf{u}=P_{0}A\mathbf{w}-\mathbf{w}\hfill \end{array}$$

(4)

$$\begin{array}{ccc}\hfill M=M_0+\frac{{\mathbf{uu}}^T}{{\mathbf{u}}^T\mathbf{w}}& & P=P_0-\frac{{\mathbf{zz}}^T}{{\mathbf{z}}^{T}A\mathbf{w}}\hfill \end{array}$$

(5)

It is easily checked that $PA\mathbf{w}=P_{0}A\mathbf{w}-\mathbf{z}=\mathbf{w}$.

Digression

This tuned preconditioner, as many others, has however an older derivation. Let us now consider the solution of a nonlinear system of equations

$$\mathbf{F}\left(\mathbf{x}\right)=0,\mathbf{F}:{\mathbb{R}}^n\to {\mathbb{R}}^n$$

by the Newton method

$$\begin{array}{ccc}\hfill F^{\prime }\left({\mathbf{x}}_k\right){\mathbf{s}}_k& =& -\mathbf{F}\left({\mathbf{x}}_k\right)\hfill \\ \hfill {\mathbf{x}}_{k+1}& =& {\mathbf{x}}_k+{\mathbf{s}}_k.\hfill \end{array}$$

(6)

Quasi-Newton methods construct a sequence of approximations of the Jacobians $B_k\approx F^{\prime }\left({\mathbf{x}}_k\right)$. Each $B_k$ is defined by a low-rank update of the previous matrix in the sequence $B_{k-1}$. Such sequences are also used to correct a given initial preconditioner to accelerate the iterative solution of systems like (6) [20,21,22,23,24].

Among those Quasi-Newton methods we recall three of the most commonly used one (note that ${\mathbf{y}}_k=\mathbf{F}\left({\mathbf{x}}_{k+1}\right)-\mathbf{F}\left({\mathbf{x}}_k\right)$):

Broyden’s method: $B_{k+1}=B_k+{\displaystyle \frac{({\mathbf{y}}_k-B_k{\mathbf{s}}_k){\mathbf{s}}_k^T}{{\mathbf{s}}_k^T{\mathbf{s}}_k}}$

SR1 (Symmetric Rank-1): $B_{k+1}=B_k+{\displaystyle \frac{({\mathbf{y}}_k-B_k{\mathbf{s}}_k){({\mathbf{y}}_k-B_k{\mathbf{s}}_k)}^T}{{({\mathbf{y}}_k-B_k{\mathbf{s}}_k)}^T{\mathbf{s}}_k}}$

BFGS (Named after Broyden, Fletcher, Goldfarb and Shanno, who all discovered the same formula independently in 1970, by different approaches) update: $B_{k+1}=B_k+{\displaystyle \frac{{\mathbf{y}}_k{\mathbf{y}}_k^T}{{\mathbf{y}}_k^T{\mathbf{s}}_k}}-{\displaystyle \frac{B_k{\mathbf{s}}_k{\mathbf{s}}_k^T{B}_k}{{\mathbf{s}}_k^T{B}_k{\mathbf{s}}_k}}$.

All these updates satisfy the so called secant condition which reads

$$B_{k+1}{\mathbf{s}}_k={\mathbf{y}}_k.$$

(7)

The first update is nonsymmetric, the SR1 and BFGS formulae define a symmetric sequence provided $B_0$ is so. Consider now the problem of updating a given preconditioner $P_0=M_0^{-1}$. We can use all the preceding equations after employing the following substitutions:

$${\mathbf{s}}_k⟶\mathbf{w},{\mathbf{y}}_k⟶A\mathbf{w},B_k⟶M_0,B_{k+1}⟶M,$$

which transform the secant condition into

$$M\mathbf{w}=A\mathbf{w},$$

(8)

that is the tuning property. We can then develop three different tuning strategies. For each of them, in

Table 1. Direct, inverse and block versions of the Broyden (nonsymmetric) update.

direct	$M=M_0+{\displaystyle \frac{(A-M_0)\mathbf{w}{\mathbf{w}}^T}{{\mathbf{w}}^T\mathbf{w}}}$	$M\mathbf{w}=A\mathbf{w}$
inverse	$P=P_0-{\displaystyle \frac{(P_{0}A\mathbf{w}-\mathbf{w}){\mathbf{w}}^T{P}_0}{{\mathbf{w}}^T{P}_{0}A\mathbf{w}}}$	$PA\mathbf{w}=\mathbf{w}$
block	$P=P_0-(P_{0}AW-W){\left(W^T{P}_{0}AW\right)}^{-1}{W}^T{P}_0$	$PAW=W$

,

Table 2. Direct, inverse and block versions of the SR1 update.

direct	$M=M_0+{\displaystyle \frac{\mathbf{u}{\mathbf{u}}^T}{{\mathbf{w}}^T\mathbf{u}}}$,	$\mathbf{u}=(A-M_0)\mathbf{w}$
inverse	$P=P_0-{\displaystyle \frac{\mathbf{z}{\mathbf{z}}^T}{{\mathbf{z}}^{T}A\mathbf{w}}}$,	$\mathbf{z}=P_{0}A\mathbf{w}-\mathbf{w}$
block	$P=P_0-Z{\left(Z^{T}AW\right)}^{-1}{Z}^T$,	$Z=P_{0}AW-W$

and

Table 3. Direct, inverse and block versions of the BFGS update.

direct	$M=M_0+{\displaystyle \frac{{(A\mathbf{w})}^{T}A\mathbf{w}}{{\mathbf{w}}^{T}A\mathbf{w}}}-{\displaystyle \frac{M_0\mathbf{w}{\mathbf{w}}^T{M}_0}{{{w}}^T{M}_0\mathbf{w}}}$
inverse	$P={\displaystyle \frac{\mathbf{w}{\mathbf{w}}^T}{{\mathbf{w}}^{T}A\mathbf{w}}}+\left(I-{\displaystyle \frac{\mathbf{w}{\mathbf{w}}^{T}A}{{\mathbf{w}}^{T}A\mathbf{w}}}\right)P_0\left(I-{\displaystyle \frac{A\mathbf{w}{\mathbf{w}}^T}{{\mathbf{w}}^{T}A\mathbf{w}}}\right)$
block	$P=W{\Pi }^{-1}{W}^T+H{P}_0{H}^T$  $\Pi =W^{T}AW$  (1) $H=I-W{\Pi }^{-1}{W}^{T}A$

, we write the “direct” formula (with $M,M_0$) the inverse representation (with $P,P_0$), by the Shermann-Morrison inversion formula, and develop the corresponding inverse block formula.

The Broyden tuning strategy must be used for nonsymmetric problems, the BFGS formula is well suited to accelerate the PCG method due to the following result:

Theorem 2.

The preconditioner P yielded by the BFGS update formula is SPD provided $P_0$ is so.

Proof. 

For every nonzero $\mathbf{x}\in {\mathbb{R}}^n$ we set $\mathbf{z}=H^T\mathbf{x}$ and $\mathbf{u}=W^T\mathbf{x}$. Then we have

$${\mathbf{x}}^{T}P\mathbf{x}={\left(W^T\mathbf{x}\right)}^T{\Pi }^{-1}\left(W^T\mathbf{x}\right)+{\mathbf{x}}^{T}H{P}_0{H}^T\mathbf{x}={\mathbf{u}}^T{\Pi }^{-1}\mathbf{u}+{\mathbf{z}}^T{P}_0\mathbf{z}\ge 0,$$

(9)

the last inequality holding since both ${\Pi }^{-1}$ and $P_0$ are SPD matrices. The inequality is strict since if $\mathbf{u}=0$ then $W^T\mathbf{x}=0$ and hence $\mathbf{z}=(I-AW{\Pi }^{-1}{W}^T)\mathbf{x}=\mathbf{x}\ne 0$. □

Finally, the SR1 update provides clearly symmetric matrices, upon symmetry of $P_0$. If in addition $P_0$ is SPD the new update P is also likely to be SPD as well, depending on the quality of the initial preconditioner $P_0$, as stated in the following result [25] (Section 3.2).

$${\lambda }_{min}\left(P\right)\ge {\lambda }_{min}\left(P_0\right)-\parallel P_0{A-I\parallel }^2\parallel {\left(W^{T}AW\right)}^{-1}\parallel ,$$

(10)

which has however a modest practical use as $\parallel {\left(W^{T}AW\right)}^{-1}\parallel$ may be large.

Whenever the columns of W approximate the smallest eigenvalues of $P_{0}A$ something more can be said. Assume that $P_{0}AW=W\Theta$, with $\Theta =\mathrm{diag}\left({\mu }_1,\dots ,{\mu }_p\right)$ the diagonal matrix with the eigenvalues of the initially preconditioned matrix. Assume further that ${\mu }_j<1,j=1,\dots ,p$. Since $P_{0}A$ is not symmetric but $P_0^{1/2}A{P}_0^{1/2}$ is so it follows that W must satisfy $W^T{P}_0^{-1}W=I$ and hence $W^{T}AW=\Theta$. In such a case we have that matrix

$$-Z^{T}AW={\left(W-P_{0}AW\right)}^{T}AW=(I-\Theta )W^{T}AW=(I-\Theta )\Theta$$

is obviously SPD and therefore so is P which is the sum of an SPD matrix ($P_0$) and a symmetric positive semidefinite matrix ($-Z{\left(Z^{T}AW\right)}^{-1}{Z}^T$).

The Broyden tuned preconditioner was used in Reference [26] to accelerate the GMRES for the inner linear systems within the Arnoldi method for eigenvalue computation. The SR1 preconditioner appeared in References [19,25] and also in Reference [27]. This low-rank update has also been employed to accelerate the PCG method in the solution of linear systems involving SPD Jacobian matrices. In Reference [28] some conditions are proved under which the SR1 update maintains the symmetric positive definiteness of a given initial preconditioner. The BFGS preconditioner was used (under the name of balancing preconditioner) in Reference [29], in Reference[30] for eigenvalue computation and also in Reference [31] to update the Constraint Preconditioner for sequences of KKT linear systems. It has been successfully employed (under the name of limited memory preconditioner) to accelerate sequences of symmetric indefinite linear systems in Reference [32].

We finally mention the work in Reference [33], where the author used this correction to update an Inexact Constraint Preconditioner, by damping the largest eigenvalues of the preconditioned matrices.

2.3. Spectral Preconditioners

A further strategy to improve the preconditioner’s efficiency by a low-rank update can be found in References [8,9,27,34,35]. A spectral preconditioner for SPD linear systems is defined as

$$P=P_0+W{\left(W^{T}AW\right)}^{-1}{W}^T,$$

with the leftmost eigenvectors of $P_{0}A$ as the columns of W. Denoted as ${\mu }_1,{\mu }_2,\dots ,{\mu }_p$ the smallest eigenvalues of $P_{0}A$ it is easily checked that matrix $PA$ has the same eigenvectors ${\mathbf{w}}_1,\dots ,{\mathbf{w}}_p$ as $P_{0}A$ with ${\mu }_1+1,\dots ,{\mu }_p+1$ as eigenvalues.

3. Implementation and Computational Complexity

We sketch below the most time-consuming computational tasks when using a low-rank update of a given preconditioner. Notation: with the symbol $O\left(f\right(n\left)\right)$ we mean that there are two constant $c_1$ and $c_2$ independent of n such that $c_{1}f\left(n\right)\le O\left(f\left(n\right)\right)\le c_{2}f\left(n\right)$.

We first evaluate the preliminary cost, which has to be done before starting he iterative process:

• All low-rank updates: Computation of $A_W=AW(O\left(np\right)$ operations).
• Broyden and SR1 only: Computation of $Z=P_0{A}_W-W$ (p applications of the initial preconditioner).
• All low-rank updates: Computation of $\Pi =W^T{A}_W(O\left(p^{2}n\right)$ operations).

The cost of this tasks is likely to be amortized during the PCG iterations of the various linear systems to be solved.

Regarding the most important cost per iteration, the deflated, Broyden, SR1 and spectral preconditioners behave in the same way, being all based on the following main tasks, in addition to the application of the initial preconditioner $P_0$ to a given vector.

• Multiplication of an $n\times p$ matrix times a vector ($O\left(np\right)$ operations)
• Solution of a small linear system with matrix $\Pi$. ($O\left(p^3\right)$ operations).
• Multiplication of a $p\times n$ matrix times a vector ($O\left(np\right)$ operations)

The BFGS preconditioner, as it needs to correct $P_0$ by a rank-$\left(2p\right)$ matrix, roughly doubles tasks 1–3 and hence the total updating cost at each iteration.

4. Choice of the Vectors $\left\{{\mathbf{w}}_{\mathit{j}}\right\}$

In principle every set of linearly independent vectors can be selected as columns of W. However, the optimal choice is represented by the eigenvectors of the preconditioned matrix $P_{0}A$ corresponding to the smallest eigenvalues. Approximation of these vectors as a by-product of the Conjugate Gradient method has been considered in various works. We will turn on this crucial issue later on.

We start by presenting some numerical results obtained computing the “true” 10 leftmost eigenvectors of $P_{0}A$ using the Matlab eigs command. To measure the sensitivity of the method to eigenvector accuracy we also use perturbed eigenvectors that is, satisfying $\parallel P_{0}A{\mathbf{w}}_j-{\mu }_j{\mathbf{w}}_j\parallel \approx \delta$. The coefficient matrix was chosen as the FD discretization of the Laplacian on an L-shaped domain obtained using the Matlab command

A = delsq(numgrid(‘L’,500));

which returns a sparse matrix of order $n=186003$. The linear system $A\mathbf{x}=\mathbf{b}$, with $\mathbf{b}$ a random uniformly distributed vector, was solved by the PCG method with various low-rank update techniques with $P_0=IC\left(0\right)$.

The results in

Table 4. Influence of accuracy of the eigenvectors of $P_{0}A$ on the efficiency of the low-rank preconditioners. PCG iterations with various updating strategies are reported.

	No Update	Tuned	Deflated	Spectral
exact	466	254	254	254
$\delta =0.01$	466	261	259	290
$\delta =0.05$	466	378	260	286

show that a very low accuracy on the eigenvectors (absolute residual of order of $\delta =0.01$) is sufficient to provide good preconditioning results. Also notice that the tuned preconditioner seems to be more sensitive to the parameter $\delta$.

It may also happen that some of the leftmost eigenpairs of the coefficient matrix are instead available. What if one uses these vectors as columns of W? The following

Table 5. Influence of accuracy of the eigenvectors of A on the efficiency of the low-rank preconditioners. PCG iterations with various updating strategies are reported.

	No Update	Tuned	Deflated	Spectral
exact	466	254	254	254
$\delta ={10}^{-3}$	466	296	296	297
$\delta =0.01$	466	362	361	369

gives a surprising answer: if they are computed to a good accuracy they provide a notable acceleration of convergence. Again we use either exact eigenvectors or vectors satisfying $\parallel A{\mathbf{w}}_j-{\lambda }_j{\mathbf{w}}_j\parallel \approx \delta$.

The conclusion is: yes we can use the leftmost eigenvectors of A instead of those of $P_{0}A$ without dramatically loosing efficiency. The reason for this behavior is connected to the action of preconditioners such as Incomplete Cholesky or approximate inverses. They usually leave almost unchanged the eigenvectors corresponding to the smallest eigenvalues (though increasing the latter ones).

Consider again the previous matrix and compare the eigenpairs $({\lambda }_j,{\mathbf{v}}_j)$ of A, with those of the preconditioned matrix with two choices of $P_0$: $IC(0)({\mu }_j^{P1},{\mathbf{v}}_j^{P1})$ and the IC factorization with droptol $=1e-2({\mu }_j^{P2},{\mathbf{v}}_j^{P2})$. As expected, in Figure 1 (left) we give evidence that the preconditioners shift bottom-up the smallest eigenvalues. However the angle between the corresponding eigenvectors, computed as

$$\angle \left(\mathbf{u},\mathbf{v}\right)=arccos\frac{{\mathbf{u}}^T\mathbf{v}}{\parallel \mathbf{u}\parallel ·\parallel \mathbf{v}\parallel },$$

and displayed in Figure 1, right, turns out to be close to zero showing that they are almost parallel.

Using Previous Solution Vectors

Computing a subset of the extremal eigenvectors of a given large and sparse matrix may reveal computationally costly. In the next Section we will develop some methods to save on this cost. However another possibility for the vectors ${\mathbf{w}}_j$ is to use a number of the solutions to the previous linear systems in the sequence. Assume now that the matrices in the sequence are allowed to (slightly) change, that is, we want to solve

$$A_k{\mathbf{x}}_k={\mathbf{b}}_k,k=1,\dots ,N.$$

In this case computing the leftmost eigenvectors of $A_k$ is inefficient since this preprocessing should be done at each linear system.

We then use, for the jth linear system ($1<j\le N$), the $j_{eff}=min\{j-1,p\}$ solutions to the previous linear systems as columns of $W\equiv W_j$, thus yielding:

$$W_j=\left[\begin{array}{ccc}{\mathbf{x}}_{j-j_{eff}}& \dots & {\mathbf{x}}_{j-1}\end{array}\right].$$

(11)

The idea beyond this choice, which is completely cost-free, comes from the reasoning that:

• A solution ${\mathbf{x}}_k$ of the k-th linear system has with high probability non-negligible components in the directions of the leftmost eigenvectors of $A_k$ since, if ${\displaystyle {\mathbf{b}}_k=\sum _{i=1}^n{\alpha }_i{\mathbf{u}}_i^{(k)}}$, then ${\displaystyle {\mathbf{x}}_k=\sum _{i=1}^n\frac{{\alpha }_i}{{\lambda }_i^{(k)}}{\mathbf{u}}_i^{(k)}}$, with the largest weights provided by the smallest eigenvalues.
• The leftmost eigenvector of $A_k,k=j_{eff},\dots j-1$ are good approximation of the leftmost eigenvectors of $A_j$.

To give experimental evidence of the efficiency of this approach we report some results in solving a sequence of 30 linear systems arising from the Finite Element/Backward Euler discretization of the branched transport equation [35] which gives raise to a sequence of $(2\times 2)$ block nonlinear systems in turn solved by the Inexact Newton method. After linearization and block Gaussian Elimination the problem reduces to the solution of a sequences of symmetric linear systems of the form

$${(A+B^{T}EB)}_k\mathbf{x}={\mathbf{b}}_k,k=1,\dots ,$$

where A is SPD, B is rectangular with full row rank, E is diagonal and possibly indefinite. We take 30 consecutive linear systems from this sequence, and solve them with the minimum residual (MINRES) iterative method using as the initial preconditioner the incomplete Cholesky factorization of A with drop tolerance ${10}^{-2}$ and the following updates

• BFGS, with update directions as in (11).
• Spectral preconditioner, with update directions as in (11).
• BFGS, with the accurate leftmost eigenvectors of $A_k$ as the update directions.

The results are reported in Figure 2 where we can appreciate reduction of the number of iterations provided by the low-rank update. The (cheaper) spectral preconditioner provides as good acceleration as the BFGS approach with update directions as in (11). Moreover the number of iterations is comparable to those obtained with the “ideal” eigenvectors (which are obtained by the Matlab eigs function).

5. Cost-Free Approximation of the Leftmost Eigenpairs

The most appropriate candidate vectors for the columns of W are anyway the eigenvectors of the preconditioned matrix corresponding to the smallest eigenvalues.

The Lanczos-PCG Connection

A simple way to recover some of the extremal eigenpairs of the preconditioned matrix is to exploit the so called Lanczos connection [36,37]. During the PCG method, preconditioned with $P_0$, it is possible to save the first m (scaled) preconditioned residuals as columns of a matrix $V_m$:

$$V_m=\left[\begin{array}{c}{\displaystyle \frac{P_0{\mathbf{r}}_0}{\sqrt{{\mathbf{r}}_0^T{P}_0{\mathbf{r}}_0}}},{\displaystyle \frac{P_0{\mathbf{r}}_1}{\sqrt{{\mathbf{r}}_1^T{P}_0{\mathbf{r}}_1}}},\dots ,{\displaystyle \frac{P_0{\mathbf{r}}_{m-1}}{\sqrt{{\mathbf{r}}_{m-1}^T{P}_0{\mathbf{r}}_{m-1}}}}\end{array}\right].$$

Matrix $V_m$ is such that $V_m^T{P}_0^{-1}{V}_m=I_m$, in view of the $P_0-$orthogonality of the residuals generated by the PCG method. The Lanczos tridiagonal matrix can be formed using the PCG coefficients ${\alpha }_k,{\beta }_k$:

$$T_m=\left[\begin{array}{cccc}{\displaystyle \frac{1}{{\alpha }_0}}& -{\displaystyle \frac{\sqrt{{\beta }_1}}{{\alpha }_0}}& & \\ -{\displaystyle \frac{\sqrt{{\beta }_1}}{{\alpha }_0}}& {\displaystyle \frac{1}{{\alpha }_1}}+{\displaystyle \frac{{\beta }_1}{{\alpha }_0}}& -{\displaystyle \frac{\sqrt{{\beta }_2}}{{\alpha }_1}}& \\ & & \ddots & \\ & & & -{\displaystyle \frac{\sqrt{{\beta }_{m-1}}}{{\alpha }_{m-2}}}\\ & & -{\displaystyle \frac{\sqrt{{\beta }_{m-1}}}{{\alpha }_{m-2}}}& {\displaystyle \frac{1}{{\alpha }_{m-1}}}+{\displaystyle \frac{{\beta }_{m-1}}{{\alpha }_{m-2}}}\end{array}\right]$$

Matrices $V_m$ and $T_m$ obey to the classical Lanczos relation that is:

$$V_m^{T}A{V}_m=T_m.$$

The eigensolution of $T_m$ provides the eigendecomposition $T_m{\Lambda }_m={\Lambda }_m{Q}_m$. Then, the eigenvectors corresponding to the p smallest eigenvalues are selected as p columns of matrix $Q_p$. Finally, the columns of $W_p=V_m{Q}_p$ are expected to approximate the p leftmost eigenvectors of $P_{0}A$.

This procedure can be implemented to a very little computational cost but it has a number of disadvantages: First, it requires the storage of m preconditioned residuals, Second, as the convergence for the Lanczos process to the smallest eigenvalues is relatively slow, it sometimes happens that PCG convergence takes place before eigenvector convergence. Third, some of the leftmost eigenpairs can be missing by the non-restarted Lanczos procedure.

To partially overcome these drawbacks in Reference [38] the authors suggest implementing a thick restart within the PCG iteration (without affecting it) and incrementally correct the eigenvector approximations during subsequent linear system solvers. The final set of directions is, in that paper, used only to deflate the initial vector and not to form a preconditioner (eigCG algorithm) in the solution of a sequence of hundreds of linear systems in lattice quantum chromodynamics.

A number of alternative ways to approximate eigenvectors during the PCG iteration can be found in the literature. In Reference [5], a technique is suggested to refine them during subsequent linear systems with the same matrix but different right-hand-sides by using the p smallest harmonic Ritz values. In this approach the vector stored within the PCG iteration are the conjugate directions ${\mathbf{p}}_k$ instead of the preconditioned residuals (they actually lie in the same Krylov subspace).

6. Sequences of Nonsymmetric Linear Systems

The techniques described so far have been especially developed to accelerate the PCG method for SPD matrices. In principle they can however be extended also to nonsymmetric linear systems. Even though eigenvalues of the preconditioned matrices are not always completely informative [39,40] on the convergence of Krylov methods for such problems, in any case shifting a small number of eigenvalues towards the middle of the spectrum is usually beneficial for iterative solvers. The spectral low-rank preconditioner (SRLU) was originally presented in Reference [9] for general matrices, with no diagonalizability assumptions, to accelerate the GMRES-DR method proposed by Morgan in Reference [41]. Regarding the tuning preconditioner, the Broyden update has been successfully used, for example, in Reference [26] to accelerate the GMRES for the inner linear systems within the Arnoldi method for eigenvalue computation. For shifted linear systems, the shift invariance property of Krylov subspace is very important for building the preconditioner (see e.g., References [7,42] for details). Regarding (non necessarily low-rank) updates of a given preconditioner, we mention Reference [43] in the framework of constrained optimization, and Reference [44], which showed that (roughly) one matrix factorization and a single Krylov subspace is enough for handling the whole sequence of shifted linear systems.

7. Numerical Results

The technique previously described are very powerful when applied to sequences of linear systems where the coefficient matrix remains the same. In this case it is possible to refine the eigenvalues of the preconditioned matrix during the first linear systems solutions up to machine precision [38]. However, in most situations this high accuracy is not really needed (see the discussion in Section 4).

In this section we will consider a particular case of sequences of linear systems with (slightly) changing matrices as

$$A_k\mathbf{x}={\mathbf{b}}_k,A_k=A_0+{\alpha }_k{A}_1,$$

which can be viewed as shifted linear systems. These sequences appear for example, in the Finite Difference or Finite Element discretization of transient PDEs or as an inner linear system within iterative eigensolvers. We will show that the low-rank techniques could be successfully employed in the acceleration of iterative methods in the solution of such linear systems.

7.1. Fe Discretization of a Parabolic Equation

Consider the parabolic PDE

$$\begin{array}{ccc}\hfill S_h\frac{\partial u(\overrightarrow{x},t)}{\partial t}& =& \nabla (K(\overrightarrow{x})\nabla u(\overrightarrow{x},t))+\mathrm{BCs}.\hfill \end{array}$$

u is the pressure of the fluid, $S_h$ the elastic storage, K the (SPD) hydraulic conductivity tensor. This equation is used to model water pumping in the Beijing plan to predict land subsidence in that area [45]. The 3D computational domain turns out to be very heterogeneous and geometrically irregular as sketched in Figure 3:

FE discretization in space needs highly irregular tetrahedra, which, combined with the high heterogeneity of the medium, gives raise to an ill-conditioned discretized diffusion matrix. After FE discretization a system of ODEs is left:

$$\begin{array}{ccc}\hfill B\frac{\partial \mathbf{u}(t)}{\partial t}& =& -A\mathbf{u}(t)+\mathbf{b}(t).\hfill \end{array}$$

where A is the SPD stiffness matrix, B is the SPD mass matrix. The right-hand-side account for source terms and Neumann boundary conditions. Using a time marching scheme (e.g., Crank-Nicolson for stiff problems) yields a sequence of linear systems

$$(A+{\sigma }_{k}B){\mathbf{u}}_k={\mathbf{b}}_k,{\sigma }_k=\frac{2}{\Delta t_k},k=1,Nstep.$$

(12)

The discretized FE matrices have both n = 589,661 and a number of nonzeros $nz$ = 8,833,795. We now report the results of a complete simulation which took $Nstep=32$ steps. The near steady state is reached after T = 36,000 days. Initially we set $\Delta t_1=60$ (days). Then we incremented the timesteps as

$$\Delta t_k=min(1.2\Delta t_{k-1},7300,T-\Delta t_{k-1})⟹{\sigma }_k={\displaystyle \frac{2}{\Delta t_k}},k=2,\dots ,Nstep.$$

This formula allows the timestep size to constantly increase and to take its largest value approaching the steady state solution. The initial preconditioner $P_0$ was evaluated employing the Matlab function ichol with

(1)

droptol$={10}^{-5}$ applied to $A_1=A+5B$ to enhance diagonal dominance; this resulted in a triangular Cholesky factor with density $5.18$ and

(2)

droptol$={10}^{-4}$ applied to $A_2=A+20B$ (density $2.83$).

We solved the first linear system to a high accuracy to assess as accurately as possible the $p=10$ leftmost eigenvectors. Using the residual ${\delta }_j={\displaystyle \frac{\parallel P_{0}A{\mathbf{w}}_j-{\mu }_j{\mathbf{w}}_j\parallel }{\parallel {\mathbf{w}}_j\parallel }}$ to measure the accuracy of the approximate eigenvectors we found:

• Case (1). 179 PCG iterations for the first linear system; ${\delta }_j\in [1.1{10}^{-3},5.8{10}^{-3}]$.
• Case (2). 370 PCG iterations for the first linear system; ${\delta }_j\in [1.1{10}^{-3},3.3{10}^{-3}]$.

This set of vectors is then used to construct deflation, tuned and spectral preconditioners for all the subsequent linear systems in the sequence (with exit test on relative residual and tolerance tol $=1e-10$).

The results are summarized in

Table 6. Number of iterations and total CPU time for solving the sequence (12) of shifted linear system. $†=$ for some of the system in the sequence convergence not achieved.

	IC (${\mathit{A}}_1,1\mathit{e}-5$)			IC (${\mathit{A}}_2,1\mathit{e}-4$)
update	# its	CPU	CPU per it.	# its	CPU	CPU per it.
no update	8231	805.8	0.098	16582	1177.2	0.071
spectral	5213	557.3	0.107	10225	817.5	0.080
SR1 tuned	5161	551.7	0.107	10152	811.4	0.080
BFGS tuned	5178	612.3	0.118	10198	924.6	0.091
deflated	†	†	†	†	†	†

. The eigenvectors estimated during the first linear system are adequate to accelerate the subsequent iterations for both the tuned and spectral approaches which provide an improvement in both iteration number and CPU time. Note also that the time per iteration is only slightly increased with respect to the “no update” case.

Instead, the PCG residuals obtained with the deflation preconditioner in some instances stagnate before reaching the exit test. The lack of robustness of the deflation approach is also observed in Reference [5] where the authors suggest a reorthogonalization of the residuals at the price of increasing costs per iteration and also Reference [26] in the framework of iterative eigensolvers. Regarding tuned preconditioners, SR1 is definitely superior to the BFGS approach due to its lower computational cost.

In Figure 4, left, we plot the number of iterations needed by the “no update” and tuned preconditioners for every system in the sequence compared also with the scaled $\sqrt{\Delta t_k}\approx {\sigma }_k^{-1/2}$ values (blue lines) as indicators of the increasing condition number of the coefficient matrices. On the right plot the behavior of the tuned vs spectral preconditioners is displayed, showing the substantial equivalence between the two approaches. Note finally that the PCG solver with the tuned preconditioner is almost insensitive to the increasing condition number of the matrices, displaying a number of iterations roughly constant throughout the timesteps.

7.2. Iterative Eigensolvers

When computing the smallest or a number of the interior eigenvalues of a large and sparse matrix, most iterative solvers require the solution of a shifted linear system of the type

$$(A-\sigma B)\mathbf{x}=\mathbf{b}.$$

This is true for instance in the shift-invert Arnoldi method, in the inverse power iteration, in the Rayleigh Quotient iteration. A related linear system is also to be solved within the correction equation in the Jacobi-Davidson method. Other gradient-based methods such as the Locally Optimal Block PCG method (LOBPCG) [46] or the Deflation-Accelerated Conjugate Gradient method (DACG) [47] implicitly solve a shifted linear system. The sequence of linear systems are characterized by a constant or a slowly varying parameter $\sigma$ so informations on matrix A or on the pencil $(A,B)$ can be profitably used for all the sequence.

We denote as $0<{\lambda }_1\le {\lambda }_2\le \dots {\lambda }_m\dots \le {\lambda }_n$ its eigenvalues and ${\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_m,\dots {\mathbf{v}}_n$ the corresponding eigenvectors. Our aim is to investigate the efficiency of the SR1 tuned preconditioner in the acceleration of the so-called correction equation in the simplified Jacobi-Davidson (JD) method [48,49].

To assess eigenpair $({\lambda }_j,{\mathbf{v}}_j)$ a number of linear systems have to be solved of the form

$$\begin{array}{ccc}\hfill J_k^{(j)}\mathbf{s}& =& -(A-{\theta }_j^{(k)}I){\mathbf{u}}_k;\mathrm{for}\mathbf{s}\perp {\mathbf{u}}_k,\mathrm{where}\hfill \end{array}$$

(13)

$$\begin{array}{ccc}\hfill J_k^{(j)}& =& (I-{\mathit{QQ}}^T)(A-{\theta }_j^{(k)}I)(I-{\mathit{QQ}}^T),{\theta }_j^{(k)}=q({\mathbf{u}}_k)\equiv \frac{{\mathbf{u}}_k^{T}A{\mathbf{u}}_k}{{\mathbf{u}}_k^T{\mathbf{u}}_k},\hfill \\ \hfill Q& =& \left[\begin{array}{ccc}{\mathbf{v}}_1& \dots {\mathbf{v}}_{j-1}& {\mathbf{u}}_k\end{array}\right].\hfill \end{array}$$

(14)

The PCG method is proved to converge if applied to systems (13) as the Jacobian $J_k^{(j)}$ is shown to be SPD in the subspace orthogonal to ${\mathbf{u}}_k$ [22,49]. Following the works in [27,50], which are based on a two-stage algorithm, the columns of W are orthonormal approximate eigenvectors of A provided by a gradient method, DACG [47], satisfying

$$\begin{array}{ccc}\hfill A{\tilde{\mathbf{v}}}_s={\tilde{\lambda }}_s{\tilde{\mathbf{v}}}_s+{\mathbf{res}}_s,\parallel {\mathbf{res}}_s\parallel & \le & \tau {\lambda }_s\hfill \\ s=1,\dots ,m.\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill {\tilde{\lambda }}_s-{\lambda }_s& \le & \tau {\lambda }_s.\hfill \end{array}$$

(16)

These vectors are also used as initial guesses for the Newton phase. To update a given initial preconditioner for $({\lambda }_j,{\mathbf{v}}_j)$ we use the next approximate eigenvectors ${\mathbf{v}}_{j+1},\dots ,{\mathbf{v}}_m$ to define an SR1 tuned preconditioner as

$$P_j=P_0-Z{\left(Z^{T}A{W}_j\right)}^{-1}{Z}^T,Z=P_{0}A{W}_j-W_j,\mathrm{with}{W}_j=\left[\begin{array}{ccc}{\tilde{\mathbf{v}}}_{j+1}& \dots & {\tilde{\mathbf{v}}}_m\end{array}\right].$$

(17)

We recall the following result stated in ([27], Lemma 3.1):

Lemma 1.

Let $P_j$ a block tuned preconditioner satisfying condition (17). In the hypothesis (15) each column of $W_j$ that is, ${\tilde{\mathbf{v}}}_s,s=j+1,\dots ,m$, is an approximate eigenvector of $P_j(A-{\theta }_{j}I)$ satisfying

$$P_j(A-{\theta }_{j}I){\tilde{\mathbf{v}}}_s=\left(1-\frac{{\theta }_j}{{\tilde{\lambda }}_s}\right){\tilde{\mathbf{v}}}_s+{\epsilon }_s,with\parallel {\epsilon }_s\parallel \le \tau {\lambda }_{j+1}\parallel P_j\parallel .$$

The effect of applying a tuned preconditioner to $A-{\theta }_{j}I$ is to set a number of eigenvalues of $P(A-{\theta }_{j}I)$ to a value that is close to one only under the conditions that the eigenvalues are well separated, that is, ${\displaystyle \frac{{\lambda }_j}{{\lambda }_{j+1}}}\ll 1$, which is not always the case in realistic problems.

In order to define a more effective preconditioner for shifted linear systems one can allow the preconditioned matrix $PA$ to have eigenvalues different from one corresponding to the columns of matrix W. In Reference [50] a generalized block tuned (GBT) preconditioner is defined:

Definition 3.

Given a preconditioner $P_0$, an $n\times p$ matrix W with full column rank, and a diagonal matrix $\Gamma =\mathrm{diag}({\gamma }_1,\dots ,{\gamma }_p)$, a generalized block tuned preconditioner for matrix A is a matrix P obtained by correcting $P_0$ with a low-rank matrix depending on $A,W$ and Γ and satisfying

$$PAW=W\Gamma .$$

(18)

The generalized SR1 preconditioner is defined as

$$P=P_0-Z{\Pi }^{-1}{Z}^T,\mathrm{where}Z=P_{0}AW-W\Gamma ,\mathrm{and}\Pi =Z^{T}AW.$$

(19)

Note that the above preconditioner is not in general symmetric as small matrix $\Pi$ is not and hence its use would prevent convergence either of the DACG eigensolver or the inner PCG iteration within the Newton method. However, this drawback can be circumvented when $W\equiv W_j$ represents the matrix of the (approximate) eigenvectors corresponding to the diagonal matrix with the eigenvalues of A: ${\Lambda }_j=\mathrm{diag}({\tilde{\lambda }}_{j+1},{\tilde{\lambda }}_{j+2},$$\dots ,{\tilde{\lambda }}_m)$. In such case we can approximate $\Pi$ as

$$\Pi =W_j^{T}APA{W}_j-\Gamma W_j^{T}A{W}_j\approx W_j^{T}APA{W}_j-\Gamma {\Lambda }_j=\tilde{\Pi },$$

(20)

so restoring symmetry. This modified preconditioner ${\tilde{P}}_j=P_0-Z{\tilde{\Pi }}^{-1}Z$ satisfies only approximately the tuning property as

$${\tilde{P}}_{j}A{\tilde{\mathbf{v}}}_s={\gamma }_s{\tilde{\mathbf{v}}}_s+{\mathcal{E}}_s,\parallel {\mathcal{E}}_s\parallel \le \tau {\lambda }_s\parallel Z{\tilde{\Pi }}^{-1}\Gamma \parallel ,s=j+1,\dots ,m.$$

(21)

The following theorem states that it is however possible to have p eigenvalues of the preconditioned matrix ${\tilde{P}}_j(A-{\theta }_{j}I)$ very close to one depending on how the columns of matrix W approximate the eigenvectors of A. We also define the reciprocal of the relative separation between pairs of eigenvalues as

$${\xi }_j={\displaystyle \frac{{\lambda }_{j+1}}{{\lambda }_{j+1}-{\lambda }_j}}.$$

(22)

Theorem 3.

Let matrix $W_j$ be as in (17), ${\tilde{P}}_j$ an approximate GBT preconditioner satisfying condition (21), with ${\gamma }_s={\tilde{\lambda }}_s/({\tilde{\lambda }}_s-{\tilde{\lambda }}_j),s=j+1,\dots ,m$, then each column of $W_j$ is an approximate eigenvector of ${\tilde{P}}_j(A-{\theta }_{j}I)$ corresponding to the approximate eigenvalue

$${\mu }_s={\displaystyle \frac{{\tilde{\lambda }}_s-{\theta }_j}{{\tilde{\lambda }}_s-{\tilde{\lambda }}_j}}.$$

satisfying

$${\tilde{P}}_j(A-{\theta }_{j}I){\tilde{\mathbf{v}}}_s={\mu }_s{\tilde{\mathbf{v}}}_s+{\epsilon }_s,with\parallel {\epsilon }_s\parallel \le \tau \left({\lambda }_s\parallel Z{\tilde{\Pi }}^{-1}\Gamma \parallel +{\lambda }_{j+1}\parallel {\tilde{P}}_j\parallel \right).$$

Proof. 

See Reference [50]. □

The eigenvalues ${\mu }_s$ can be very close to one depending on the initial tolerance $\tau$ as stated in Corollary 1, under reasonable additional hypotheses

Corollary 1.

Let ${\theta }_j,\tau$ such that, for all $j=1,\dots ,m-1$:

$$\begin{array}{ccc}\hfill {\lambda }_j& \le & {\theta }_j\le {\tilde{\lambda }}_j,\hfill \end{array}$$

(23)

$$\begin{array}{ccc}\hfill \tau & <& {(2{\xi }_j)}^{-1},\hfill \end{array}$$

(24)

then

$$1\le {\mu }_s\le 1+2\tau ({\xi }_j-1),s=j+1,\dots ,m.$$

Proof. 

See Reference [50]. □

From Corollary 1 it is clear that ${\mu }_s$ can be made arbitrarily close to one by appropriately reducing the tolerance $\tau$. As an example, if ${\xi }_j={10}^2$, $\tau ={10}^{-3}$ then all ${\mu }_s$ are expected to be in $(1,1.2).$

In Reference [50] (Theorem 2) is also stated that the eigenvalues of the preconditioned projected Jacobian matrix (13) are characterized in the same way as stated in Theorem 3, that is, that for a suitable function $C(\tau )$, increasing in $\tau$

$$P_Q{J}_k^{(j)}{\tilde{\mathbf{v}}}_s={\mu }_s{\tilde{\mathbf{v}}}_s+\mathbf{err},\parallel \mathbf{err}\parallel \le \tau C,$$

(25)

where ${\tilde{P}}_j$ a generalized block tuned preconditioner, $P_Q=(I-{\mathit{QQ}}^T){\tilde{P}}_j(I-{\mathit{QQ}}^T)$.

To limit the cost of the application of the low-rank correction it is customary to fix the maximum number of columns of matrix $W_j$, parameter $l_{max}$. Conversely, it is required to enlarge it when assessing eigenpairs with $j\approx m$. The first DACG stage is hence used to compute an extra number $(\mathrm{win})$ of approximated eigenpairs. In this way the number of columns of $W_j$ is $m_j=min\{l_{max},m+\mathrm{win}-j\}$.

The construction (C) of ${\tilde{P}}_j$ and its application (A) as ${\tilde{P}}_j\mathbf{r}$ are sketched below. MVP = matrix-vector products, $Z_j=Z_0($:,j+1:j+$m_j$) and ${\Pi }_j={\Pi }_0($j+1:j+$m_j$, j+1:j+$m_j$).

Phase	When	What	Relevant Cost
C	once and for all	$Z_0=P_{0}A{V}_0$ ${\Pi }_0=Z_0^{T}A{V}_0$	$m+\mathrm{win}$ MVPs and applications of $P_0$ ${(m+\mathrm{win})}^2/2$ dot products.
C	for every eigenpair	$Z=Z_j-V_j\Gamma$ $\tilde{\Pi }={\Pi }_j-\Gamma {\Lambda }_j$	$m_j$ daxpys
A	at each iteration	$\mathbf{h}=Z^T\mathbf{r}$ $\mathbf{g}=\tilde{\Pi }\backslash \mathbf{h}$ $\mathbf{w}=P_0\mathbf{r}-Z\mathbf{g}$	$m_j$ dot products 1 system solve of size $m_j$ 1 application of $P_0$, $m_j$ daxpys

We report from Reference [50] the results of the described preconditioner in the computation of the 20 smallest eigenpairs of matrix thermomec, which is available in the SuiteSparse Matrix Collection at https://sparse.tamu.edu/. This is a challenging test due to the high clustering of its smallest eigenvalues. The CPU times (in seconds) refer to running a Fortran 90 code on a 2 x Intel Xeon CPU E5645 at 2.40 GHz (six core) and with 4 GB RAM for each core. The exit test is on the relative eigenresidual: ${\displaystyle \frac{\parallel A\mathbf{u}-q(\mathbf{u})\mathbf{u}\parallel }{q(\mathbf{u})}\le \epsilon ={10}^{-8}}$, with $q(\mathbf{u})$ defined in (14).

The results of the runs are displayed in

Table 7. Timings and iterations for the DACG -Newton method for the computation of $m=20$ eigenpairs of matrix thermomec. In boldface the smallest overall number of iterations and CPU time.

			DACG			Newton			Total
						Iterations		CPU
Prec.	win	${\mathit{l}}_{max}$	$\mathit{\tau }$	Its.	CPU	OUT	Inner		MVP	CPU
Fixed	0	0	${10}^{-4}$	1510	15.86	153	2628	34.12	4291	52.97
SR1 Tuned	10	10	${10}^{-3}$	1335	14.89	137	2187	32.19	3659	51.48
GBT	5	10	${10}^{-3}$	777	11.16	44	607	9.42	1428	20.74

where the number of iterations (and CPU time) of the initial DACG method are reported as well as the cumulative number of iterations (and CPU time) needed for solving all the linear systems (13) within the Newton iteration for the 20 eigenpairs. Especially the second (Newton) phase takes great advantage by the GBT preconditioner as also accounted for by Figure 5 shows the number of iterations taken by the Newton phase with the fixed IC, SR1 tuned and GBT preconditioners, together the (scaled) $log{\xi }_j$, which is a measure of the condition number of the inner linear systems.

The almost constant GBT curve confirms the property of this preconditioner which makes the number of iterations nearly independent on the relative separation between consecutive eigenvalues.

8. Conclusions

We have presented a general framework of updating a given preconditioner $P_0$ with a low-rank matrix with the aim of further clustering eigenvalues of the preconditioned matrix $P_{0}A$. We have shown that this approach is particularly efficient when a sequence of linear systems has to be solved. In this case the cost of computation of the leftmost eigenvector of the preconditioned matrices is payed for by the number of system solutions. Alternatively, the vector forming the low-rank updates can be chosen as the previous solutions in the sequence. In summary we have reviewed three updating processes:

• Deflation, aimed at shifting to zero a number of approximate eigenvalues of $P_{0}A$.
• Tuning, aimed at shifting to one a number of approximate eigenvalues of $P_{0}A$.
• Spectral preconditioners, aimed at adding one to a number of approximate eigenvalues of $P_{0}A$.

The most popular choices of the vectors (such as leftmost eigenvectors of either A or $P_{0}A$) have been considered together with some techniques to cheaply assess them.

Finally we have considered a number of applications, such as discretization of evolutionary PDEs and solution of large eigenvalue problems. In all these applications the low-rank correction is much beneficial to reduce the number of iterations of the iterative solver of choice.

As a general recommendation we can say that accelerating a given preconditioner by a low-rank matrix can be useful when both the following situations occur:

• A sequence of linear systems with constant or slightly varying matrices has to be solved.
• Either the smallest or the largest eigenvalues do not form a cluster.