Low-rank Updates of Preconditioners for sequences of symmetric linear systems

The aim of this survey is to review some recent developements in devising efficient preconditioners for sequences of linear systems A x = b. Such a problem arise in many scientific applications, such as discretization of transient PDEs, solution of eigenvalue problems, (Inexact) Newton method applied to nonlinear systems, rational Krylov methods for computing a function of a matrix. Full purpose preconditioners such as the Incomplete Cholesky (IC) factorization or approximate inverses 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 P0 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) 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 P0 A. 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.


Introduction
The aim of this survey is to review some recent developements in devising efficient preconditioners for sequences of linear systems Ax = b. Such a problem arise in many scientific applications, such as discretization of transient PDEs, solution of eigenvalue problems, (Inexact) Newton method applied to nonlinear systems, rational Krylov methods for computing a function of a matrix. Full purpose preconditioners such as the Incomplete Cholesky (IC) factorization or approximate inverses 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) 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. 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.

Synopsis
In Section 1 we will discuss various low-rank updates, namely, the deflation technique, the tuning strategy, developed based on a revisited secant condition. and the spectral preconditioners. In Section 2 we will concerned with the computational complexity of the low-rank updates. Section 3 discusses the choice of the vectors which form the low-rank matrices while Section 4 will address the issue to cheaply assessing some of the leftmost eigenpairs of the coefficient matrix. Section 5 is devoted to reporting some numerical experiments in the framework of discretization of transient PDES and solution eigenvalue problems. The conclusions are drawn in Section 6.

Deflation
The idea to use deflation to accelerate the Conjugate Gradient method goes back to the middle of the 80's. The most cited article on this subject is a paper by Nicolaides [1] who developed a deflated CG based on a three-term recurrence. Other similar formulations are present in [2] and [3]. Some years later in [4] a deflated GMRES method has been proposed by using an augmented basis. Rigorous bounds on the eigenvalues of the deflated-preconditioned matrices are given in [5]. We finally mention recent interesting survey on deflation methods in [6].
We will follow here the formulation of the Deflated-CG algorithm developed in [7]. Given a full column rank n × p matrix W = w 1 . . . w p the A-orthogonal projector H is defined as which maintains the CG residuals orthogonal to the subspace spanned by {w 1 , . . . , w p }, provided that also r 0 satisfies W T r 0 = 0. The Deflated-PCG algorithm is described below.
Note that this algorithm is mathematically equivalent to the PCG method with preconditioner P = HP 0 H T since r T k z k = r T k P 0 r k = r T k HP 0 H T r k , being W T r k = 0 and therefore H T r k = r k . P is only symmetric positive semidefinite, however PCG could not breakdown as stated in the following result [7,Theorem 4.2].
Theorem 1. Let A be an SPD matrix and W an n × p matrix with full column rank. The Deflated-PCG algorithm applied to the linear system Ax = b will not break down at any step. The approximate solution x (k) is the unique minimizer of x (k) − x ( * ) A over the affine space x 0 + span(W, K k (r 0 )).
The action of this deflated preconditioner is to put some of the eigenvalues of the preconditioner matrix to zero in fact: i.e. 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 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.

Tuning
The adjective tuned associated with a preconditioner has been introduced in [8] 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 w, a tuned preconditioner for matrix A is a matrix P ≡ M −1 obtained by correcting P 0 by a rank one or rank two matrix depending on A and w and satisfying In this way matrix M agrees with A in the direction w and also the preconditioned matrix has the eigenvalue 1 corresponding to the eigenvector w. This definition can be easily extended to multiple vectors: Definition 2. Given a preconditioner P 0 and an n × 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 In [9] the following single-vector tuned preconditioned for SPD linear systems was proposed (here M 0 = LL T ≈ A and P = M −1 is computed by the Sherman-Morrison formula).
It is easily checked that PAw = P 0 Aw − z = w.

Digression
This tuned preconditioner, and many others, has however an older derivation. Let us now consider the solution of a nonlinear system of equations Quasi-Newton methods construct a sequence of approximations of the Jacobians B k ≈ F (x k ). 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) [10][11][12][13][14] Among those Quasi-Newton methods we recall three of the most commonly used (note that Broyden's method: All these updates satisfy the so called secant condition which reads 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 −1 0 . We can use all the preceding equations after employing the following substitutions: which transform the secant condition into that is the tuning property. We can then develop three different tuning strategies. For each of them we write the "direct" formula (with M, M 0 ) the inverse formulation (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 x ∈ R n we set z = H T x and u = W T x. Then we have the last inequality holding since both Π −1 and P 0 are SPD matrices. The inequality is strict since if u = 0 then W T x = 0 and hence z = ( 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 [15,Sec. 3.2]. which has however a modest practical use as (W T AW) −1 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Θ, with Θ = diag µ 1 , . . . , µ p the diagonal matrix with the eigenvalues of the initially preconditioned matrix. Assume further that µ j < 1, j = 1, . . . , p. Since P 0 A is not symmetric but P 1/2 0 AP 1/2 0 is so it follows that W must satisfy W T P −1 0 W = I and hence W T AW = Θ (see the discussion on the forthcoming Section 5.1). In such a case we have that matrix is obviously SPD and therefore so is P which is the sum of an SPD matrix (P 0 ) and a symmetric positive The Broyden tuned preconditioner has been used in [16] to accelerate the GMRES for the inner linear systems within the Arnoldi method for eigenvalue computation. The SR1 preconditioner appeared in the cited works [9,15] and also in [17]. This low-rank update has also been employed to accelerate the PCG method in the solution of linear systems involving SPD Jacobian matrices. In [18] some conditions are proved under which the SR1 update maintains the symmetric positive definiteness of a given initial preconditioner. The BFGS preconditioner has been used (under the name of balancing preconditioner) in [19], in [20] for eigenvalue computation and also in [21] 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 [22].
We finally mention the work in [23] where the author use this correction to update an Inexact Constraint Preconditioner, by damping the largest eigenvalues of the preconditioned matrices.

Spectral preconditioners
A further strategy to improve the preconditioner's efficiency by a low-rank update can be found in [17,[24][25][26]. A spectral preconditioner for SPD linear systems is defined as with the leftmost eigenvectors of P 0 A as the columns of W. Denoted as µ 1 , µ 2 , . . . , µ p the smallest eigenvalues of P 0 A it is easily checked that matrix PA has the same eigenvectors w 1 , . . . , w p as P 0 A with µ 1 + 1, . . . , µ p + 1 as eigenvalues.

Implementation and computational complexity
We sketch below the most time-consuming computational tasks when using a low-rank update of a given preconditioner. We first evaluate the preliminary cost, which has to be done before starting of the iterative process: 1. All low-rank updates: Computation of A W = AW (O(np) operations). 2. Broyden and SR1 only: Computation of Z = P 0 A W − W (p applications of the initial preconditioner).

All low-rank updates: Computation of
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 the same, being based on the following main tasks, in addition to the application of the initial preconditioner P 0 to a given vector.

Multiplication of a p × n matrix times a vector (O(np) operations)
The BFGS preconditioner, as it corrects P 0 by a rank-(2p) matrix, roughly doubles tasks 1-3 and hence the total updating cost at each iteration.

Choice of the vectors {w j }.
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 i.e. satisfying P 0 Aw j − µ j w j ≈ δ. The coefficient matrix has been chosen as the FD discretization of the Laplacian on an L-shaped domain obtained using the MatLAB command A = d e l s q ( numgrid ( ' L ' , 5 0 0 ) ) ; which returns a sparse matrix of order n = 186003. The linear system Ax = b, with b a random uniformly distributed vector, has been solved by the PCG method with various low-rank update techniques with P 0 = IC(0).
The results in Table 4 show that a very low accuracy on the eigenvectors (absolute residual of order of δ = 0.01) is sufficient to provide good preconditioning results. Also notice that the tuned preconditioner seems to be more sensitive to the parameter δ. 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 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 Aw j − λ j w j ≈ δ. 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).

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 w j is to use some of the solutions to the previous linear systems in the sequence. Assume now that the matrices in the sequence are allowed to (slightly) change, i.e. we want to solve In this case computing 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 (j > 1), the j e f f = min{j − 1, p} solutions to the previous linear systems as columns of W ≡ W j , thus yielding: The idea beyond this choice, which is completely cost-free, comes from the reasoning that a solution x k of the k-th linear system has with high probability non-negligible components in the weights provided by the smallest eigenvalues.
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 [26] which gives raise to a sequence of (2 × 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 linear systems of the form where A is SPD, B is rectangular with full row rank, D is diagonal and possibly indefinite. We take 30 consecutive linear systems from this sequence, and solve them with the MINRES iterative method using as the initial preconditioner the incomplete Cholesky factorization of A with drop tolerance 10 −2 and the following updates 1. BFGS, with update directions as in (10).
3. BFGS, with as update directions the accurate leftmost eigenvectors of A k .
The results are reported in Figure 2 where we can appreciate the number of iterations reduction provided by the low-rank update. The (cheaper) spectral preconditioner provides as good acceleration as the BFGS approach with update directions as in (10). Moreover the number of iterations is comparable to those obtained with the "exact" eigenvectors (which are obtained by the Matlab eigs function).

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 [27,28]. 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 : Matrix V m is such that V T m P −1 0 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 α k , β k : Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 March 2020 doi:10.20944/preprints202003.0094.v1 Matrices V m and T m obey to the classical Lanczos relation i.e.: After eigensolving T m , thus obtaining T m = QΛ m Q T , the eigenvectors corresponding to the p smallest eigenvalues are selected as columns of matrix Q p . Finally, the columns of W p = V m Q p are expected to approximate p of the 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 [29] the authors suggest to implement 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 have been can be found in the literature. In [7] 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 p k instead of the preconditioned residuals (they actually lie in the same Krylov subspace).

Numerical results
The technique previously describe 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 [29]. However, in most of the situations this high accuracy is not really needed (see the discussion in Section 6.5).
In this section we will consider a particular case of sequences of linear systems with (slightly) changing matrices as which will be denoted as shifted linear systems. These sequences appear e.g. 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.

FE Discretization of a parabolic equation
Consider the parabolic PDE 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 [30]. 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: 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 The discretized FE matrices have both n = 589661 and a nonzero number nz = 8833795. We now report the results of a complete simulation which took Nstep = 32 steps. The near steady state is reached after T = 36000 days. Initially we set ∆t 1 = 60 (days). Then we incremented the timesteps as The initial preconditioner P 0 has been 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 δ j = P 0 Aw j − µ j w j w j to measure the accuracy of the approximate eigenvectors we found: This set of vectors is used them 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. 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 [7] where the authors suggest a reorthogonalization of the residuals at the price of increasing costs per Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 March 2020 doi:10.20944/preprints202003.0094.v1 iteration and also in [16] in the framework of iterative eigensolvers. Regarding tuned preconditioners, SR1 is definitely superior to the BFGS approach due do 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 √ ∆t k ≈ σ −1/2 k values 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.

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 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 LOBPCG [31] or DACG [32] implicitly solve a shifted linear system. The sequence of linear systems are characterized by a constant or a slowly varying parameter σ so informations on matrix A or on the pencil (A, B) can be profitably used for all the sequence.
We denote as 0 < λ 1 ≤ λ 2 ≤ . . . λ m . . . ≤ λ n its eigenvalues and v 1 , v 2 , . . . , v m , . . . 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 [33,34]. To assess eigenpair (λ j , v j ) a number of linear systems have to be solved of the form The PCG method is proved to converge if applied to systems (12) as the Jacobian J (j) k is shown to be SPD in the subspace orthogonal to u k [12,34]. Following the works in [17,35], which are based on a two-stage algorithm, the columns of W are orthonormal approximate eigenvectors of A provided by a gradient method, DACG [32], satisfying These vectors are also used as initial guesses for the Newton phase. To update a given initial preconditioner for (λ j , v j ) we use the next approximate eigenvectors v j+1 , . . . , v m to define an SR1 tuned preconditioner as We recall the following result stated in [ Lemma 1. Let P j a block tuned preconditioner satisfying condition (15). In the hypothesis (13) each column of W j i.e.ṽ s , s = j + 1, . . . , m, is an approximate eigenvector of P j (A − θ j I) satisfying The effect of applying a tuned preconditioner to A − θ j I is to set a number of eigenvalues of P(A − θ j I) to a value that is close to one only under the conditions that the eigenvalues are well separated, i.e. λ j λ j+1 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 [35] a generalized block tuned (GBT) preconditioner is defined: Definition 3. Given a preconditioner P 0 , an n × p matrix W with full column rank, and a diagonal matrix Γ = diag(γ 1 , . . . , γ 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 The generalized SR1 preconditioner is defined as Note that the above preconditioner is not in general symmetric as small matrix Π 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 ≡ W j represents the matrix of the (approximate) eigenvectors corresponding to the diagonal matrix with the eigenvalues of A: Λ j = diag(λ j+1 ,λ j+2 , . . . ,λ m ). In such case we can approximate Π as so restoring symmetry. This modified preconditioner P j = P 0 − ZΠ −1 Z T satisfies only approximately the tuning property as The following theorem states that it is however possible to have p eigenvalues of the preconditioned matrix P j (A − θ 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 Theorem 3. Let matrix W j be as in (15), P j an approximate GBT preconditioner satisfying condition (19), with γ s =λ s /(λ s −λ j ), s = j + 1, . . . , m, then each column of W j is an approximate eigenvector of P j (A − θ j I) corresponding to the approximate eigenvalue satisfying P j (A − θ j I)ṽ s = µ sṽs + ε s , with ε s ≤ τ λ s ZΠ −1 Γ + λ j+1 P j .
The eigenvalues µ s can be very close to one depending on the initial tolerance τ as stated in Corollary 1, under reasonable additional hypotheses Corollary 1. Let θ j , τ such that, for all j = 1, m: Proof. See [35].
From Corollary 1 it is clear that µ s can be made arbitrarily close to one by appropriately reducing the tolerance τ. As an example, if ξ j = 10 2 , τ = 10 −3 then all µ s are expected to be in (1, 1.2).
In [35,Theorem 2] is also stated that the eigenvalues of the preconditioned projected Jacobian matrix (12) are characterized in the same way as stated in Theorem 3, i.e. that for a suitable function C(τ), increasing in τ where P j a generalized block tuned preconditioner, P Q = (I − QQ T )P j (I − 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 ≈ m. The first DACG stage is hence used to compute an extra number (win) of approximated eigenpairs. In this way the number of columns of W j is m j = min{l max , m + win − j}.
The construction (C) of P j and its application (A) as P j r are sketched below. MVP = matrix-vector products, Z j = Z 0 (:,j+1,j+m j ) and Π j = Π 0 (j+1:j+m j , j+1:j+m j ). phase when what relevant cost C once and • Z 0 = P 0 AV 0 m + win MVPs and applications of P 0 for all We report from [35] 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.40GHz (six core) and with 4GB RAM for each core. The exit test is on the relative eigenresidual: Au − q(u)u q(u) ≤ ε = 10 −8 .
The results of the runs are displayed in Table 7 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 (12) within the Newton iteration for the 20 eigenpairs.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 March 2020 doi:10.20944/preprints202003.0094.v1 Especially the second (Newton) phase takes great advantage by the GBT preconditioner as also accounted for by Figure 5 which shows the number of iterations taken by the Newton phase with the fixed IC, SR1 tuned and GBT preconditioners. 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.

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. Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 5 March 2020 doi:10.20944/preprints202003.0094.v1 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: 1. A sequence of linear systems with constant or slightly varying matrices has to be solved.
2. Either the smallest or the largest eigenvalues do not form a cluster.