Hussam Al Daas

Laura Grigori

Pascal Hénon

Philippe Ricoux

Recyclage de Sous-Espaces de Krylov et Troncature de Sous-Espaces de Déflation pour Résoudre Séquence de Systèmes Linéaires

Keywords: espaces de Krylov, solveur linéaires itératifs, recyclage, déflation, réduction, svd Krylov subspaces, iterative linear solvers, recyling, deflation, reduction, svd

Ce papier présente plusieures stratégies de déflation liées aux méthodes de recyclage de sous-espaces de Krylov pour résoudre une séquence de systèmes linéaires. À côté de stratégies de déflation très connues qui sont basées sur la déflation des vecteurs de Ritz et Ritz harmonique, on introduit une technique de déflation basée sur la décomposition en valeurs singulières. On considère deux contextes du recyclage, le recyclage de l'espace de Krylov entre les cycles de resart et le recylcage de l'espaces de déflation quand la matrice change dans la séquence. L'efficacité de la méthode proposée est étudiée sur des séquence de systèmes linéaires issues de la modélisation de réservoirs.

Introduction

We consider the system Ax = b

where A ∈ C n×n is non-singular, b ∈ C n is a right-hand side vector and x ∈ C n is the vector of unknowns. We are interested in solving (1) by using a Krylov iterative method. GMRES [START_REF] Saad | GMRES : A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF], CG [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF], GCR [START_REF] Eisenstat | Variational iterative methods for nonsymmetric systems of linear equations[END_REF], and BiCG [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] are widely used Krylov subspace methods. Several works have studied the impact of the spectrum of the matrix A on the convergence rate of Krylov methods [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF][START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Morgan | GMRES with deflated restarting[END_REF][START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF][START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF][START_REF] Sacchi | A GMRES convergence analysis for localized invariant subspace ill-conditioning[END_REF]. It was noticed that the small and nearly isolated eigenvalues induce a slow rate of convergence. In order to have robustness and fast convergence, a deflation subspace correction can annihilate the impact of this part of the spectrum. Different types of deflation techniques depending on the Hermitian structure of the matrix A are described in the literature [START_REF] Nicolaides | Deflation of conjugate gradients with applications to boundary value problems[END_REF][START_REF] Frank | On the construction of deflation-based preconditioners[END_REF][START_REF] Yogi | Deflation and balancing preconditioners for krylov subspace methods applied to nonsymmetric matrices[END_REF][START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF][START_REF] Saad | h. A deflated version of the conjugate gradient algorithm[END_REF][START_REF] Morgan | GMRES with deflated restarting[END_REF][START_REF] Erhel | Restarted GMRES preconditioned by deflation[END_REF][START_REF] Morgan | A restarted GMRES method augmented with eigenvectors[END_REF][START_REF] André | A framework for deflated and augmented krylov subspace methods[END_REF]. Some of these strategies construct a deflation subspace during the iterative solve of [START_REF] Kapil | Recycling bicgstab with an application to parametric model order reduction[END_REF].

In several applications, we have to solve successive linear systems A i x = b i where A i ≈ A i+1 , e.g., large-scale eigenvalue problems and Newton's method for solving nonlinear systems. Furthermore, Newton's method combined with a Constrained Pressure Residual (CPR [START_REF] Wallis | Constrained residual acceleration of conjugate residual methods[END_REF][START_REF] Cao | Parallel scalable unstructured CPR-Type Linear Solver for Resevoir Simulation[END_REF]) solver have both cases, A i ≈ A i+1 and A i = A i+1 . Thus, the deflation subspace that is built during the solution of A i x = b i can be useful to solve A i+1 x = b i+1 . However, the deflation subspaces introduced in the previously mentioned work may not be directly used to solve successive linear systems with different right-hand sides even if the matrix does not change.

In [START_REF] Rey | A rayleigh-ritz preconditioner for the iterative solution to large scale nonlinear problems[END_REF], the authors introduced a recycling strategy based on Rayleigh-Ritz preconditioner to solve sequence of linear systems of equations with symmetric positive definite matrices. In the same time, De Sturler presented in [START_REF] Sturler | Truncation strategies for optimal krylov subspace methods[END_REF] a truncation strategy based on orthogonality constraints to select vectors from the search subspace and used them to accelerate convergence after restarting the method or when solving a new linear system in a sequence. Later, Parks et al. introduced the method GCRO-DR in [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF]. This method is a variant of the restarted GMRES in which an abstract subspace can be used as a deflation subspace. Several papers followed based on this work [START_REF] Carvalho | A flexible generalized conjugate residual method with inner orthogonalization and deflated restarting[END_REF][START_REF] Gosselet | Total and selective reuse of krylov subspaces for the resolution of sequences of nonlinear structural problems[END_REF][START_REF] André | A framework for deflated and augmented krylov subspace methods[END_REF][START_REF] Kirk | Krylov subspace recycling for sequences of shifted linear systems[END_REF][START_REF] Kapil | Recycling bicgstab with an application to parametric model order reduction[END_REF][START_REF] Carlberg | Krylov-subspace recycling via the pod-augmented conjugate-gradient method[END_REF][START_REF] Jolivet | Block iterative methods and recycling for improved scalability of linear solvers[END_REF].

Usually, deflation subspaces in recycling methods are approximations of eigenspaces. In [START_REF] Erhel | Restarted GMRES preconditioned by deflation[END_REF], the authors proposed an approximation of the eigenspace by using the Ritz vectors. In order to preserve the structure of the Arnoldi procedure, Morgan presented the deflation of harmonic Ritz vectors [START_REF] Morgan | GMRES with deflated restarting[END_REF]. A recent survey on recycling subspaces covers most developments in this area [START_REF] Kirk | A survey of subspace recycling iterative methods[END_REF]. Few papers discuss the impact of singular vectors on the convergence of Krylov methods [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF][START_REF] Meza | Deflated krylov subspace methods for nearly singular linear systems[END_REF][START_REF] Meza | A modification to the GMRES method for ill-conditioned linear systems[END_REF]. In [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF], Simoncini shows the impact of very small singular values on the convergence of several restarted variants of Krylov methods when the right-hand side has a large coefficient in the direction of the left singular vectors associated with the smallest singular values. Simoncini's work inspired us to consider deflation of approximate singular vectors rather than eigenvectors.

In the Hermetian case, Vuik et al. [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF] prove that for any deflation subspace the effective condition number of the deflated matrix is smaller than the effective condition number of the matrix without deflation. However, this is not true for the non-symmetric case. Moreover, we will see in the numerical experiments 8 how the deflation subspace, in several cases might deteriorate the convergence of the method. Thus, when considering a new matrix in the sequence of linear systems, it is important to check whether the deflation subspace is appropriate or not. For this aim, we introduce a new strategy to truncate the deflation subspace based on the approximation type of the deflation. This truncation strategy allows to select the useful part of the deflation subspace and removes the rest.

We introduce the deflation based on singular vectors. We derive the deflation formulation of singular vectors for a general linear systems. It turns out to be a special case of the general Inria deflation formula presented previously in [START_REF] André | A framework for deflated and augmented krylov subspace methods[END_REF].

Then we derive the approximation strategy for deflating the singular vectors. We compare the combination of the truncation of the deflation subspace and the deflation of singular vectors to the existing method GCRO-DR on a sequence of linear systems arising from reservoir simulation. Since we are interested in solving a sequence of linear systems, we follow the presentation of Parks et al. in [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF].

The paper is organized as follows. In the following section, we give a brief review about the Krylov subspace and the Arnoldi procedure, and then we present the deflation subspace strategy. Afterwards, we review the method GCRO-DR [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF] with an abstract deflation subspace related to the previous search subspace. To motivate the deflation of the singular vectors, an extension of a result obtained in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] on the influence of the singular values on the convergence of GMRES is extended to GCRO-DR and presented in Section 3 In Section 4, we introduce the deflation of singular vectors. We derive a formulation of deflation of exact singular vectors. Section 5 presents strategies of computing deflation vectors for GCRO-DR. At first, we recall the Ritz pairs and the harmonic Ritz pairs [START_REF] Paige | Approximate solutions and eigenvalue bounds from krylov subspaces[END_REF], and then we introduce a strategy of deflating approximate singular vectors. In Section 6, we present the truncation of the deflation subspace in recycling Krylov subspace methods. For each type of recycled vectors we propose a truncation strategy. Section 7 details the parallel implementation of the method and analyzes its computation and communication cost. In Section 8, we compare the three variants of deflation techniques combined with the truncation of the deflation subspace on a sequence of linear systems arising from reservoir simulation.

Background

In this section, we review the Krylov subspaces, the Arnoldi procedure, and the deflation of eigenvectors. We review the GCRO-DR algorithm to which we propose a modification in the following sections.

Notation

For an integer m > 0 we refer to the m × m identity matrix as I m . Let the integers m > 0, p > 0, and let h ij ∈ C, for i = 1, . . . , m, j = 1, . . . , p, we refer to the matrix whose (i, j) coefficient is h i,j for i = 1, . . . , m, j = 1, . . . , p as (h ij ) 1≤i≤m 1≤j≤p . If for some pair (i, j) such that 1 ≤ i ≤ m and 1 ≤ j ≤ p the element h i,j is not defined, it is set to zero. Let H = (h i,j ) 1≤i≤m 1≤j≤p ∈ C m×p . we refer to the matrix (h i,j ) 1≤i≤m-1 1≤j≤p ∈ C (m-1)×p as H. Unless otherwise stated, the matrix 1)×p is referred to as H. The m × p matrix with zeros elements is referred to as 0 m,p and the subscript is omitted when the size is clear from the context.

(h i,j ) 2≤i≤m 1≤j≤p ∈ C (m-

Krylov subspaces and Arnoldi procedure

The kth Krylov subspace associated to A and b is defined as for i = 1 : j do 4:

K k (A, b) = span b, Ab, . . . , A k-1 b . (2)
h ij = v H i w 5:
end for

6: w = w - j i=1 h ij v i 7:
h j+1,j = w 2

8:

if h j+1,j = 0 then 9:

Set m = j, v m+1 = 0

v j+1 = w
hj+1,j

12:

end if 13:

V j = [v 1 , • • • , v j ], V j+1 = [V j , v j+1 ], H j = (h ik ) 1≤i≤j+1 1≤k≤j
14: end for There are multiple orthogonalization procedures that can be used in the Arnoldi procedure. We refer the reader to [START_REF] Langou | Iterative methods for solving linear systems with multiple right-hand sides[END_REF][START_REF] Giraud | The loss of orthogonality in the gram-schmidt orthogonalization process[END_REF][START_REF] Giraud | Rounding error analysis of the classical gram-schmidt orthogonalization process[END_REF] for analysis on their behavior.

According to the notations in Algorithm 1 the following relation holds at each iteration k

AV k = V k+1 H k , (3)
where

V k = {v 1 , . . . , v k }, V k+1 = [V k , v k+1 ],
and

H k = (h i,j ) 1≤i≤k+1 1≤j≤k .

Deflation of eigenvectors

Let S be an A-invariant subspace of dimension k ≥ 0 related to the smallest eigenvalues in magnitude of the matrix A. Let Z ∈ C n×k be a matrix whose columns form an orthonormal basis of the subspace S, and let T be the projection of A on the subspace S, i.e., T = Z H AZ. We note that the following relation holds

AZ = ZT.
Let {λ 1 , . . . , λ n } be the eigenvalues of A ordered increasingly by the magnitude of their absolute value, i.e., |λ 1 | ≤ . . . ≤ |λ n |. Consider the preconditioned system

(I + Z(T -1 -I)Z H )Ax = (I + Z(T -1 -I)Z H )b (4)
The following Theorem 1 can be found in [START_REF] Erhel | Restarted GMRES preconditioned by deflation[END_REF]. It describes how to deflate the eigenvalues of the matrix A by adding a low-rank correction.

Theorem 1 The eigenvalues of the matrix A + Z(T -1 -I)Z H A are {1, . . . , 1, λ k+1 , . . . , λ n }.

Inria

The equation ( 4) can be rewritten as :

A -ZZ H A + ZT -1 Z H A x = b -ZZ H b + ZT -1 Z H b
In other words, to solve (1), we look for the solution in two subspaces, the first is S, the second is the orthogonal subspace to S. The first finds the part of the solution that is generated by the basis vectors V . This is done by solving the deflation subspace system with the matrix T . The second finds the part that is orthogonal to S. This part is left to the Krylov method. Since the eigenvectors of the matrix are not orthogonal in general, we need at each iteration of the Krylov method to correct the new computed basis vector. In the Hermitian case, this correction can be avoided by choosing a specific initial guess, see [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF].

Ritz pairs

In this paragraph we review the definition of Ritz pairs of a matrix related to a subspace.

Definition 1 Let B ∈ C n×n and let S be a subspace of dimension k < n. We say that the pair

(u, θ) ∈ C n × C is a Ritz pair of B related to the subspace S if u ∈ S, y H (Bu -θu) = 0, ∀y ∈ S. (5)
In the case where the subspace S = K k (A, b), the Ritz pairs of the matrix A correspond to

R = {(V k u, θ), H k u = θu}, see Algorithm 1. If we consider the subspace S = AK k (A, b), the Ritz pairs of the matrix A -1 correspond to R = {(V k u, θ), H H k u = θH H k H k u}.
Lemma 1 Let B ∈ C m×m and let S be a subspace of dimension k < m that contains an eigenpair (λ, u) of the matrix B. Let V ∈ C m×k be a matrix whose columns form an orthonormal basis of S. Then, there exists a vector w ∈ C k such that (λ, w) is an eigenpair of V H BV , the restriction of B to the subspace S. Furthermore, w = V H u.

Proof 1 The proof is simple and direct.

In Section 5, we will review in detail the derivation of these Ritz pairs. Furthermore, we derive the Ritz pairs relative to the singular vectors and values of A.

In the following, we review the GCRO-DR method and the construction of an abstract deflation subspace related to it.

GCRO-DR

We recall that GCRO-DR is an iterative method for solving a sequence of linear systems which takes advantage of the history of previous search subspaces by performing recycling. In [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF], the authors present the algorithm with a deflation subspace related to the harmonic Ritz vectors [START_REF] Paige | Approximate solutions and eigenvalue bounds from krylov subspaces[END_REF] as well as another deflation subspace that was previously proposed in [START_REF] Sturler | Truncation strategies for optimal krylov subspace methods[END_REF]. In the following, we explain briefly the method. We note that most of the following algebraic results can be found in [START_REF] De Sturler | Nested krylov methods based on gcr[END_REF]. Given the sequence of linear systems A i x = b i , for i = 1, . . . , p, where p > 1, the method GCRO-DR(m, k) works as the following. To solve the linear system A 1 x = b 1 , the first cycle, which is a sequence of iterations between two restarts, is performed as a usual GMRES cycle. The subspace relation that holds at the end of the cycle is the usual Arnoldi relation :

A 1 V m = V m+1 H m .
The solution x 1 and the residual r 1 are computed as usual. Based on a given strategy, a column matrix P k ∈ C m×k is computed (Section 5 details the computation of this matrix following different strategies.) This matrix is expanded by V m to form the basis of the deflation subspace Ỹk = V m P k . The image of the deflation subspace by A 1 is computed

A 1 V m P k = V m+1 H m P k , then orthonormalized, C k = V m+1 Q,
where Q stands for the unitary factor of the QR factorization of H m P k = QR. The relation that holds at this stage is :

A 1 Ỹk R -1 = C k .
We set U k = Ỹk R -1 . In order to avoid round-off errors, the vectors of U k are normalized. We define the diagonal matrix D k such that the following holds

Ũk = U k D -1 k , A 1 Ũk = C k D k .
At that moment, the method is ready to carry on the next cycle, m -k Arnoldi iterations are performed with (I -C k C H k )A 1 as an operator and the normalized residual as the starting vector. The following relation holds

A 1 [ Ũk , V m-k ] = [C k , V m-k+1 ]G m ,
where V m-k+1 is the resulting orthonormal basis constructed by the Arnoldi procedure,

G m = D k B m-k 0 H m-k , B m-k = C H k A 1 V m-k
and H m-k is the resulting Hessenberg matrix of the Arnoldi procedure. The solution x and the residual r are computed following the constraint r 2 2 = r 1 -A 1 x 2 has a minimal norm over the subspace spanned by the columns of [ Ũk , V m-k ]. The extraction of the deflation subspace is performed as before by replacing the matrix H m by G m . This cycle is repeated until the convergence is achieved. In order to solve the second linear system, all what is necessary is to update the image of the deflation subspace, i.e., C k , by the image of the deflation subspace by the new matrix A 2 . Algorithm 2 presents the GCRO-DR method with an abstract deflation subspace. We will discuss the deflation variants in Section 5.

Considering the notations in Algorithm 2, if the matrix Ỹk is not defined, then (3) holds. Otherwise, the following relation holds

A[ Ũk , V m-k ] = [C k , V m-k+1 ] D k B m-k 0 H m-k . (6)
The matrix [C k , V m-k+1 ] is unitary, but this is not true, in general, for the matrix [ Ũk , V m-k ].

In this case the subspace generated by the columns of Ũk , V m-k is not a Krylov subspace.

Inria

Algorithm 2 GCRO-DR Require: the maximum dimension of the search subspace m, the dimension of the deflation subspace k, the convergence tolerance ε, the initial guess x 0 . Ensure: approximate solution x of Ax = b.

1: initialization r 0 = b -Ax 0 , i = 1 2: if Ỹk is defined (from solving a previous linear system) then

3: let [Q, R] be the QR factorization of A Ỹk 4: C k = Q 5: U k = Ỹk R -1
6:

x 1 = x 0 + U k C H k r 0 7: r 1 = r 0 -C k C H k r 0 8: else 9: v 1 = r 0 / r 0 10: c = r 0 2 e 1 11:
perform m steps of GMRES, solving min c -H m y 2 for y and generating V m+1 and H m 12:

x 1 = x 0 + V m y 13: r 1 = V m+1 (c -H m y) 14:
compute k deflation vectors P k 15:

Ỹk = V m P k 16: let [Q, R] be the QR factorization of H m P k 17: C k = V m+1 Q 18: U k = Ỹk R -
Vm = [ Ũk , V m-k ], Ŵm+1 = [C k , V m-k+1 ] 25: G m = D k B m-k 0 H m-k 26:
solve min G m y -Ŵ H m+1 r i-1 2 for y 27:

x i = x i-1 + Vm y 28:

r i = r i-1 -Ŵm+1 G m y 29:
compute k deflation vectors P k 30:

Ỹk = Vm P k 31: let [Q, R] be the QR factorization of G m P k 32: C k = Ŵm+1 Q 33: U k = Ỹk R -1 34: end while 35: let Ỹk = U k 3 On the convergence of GCRO-DR
Here, we extend a theoretical result obtained in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] on the convergence of GMRES to GCRO-DR with an abstract recycling subspace. We suppose in the following that the recycling subspace is defined. Without loss of generality, we suppose that the residual is normalized. At iteration j of the GCRO-DR algorithm, the minimization problem in the search subspace is :

y j = arg min y∈C j+k e k+1 -G j+k y 2 .

The following proposition generalizes [36, Proposition 4.1]

Proposition 1 At iteration j of the GCRO-DR algorithm (the inner GMRES iteration j with the the matrix coefficient

(I -C k C H k )A)
, the following holds : -the solution of the minimization problem y j is :

y j = (1 + H -H j h j 2 2 ) -1 G -1 j+k G -H j+k hj . (7)
-the projection of r j on v 1 is

v H 1 r j = (1 + H -H j h j 2 2 ) -1 (8)
-the projection of

r j on [C k , v 2 , . . . , v j+1 ] is [C k , v 2 , . . . , v j+1 ] H r j = (1 + H -H j h j 2 2 ) -1 0 H -H j h j (9)
-cos 2 (θ(v 1 , r j )) = r j 2 2 = (1 + H -H j h j 2
2 ) -1 where v 1 , . . . , v j+1 are the basis vectors obtained by the Arnoldi procedure,

h j = [0, h j ] ∈ R k+j , G j+k = D k B j 0 H j
, and θ(v 1 , r j ) is the angle between the vectors v 1 and r j . 1 .

Proof 2

y j = (G H j+k G j+k ) -1 G H j+k e k+1 , = (G H j+k G j+k + hj hH j ) -1
hj ,

Using the Sherman-Morrison formula, we have

y j = (1 + G -H j+k hj 2 2 ) -1 G -1 j+k G -H j+k hj ,
Exploiting the structure of hj and G j+k we can write

(G -H j+k hj ) = [0, H -H j h j ]
. Then, we have the relation required,

y j = (1 + H -H j h j 2 2 ) -1 G -1 j+k G -H j+k
hj . Furthermore, we can write

y j = (1 + H -H j h j 2 2 ) -1 -D -1 k B j H -1 j H -H j h j H -1 j H -H j h j .
The residual in GCRO-DR can be expressed as :

r j = b -A[ Ũk , V j ]y j , = [ Ck , V j+1 ][e k+1 , G j+k ] 1 -y j , = [C k , v 1 , v 2 , . . . , v j+1 ]   0 D k B j 1 0 h H j 0 0 H j   1 -y j , = [v 1 , C k , v 2 , . . . , v j+1 ]   1 0 h H j 0 D k B j 0 0 H j   1 -y j , = [v 1 , C k , v 2 , . . . , v j+1 ] 1 hH j 0 G j+k 1 -y j , = v 1 (1 -hH j y j ) -[C k , v 2 , . . . , v j+1 ]G j+k y j .
Combining the last relation with (7), we have

v H 1 r j = 1 -(1 + H -H j h j 2 2 ) -1 G -H j+k hj 2 , = 1 -(1 + H -H j h j 2 2 ) -1 H -H j h j 2 , = (1 + H -H j h j 2 2 ) -1 .
and

[C k , v 2 , . . . , v j+1 ] H r j = G j+k y j , = (1 + H -H j h j 2 2 ) -1 G -H j+k hj , = (1 + H -H j h j 2 2 ) -1 0 H -H j h j .
The norm of the residual r j is given as :

r j 2 2 = (1 + H -H j h j 2 2 ) -2 + (1 + H -H j h j 2 2 ) -2 H -H j h j 2 2 = (1 + H -H j h j 2
2 ) -1 .

To prove the last relation, it is sufficient to note that :

v H 1 r j = cos(θ(v 1 , r j )) r j 2
Proposition 1 shows that the norm of the residual starts to decrease when the value

H -1 j h 2 2
is large. Moreover, if h j 2 is very small and H -1 j 2 is not sufficiently large the method will stagnate. As discussed in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF], the value of h j 2 may depend on different factors. A focus in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] is on the case when the right-hand side has large components corresponding to the left singular vectors associated with small singular values of the coefficient matrix. For example, let

(u k , σ k , w k ) be the kth singular triplet of A ; if v 1 = u n , then h H j = σ n w H n V j and h j 2 ≤ σ n . This can be easily generalized to the case v 1 = α n u n + . . . + α k u k for some k ∈ 1; n -1 , where n i=k α 2 k = 1
, and α j α n , for k ≤ j < n. Suppressing these components by deflating the right singular vectors associated with the smallest singular values avoids the previous situation.

In the following section we discuss the deflation of exact singular vectors.

Deflation based on singular vectors

In the previous section we extended a result obtained in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] on the convergence of GMRES to the GCRO-DR method. As previously mentioned, the work of Simoncini on the convergence of several restarted Krylov methods and its relation with the singular vectors associated to the smallest singular values [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] motivated us to investigate the deflation based on singular vectors approximation. This section introduces the deflation of singular vectors. Given a set of singular vectors (right or left), we show how to deflate these vectors during the Krylov method.

Theorem 2 Let x * be the exact solution of (1). Let A = U ΣV H be the singular value decomposition of A such that the singular values are ordered increasingly. Let k τ be the number of singular values smaller than a given threshold τ > 0. Consider the splitting of the SVD decomposition

A = U 1 Σ 1 V H 1 + U 2 Σ 2 V H 2 ,
where Σ 2 ∈ C kτ ×kτ is a diagonal matrix whose diagonal elements are the singular values smaller than τ . Consider x an approximate solution of the following linear system of equations

(I -U 2 U H 2 )Ax = (I -U 2 U H 2 )b, (10)
such that x -x 2 ≤ ε, where x is an exact solution of (10) and ε > 0. Then, the error norm in the solution of the original system, corresponding to the approximate solution of the associated projected system (10) and the correction term, is bounded by the error norm in the solution of the projected system (10), i.e., the following inequality holds

x * -(I -V 2 V H 2 )x -V 2 Σ -1 2 U H 2 b 2 ≤ ε.
Proof 3 First, we remark that x * is a solution of (10) and the set of solutions of (10) can be written as S = {x = x * + V 2 u, u ∈ C kτ }. Indeed, let x be a solution of [START_REF] Yogi | Deflation and balancing preconditioners for krylov subspace methods applied to nonsymmetric matrices[END_REF]. Then,

Ax = (I -U 2 U H 2 )b + U 2 U H 2 Ax, Ax = Ax * -U 2 U H 2 (b -Ax).
We multiply by A -1 both sides. We obtain

x = x * + V 2 Σ -1 2 U H 2 (Ax -b) = x * + V 2 Σ -1 2 (Σ 2 V H 2 x -U H 2 b) = x * + V 2 V H 2 x -V 2 Σ -1 2 U H 2 b.
On the other side, each element of S is a solution of (10).

We have

V 2 V H 2 x * = V 2 V H 2 A -1 b = V 2 Σ -1 2 U H 2 b.
Hence, we can write

x * = (I -V 2 V H 2 )x * + V 2 V H 2 x * , = (I -V 2 V H 2 )x + V 2 Σ -1 2 U H 2 b, Inria thus, x * -(I -V 2 V H 2 )x -V 2 Σ -1 2 U H 2 b = (I -V 2 V H 2 )(x -x).
Finally, we obtain

x * -(I -V 2 V H 2 )x -V 2 Σ -1 2 U H 2 b 2 ≤ (I -V 2 V H 2 )(x -x) 2 , ≤ x -x 2 , ≤ ε. Note that the matrix (I -U 2 U H
2 )A is singular. However, solving the linear system of equations ( 10) with a Krylov method is possible. The approximate solution of (1) can be computed as

x = (I -V 2 V H 2 )x + V 2 x 2 ,
where

x 2 = Σ -1 2 U H 2 b
and x is the approximate solution given by the Krylov method. In [START_REF] André | A framework for deflated and augmented krylov subspace methods[END_REF], the authors presented a general formulation of deflation by using an abstract deflation subspace. Theorem 2 is a special case for a subspace associated with singular vectors.

The following section discusses different approaches for computing deflation information from a Krylov subspace. Two approaches are based on approximating eigenvectors of the matrix, and these approaches can be found in literature [START_REF] Morgan | GMRES with deflated restarting[END_REF][START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF]. A third approach is introduced here and it is based on approximating the right singular vectors of the matrix.

Computing deflation vectors

In this section, we derive three different types of deflation subspaces. These subspaces are considered in the context of GCRO-DR [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF]. In Section 2, we reviewed the method GCRO-DR for an abstract deflation subspace issued from the last Krylov subspace. Even though any deflation vectors can be chosen for GCRO-DR, this choice has an impact on the convergence of the method. Several types of deflation vectors have been presented in previous works [START_REF] Morgan | GMRES with deflated restarting[END_REF][START_REF] Erhel | Restarted GMRES preconditioned by deflation[END_REF][START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF]. We focus on two types which are widely used and based on the smallest and largest Ritz values of the matrices A, A -1 , respectively. During the GCRO-DR method, one of the relations (3) and ( 6) holds. In the following, we present how to compute the deflation vectors mentioned above. These strategies were introduced in previous works [START_REF] Morgan | GMRES with deflated restarting[END_REF][START_REF] Erhel | Restarted GMRES preconditioned by deflation[END_REF][START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF]. We present them for the sake of completeness. Afterwards, we introduce the deflation vectors that are related to the right singular vectors.

We consider the following environment. Let K be a search subspace (Krylov subspace or else) of dimension m associated to the matrix A, let V ∈ C n×m be a matrix whose columns form a basis of K. Let W ∈ C n×(m+1) be a unitary matrix such that the following relation holds

AV = W G, (11)
where G ∈ C (m+1)×m is a Hessenberg matrix. We will derive the eigenvalue problem that is solved in order to compute the deflation vectors. For each approximation strategy, we formulate the computation following the notations in Algorithm 2. Thus, we will distinguish two cases in each strategy, depending on whether the deflation subspace already exists or not.

Approximation based on the smallest Ritz values of A

In this case, the Ritz pairs (u, θ) verify

u ∈ K, y H (Au -θu) = 0, ∀y ∈ K.
Since the columns of V span the subspace K, we have

V H (AV w -θV w) = 0,
where u = V w. This can be written as

V H W Gw = θV H V w.
By adapting the notations in Algorithm 2, the deflation vectors can be computed by solving the following eigenvalue problem

     H m w = θw, if Ỹk is not defined Ũ H k C k Ũ H k V m-k+1 0 (I m-k , 0) G m w = θ Ũ H k Ũk Ũ H k V m-k V H m-k Ũk I m-k w, if Ỹk is defined. ( 12
)

Approximation based on the largest Ritz values of A -1

To approximate the eigenvectors of A -1 we look for the approximate vectors in the subspace AK. We write

u ∈ AK, y ⊥ (A -1 u -θu), ∀y ∈ AK,
The columns of AV span the subspace AK. Hence, we have

V H A H (V w -θAV w) = 0,
this can be written as

G H W H V w = θG H Gw.
By using the notations of Algorithm 2, this requires solving the following generalized eigenvalue problem

           H H m w = θH H m H m w, if Ỹk is not defined G H m    C H k Ũk 0 V H m-k+1 Ũk I m-k 0    w = θG H m G m w, if Ỹk is defined. ( 13
)
We note that Morgan, in his paper [START_REF] Morgan | GMRES with deflated restarting[END_REF], refers to the values θ -1 as the harmonic Ritz values.

Inria

Approximation based on the smallest Ritz values of A H A

We follow the definition of the Ritz pairs given in Definition 1. In order to approximate the right singular vectors of the matrix A, we consider an approximation of the eigenvectors of the matrix A H A. Thus, we are looking for the pair (u, θ) such that

u ∈ K, y ⊥ (A H Au -θu), ∀y ∈ K
The columns of V span the subspace K, the previous relation can be written as

V H (A H AV w -θV w) = 0,
where u = V w. This can be written as

G H Gw = θV H V w.
The formulation of this eigenvalue problem in the context of Algorithm 2 is as follows

     H H m H m w = θw, if Ỹk is not defined G H m G m w = θ Ũ H k Ũk 0 0 I m-k w, if Ỹk is defined. ( 14
)
6 Deflation subspace truncation

In this section, we introduce a strategy to reduce the dimension of the deflation subspace. Consider a sequence of linear systems solved by GCRO-DR. This method builds or updates a deflation subspace during the solution of each linear system in the sequence. This deflation subspace is used to potentially enhance the convergence of the iterative method while solving the following system. However, part of the deflation subspace may not be useful for improving the convergence when the matrix changes. When the matrices in the sequence are Hermitian, the worst case that can happen is that the deflation subspace does not enhance the convergence ; see [START_REF] Tang | Comparison of two-level preconditioners derived from deflation, domain decomposition and multigrid methods[END_REF]. Nonetheless, this is not the case for non-Hermitian matrices, where an improper choice of the deflation subspace may even deteriorate the convergence. For this reason, one needs a strategy to extract the useful and relevant information from the deflation subspace and removes the improper part of it. Such a strategy should keep in the truncated deflation subspace a vector that is a good approximate of an eigenvector associated with a small eigenvalue in magnitude of A, A H A, or large eigenvalue in magnitude of A -1 depending on the deflation strategy.

One approach can be based on the accuracy of the approximated eigenvectors. Nevertheless, it was noticed in [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF][START_REF] Michael | The Iterative Solution of a Sequence of Linear Systems Arising from Nonlinear Finite Element Analysis[END_REF] that GCRO-DR does not require high accuracy of the approximated eigenvectors in order to improve the convergence. A practical selection criterion for an accuracybased approach needs the estimation of the residual. This may be a very strong criterion since it requires not only a good approximation of the eigenvectors but the eigenvalues as well. Hence, this can lead to remove the deflation subspace entirely even though part of it might be useful.

The strategy that we propose is based on Ritz values of the deflation strategy matrix (A, A H A, or A -1 ) with respect to the deflation subspace.

Let B be the matrix which defines the deflation strategy, B = A, A H A, or A -1 . Suppose that the deflation subspace contains an exact eigenvector u of B associated with the eigenvalue θ that is smaller in magnitude than τ, τ 2 , or larger in magnitude than 1 τ , respectively, where τ > 0 is a predefined value. Lemma 1 proves that the projection of the vector u to the deflation subspace is an eigenvector of the projected matrix to the deflation subspace. Therefore, on one hand, this strategy allows to keep in the truncated deflation subspace a vector that is a good approximate of an eigenvector associated with a small eigenvalue in magnitude of A, A H A, or large eigenvalue in magnitude of A -1 depending on the deflation strategy. On the other hand, this strategy relies only on the Ritz values.

The threshold τ can be chosen to be a relative value to either the magnitude of the largest approximated eigenvalue or the largest approximated singular value depending on the deflation strategy. The truncation strategy has to be associated with the recycling strategy, i.e., if the recycling strategy deflates the Ritz vectors associated with the smallest (resp. largest) Ritz values, the truncation strategy disregards the Ritz vectors associated with the largest (resp. smallest) Ritz values. The truncated subspace is then a subset of the original deflation subspace.

In the following, we derive the generalized eigenvalue problems that are necessary to perform the truncation of the deflation subspace. We suppose that the columns of the unitary matrix Z span the deflation subspace S. Table 1 presents the eigenvalue problems to be solved in order to truncate the deflation subspace corresponding to each strategy of deflation.

Reference matrix Eigenvalue problem Disregarded part

A Z H AZw = θw largest magnitude A -1 Z H A H Zw = θZ H A H AZw smallest magnitude A H A Z H A H AZw = θw largest magnitude
Table 1 -Generalized eigenvalue problems truncating the deflation subspace, Reference matrix stands for the matrix whose eigenvectors are approximated, Eigenvalue problem refers to the problem to solve in order to truncate the deflation subspace, Disregarded part refers to the part of the spectrum to which are associated the disregarded vectors.

We note that for the case of an approximation of singular vectors it is more robust to perform an SVD factorization of AZ rather than solving the eigenvalue problem Z H A H AZw = θw.

Algorithm 3 Deflation subspace truncation

Require: matrix A, deflation basis vectors Ỹ , threshold τ , def the matrix to which the deflation subspace is related Ensure: truncated deflation basis vectors Ỹ

1: compute the QR factorization of Ỹ = ZR 2: if def = A then 3: solve Z H AZw = θw 4: else if def = A -1 then 5: solve Z H A H AZw = θZ H A H Zw 6: else if def = A H A then 7:
solve Z H A H AZw = θ 2 w 8: end if 9: form Ỹ whose columns are {Zw | w is an eigenvector associated with θ < τ } Here, we propose the method Generalized Minimal Residual with Modified Deflated Restarting (GMRES-MDR), Algorithm 4. This method is based on the GCRO-DR method. The main difference is that an estimation of the previous deflation subspace (when the matrix in the sequence changes) is performed in order to keep only relevant vectors.

Inria

Algorithm 4 GMRES-MDR Require: the maximum dimension of the search subspace m, the maximum dimension of the deflation subspace k max , the convergence tolerance ε, the initial guess x 0 , the matrix to which the deflation subspace is related def , τ threshold of truncation of the deflation subspace. Ensure: approximate solution x of Ax = b.

1: initialization r 0 = b -Ax 0 , i = 1, k = k max 2:
if Ỹk is defined (from solving a previous linear system) then

C k = Q 6: U k = Ỹk R -1
7:

x 1 = x 0 + U k C H k r 0 8: r 1 = r 0 -C k C H k r 0 9: else 10: v 1 = r 0 / r 0 11: c = r 0 2 e 1 12:
perform m steps of GMRES, solving min c -H m y 2 for y and generating V m+1 and H m 13:

x 1 = x 0 + V m y 14: r 1 = V m+1 (c -H m y) 15:
compute k deflation vectors P k following def and Section 5 16:

Ỹk = V m P k 17: let [Q, R] be the QR factorization of H m P k 18: C k = V m+1 Q 19: U k = Ỹk R -
Vm = [ Ũk , V m-k ], Ŵm+1 = [C k , V m-k+1 ] 26: G m = D k B m-k 0 H m-k 27:
solve min G m y -Ŵ H m+1 r i-1 2 for y 28:

x i = x i-1 + Vm y 29:

r i = r i-1 -Ŵm+1 G m y 30:
set k = k max compute k deflation vectors P k following def and Section 5

31: In this section, we describe the parallel design of the method GMRES-MDR. At first, we describe the distribution of data across processors. Then, we analyze the cost of the method and compare the three deflation variants which we use in our implementation.

Ỹk = Vm P k 32: let [Q, R] be the QR factorization of G m P k 33: C k = Ŵm+1 Q 34: U k = Ỹk R -1 35: end while 36: let Ỹk = U k V m-k C k Ũk

Data distribution

Matrices and vectors in our implementation can be separated into two types, distributed and redundant. The matrix A and the search subspace (Krylov and deflation) vectors are distributed in block row wise. All matrices and vectors whose both dimensions are less than or equal to the dimension of the search subspace are stored redundantly on each processor. Using this distribution, the application of the sparse operator requires communication between neighbour processors. Only inner product operations require global communication.

In the following, we describe how data is stored in memory.

Memory management

Dense matrices in iterative solvers are usually stored in one of two ways, either in column major or in row major. We use column major format due to constraints in the reservoir simulator. The vectors of the search subspace are stored in full storage scheme since they do not have a special structure. We refer the reader to [START_REF]Math kernel library[END_REF] for more details on storage scheme structures. In our implementation, all dense matrices are stored in full storage scheme. Figure 1 shows how the vectors of the basis are stored in memory. This order Ũk , C k , V m-k , in column ordering, has the following advantage. The orthogonalization steps in the Arnoldi procedure and the deflation correction can be done simultaneously. Sparse matrices are stored in the Compressed Sparse Row (CSR) storage scheme.

Parallel interaction environment and implementation language

We use message passing interface MPI [START_REF] Snir | MPI-The Complete Reference[END_REF] for communication between processors and we use OpenMP for multithreading. These two libraries are compatible with the language C that we chose to use.

Since GMRES-MDR, does not need an explicit coefficient matrix, we use a reverse communication interface (RCI) [START_REF] Dongarra | Reverse communication interface for linear algebra templates for iterative methods[END_REF]. Algorithms 5, 6 present the schemes of the RCI loop and the RCI routine, respectively. if ido == 1 then 5:

apply the preconditioned operator

+ k + 3) + m 2 k + m( 5 2 k 2 + k) -2 3 k 3

Cost analysis

Here, we compare the cost of one cycle of presented variations of GMRES-MDR to the cost of the classical GMRES method. Since the cost of the sparse matrix vector multiplication depends on the sparsity of the matrix A, we do not take it into account in our analysis. In the analysis of the cost of GMRES-MDR, we suppose that the deflation vectors exist i.e., starting from the second cycle. The computation costs consider the real-valued matrices and vectors. The computation costs of the generalized eigenvalue problems can be found in [START_REF] Golub | Matrix Computations[END_REF]. In Table 2 (resp. 4), we detail the computation (resp. communication) cost of one cycle of GMRES-MDR(m, k). Table 3 (resp. 5) shows a comparison between the computation (resp. communication) costs of the variants of GMRES-MDR and the classical GMRES. It is clear that the computation cost of the SVD variation is the smallest compared to the two other deflation strategies. This cost can be written as :

C(n, m, k) = n(m 2 +(2k +5)m+3k +3)+15m 3 +m 2 (k + 7 2 )+m( 5k 2 + 11 2 -3k)- 2 3 k 3 +2k 2 -5k.
To give an illustration of the difference between the computation cost of the SVD variant of GMRES-MDR and the classical GMRES, we can set m = 2k. Then the difference is

C(n, 2k, k) -C(n, 2k, 0) = n(5k 2 + 3k) + (9 - 2 3 )k 3 -4k 2 -5k.
However, it is sufficient that k > 1 to make the communication cost of GMRES-MDR smaller than the communication cost of the classical GMRES.

Algorithm 6 RCI routine

Require: ido the interface indicator, parameters of the solver 1: static variable status = 0 2: if ido == 0 then 3:

set up the solver using the parameters 4:

if deflation vectors exist then 5:

set ido = 1, status = -1 6:
return to apply the preconditioned operator on the deflation vectors prepare for the Arnoldi procedure return to apply the preconditioned operator on the last basis vector return to apply the preconditioned operator on the solution set ido = 99, status = 0 and return 43: end if

Inria Method cost(flops) GMRES(m) C(n, m, 0) -2n -15m 3 -m 2 GMRES-MDR(m, k) (Ritz A) C(n, m, k) + nk(m + 1) + 53m 3 + m 2 GMRES-MDR(m, k) (Ritz A -1 ) C(n, m, k) + nk(m -k) + 53m 3 + m 2 GMRES-MDR(m, k) (Ritz A H A) C(n, m, k)
Table 3 -Comparison of computation cost of one cycle of GMRES-MDR (supposing that k deflation vector exist) and classical GMRES. The maximum dimension of the search subspace is m. GEVP stands for generalized eigenvalue problem. Ritz A, Ritz A -1 , and Ritz A H A refer to the deflation strategy.

C(n, m, k) = n(m 2 + (2k + 5)m + 3k + 3) + 15m 3 + m 2 (k + 7 2 ) + m( 5k 2 +11 2 - 3k) -2 3 k 3 + 2k 2 -5k.
Operations messages words Preparation of the residual, the solution, and the starting vector 2 (k, 1) Arnoldi orthonormalization at iteration j 2 (j + k, 1) Prepare matrices for GEVP (Ritz A)

1 k(m + k + 1) Prepare matrices for GEVP (Ritz A -1 ) 1 k(m + 1)

Prepare matrices for GEVP (Ritz

A H A) 1 k 2 Prepare Ũk 1 k
Table 4 -Detailed communication costs of one cycle of GMRES-MDR supposing that k deflation vectors exist. The maximum dimension of the search subspace is m. GEVP stands for generalized eigenvalue problem. Ritz A, Ritz A -1 , and Ritz A H A refer to the deflation strategy. When the operation counts 2 messages the number of words is separated into two quantities corresponding to each message. The Arnoldi procedure uses the classical Gram-Schmidt variation.

Method messages GMRES(m)

2m

+ 1 GMRES-MDR(m, k) 2m -2k + 4
Table 5 -Comparison of communication cost of one cycle of GMRES-MDR (supposing that k deflation vector exist) and classical GMRES. The maximum dimension of the search subspace is m.

Numerical experiments

In this section, we illustrate the impact of the proposed method on the convergence of sequences of linear systems arising from reservoir simulation. We present sequential experiments performed in MATLAB as well as parallel experiments. Our set of test cases are obtained from the in-house prototype code at Total, which simulates a complex enhanced oil recovery (EOR) mechanism. This simulator relies on a finite volume discretization and a two-point flux approximation. We have two sequences of linear systems denoted with initials BIGCO24 and BIGP1. BIGCO24 systems arise from the simulation of water and gas injection using a compositional model (8 hydrocarbon components). The permeability field is heterogeneous. The grid has 83587 active cells. BIGP1 systems arise from the simulation of water injection using a black-oil model. The permeability field is heterogeneous (sector model from a real field case). The grid has 42332 active cells.

We are interested in applying recycling techniques in the context of the constrained pressure residual (CPR) preconditioner. At each Newton iteration of the reservoir simulation, a system of linear equations has to be solved. Typical preconditioner is the CPR preconditioner [START_REF] Cao | Parallel scalable unstructured CPR-Type Linear Solver for Resevoir Simulation[END_REF][START_REF] Wallis | Constrained residual acceleration of conjugate residual methods[END_REF]. This preconditioner is based on a two-level approach. For each application of the preconditioner on a vector, the restricted system associated with the pressure variable needs to be solved. Therefore, at each Newton iteration, several systems of linear equations associated with the pressure variable requires the solution. The right-hand side for each of them corresponds to the the restriction of a vector on the global level to the pressure level. This vector can be a Krylov basis vector, a solution, or any vector on which the user wants to apply the CPR preconditioner. On the other hand, the matrix is the same during the whole Newton iteration. It changes, however, from one Newton iteration to the following one. Typically, a low accuracy is sufficient for the restricted pressure system.

In our experiments we use the double classical Gram-Schmidt for the orthogonalization process in the Arnoldi procedure. Figure 2 shows the impact of deflating the singular vectors associated with the smallest singular values. In both variants, 20 vectors are computed approximately to a tolerance of 10 -6 by using the routines eigs and svds. The deflation of the eigenvectors corresponds to (4) which we rewrite here,

(I + Z(T -1 -I)Z H )Ax = (I + Z(T -1 -I)Z H )b,
where Z is a unitary matrix whose columns span the approximated eigenspace. The deflation of the singular vectors follows the Theorem 2. We write the preconditioned linear system

(I -U 2 U H 2 )Ax = (I -U 2 U H 2 )b,
where U 2 = AV 2 Σ -1 2 stands for the approximated left singular vectors. The solution is computed as

x = (I -V 2 V H 2 )x + V 2 x 2
, where x is the solution obtained by GMRES and x 2 = Σ -1 2 U H 2 b. We remark that the speed of convergence for both deflated variants are approximately the same with a small advantage for the SVD deflation. Furthermore, in the first iterations the SVD deflated GMRES has a sharp convergence curve while the curve corresponding to the eigenvectors deflated GMRES stagnates for a number of iterations. This phenomenon is important for our application. In each Newton iteration during the simulation of a reservoir, the linearized system is solved by using the CPR solver. An iteration of the latter solves the linear system associated with the pressure variable. This matrix does not change during the CPR solve, only the right-hand side changes. In practice, a large threshold of convergence for the pressure level is sufficient. Thus, having stagnation at the beginning of the method, even for a small number of iterations, can impact the global number of iterations over the entire simulation. The relative tolerance of convergence in the experiments related to the sequences of linear systems is fixed to 10 -1 .

As it was noted in the thesis of Parks [START_REF] Michael | The Iterative Solution of a Sequence of Linear Systems Arising from Nonlinear Finite Element Analysis[END_REF], the choice of the deflation subspace has an impact on the convergence of sequences of linear systems. Figure 3 does not only confirm this point but , respectively. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Fifteen linear systems in the sequence BIGP1 are solved, with three matrices, each solved for five different right-hand sides (for index i the three methods share the same right-hand side). In the figure on the left, the iteration counts are obtained by using a solver that does not truncate the deflation subspace. In the figure on the right, the iteration counts are obtained by using a solver that truncates the deflation subspace. The truncation threshold is set to 10 -3 . also shows how it might be possible that in some configurations the deflation subspace becomes an issue and leads to the stagnation of the residual norm. For each deflation variant, we solve a sequence of linear systems related to the pressure variable. We have three matrices arising from three contiguous Newton steps. For each matrix, we solve five (resp. ten) linear systems with initial BIGP1 (resp. BIGCO24) with different right-hand sides each given in a time. This mimics the CPR solver. In order to make a fair comparison and avoid eventual large components of the right-hand sides on special eigenvectors or singular vectors, we generate random right-hand sides with seeds corresponding to the index of the linear system. In the same figure we show the impact of the truncation of the deflation subspace. Figure 4 shows how the truncation strategy might lead to lose information that can enhance the convergence. We realize an increase in the number of iterations due to truncating the deflation subspace. One way to avoid such a case is to perform the truncation strategy after the first restart cycle. An activation flag to initiate the truncation strategy can be set if the residual norm reduction during the first cycle is larger than a predefined threshold that the user has to precise. Our numerical experiments lead to two remarks. The first is that deflating the Ritz vectors related to A is efficient when no deflation subspace exists. The second is that deflating the Ritz vectors related to A H A needs more initial iterations than the one related to A in order to be efficient. To benefit from both subspaces, we propose a simple criterion based on the reduction of the relative residual norm during the restart cycle to choose the deflation subspace adaptively. In terms of the Algorithm 4, if ri 2 ri-1 2 ≤ ε 0 , where ε 0 is a predefined threshold, then the SVD-based deflation strategy is chosen. Otherwise, the Ritz-based deflation strategy is chosen. In our experiments we set the value of ε 0 to 10 -1 . Figures 5 and6 present numerical experiments on the same two previous sequences. Figures 7 and8 present the history of the residual norm corresponding to each strategy of deflation. As we noted previously, we notice that the decrease of the residual norm starting from the first iteration makes the SVD variant have more advantage against other strategies. The total iteration numbers in the BIGCO24 sequence corresponding to the deflation related to A, A -1 , and the alternating (A, A H A) is 418, 374, and 163, respectively. We also notice that using the adaptive alternating deflation between the Ritz pairs of A and A H A has the benefit of both strategies, quick useful deflation subspace and quick reduction of residual norm. Moreover, in this example we see how , respectively. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Thirty linear systems in the sequence BIGCO24 are solved, with three matrices, each solved for ten different right-hand sides (for index i the three methods share the same righthand side). In the figure on the left, the iteration counts are obtained by using a solver that does not truncate the deflation subspace. In the figure on the right, the iteration counts are obtained by using a solver that truncates the deflation subspace. The truncation threshold is set to 10 -3 .

V svds (:, i) -QQ H V svds (:, i) 2 V eigs (:, i) -Q HR Q H HR V eigs (:, i) 2 V eigs (:, i) -Q R Q H R V eigs (:, i) 2 1 9.8 × 10 -2 7.6 × 10 -3 2.6 × 10 -2 2 2.0 × 10 -1 1.2 × 10 -2 2.3 × 10 -2 3 1.2 × 10 -1
3.8 × 10 -2 5.6 × 10 -2 4 1.5 × 10 -1 9.9 × 10 -1 9.9 × 10 -1 5 2.6 × 10 -1

1.0 × 10 -1 1.2 × 10 -1 the residual norm stagnates starting from the third matrix systems with the harmonic Ritzbased deflation subspace. The total iteration numbers in the BIGP1 sequence corresponding to the deflation strategies related to A, and the alternating (A, A H A) is 1062 and 686, respectively.

To measure the effectiveness of approximating the singular vectors and eigenvectors, we perform the following experiment :

-Compute the singular vectors (V svds ) of the operator that are associated with the smallest five singular values by using the routine svds. -Solve a linear system by using the SVD recycling strategy with GMRES-MDR [START_REF] Nicolaides | Deflation of conjugate gradients with applications to boundary value problems[END_REF][START_REF] De Sturler | Nested krylov methods based on gcr[END_REF] with convergence threshold 10 -6 . -Compute QR decomposition of the recycling subspace basis Ũk = QR,.

-Measure the error norm v -QQ H v 2 , where v is a column vector of V svds . We perform the same experiment with the two other variants and measure the error with respect to the eigenvector (V eigs ) associated with the smallest five eigenvalues in magnitude of the operator computed by using the routine eigs. Table 6 shows the results of this experiment. We observe that low approximation accuracy for each variant is sufficient. SVD, REV, and HEV stand for the deflation method related to the matrices A H A, A, and A -1 , respectively. The SVD variant uses an adaptive strategy that alternates the deflation between SVD and REV based on the relative residual norm during the cycle. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Fifteen linear systems in the sequence BIGP1 are solved, with three matrices, each solved for five different right-hand sides (for index i the three methods share the same right-hand side).

Inria i V svds (:, i) -QQ H V svds (:, i) 2 V eigs (:, i) -Q HR Q H HR V eigs (:, i) 2 V eigs (:, i) -Q R Q H R V eigs (:, i) 2 1 1.0 × 10 -6 1.2 × 10 -2 2.5 × 10 -2 2 2.6 × 10 -1 3.9 × 10 -3 7.1 × 10 -3 3 8.4 × 10 -2 7.9 × 10 -2 8.7 × 10 -2 4 2.0 × 10 -1 1.0 1.0 5 2.2 × 10 -1 3.3 × 10 -2 4.5 × 10 -2
Table 7 -Approximation measurement of recycling vectors to eigenspaces associated with the operator associated with the recycling strategy. Experiment is run on system with initial BIGCO24. The right-hand side is chosen to be the left singular vector associated with the smallest singular value. The number of iterations required to reach convergence is 185, 191, and 180 by using the SVD, harmonic Ritz, and Ritz recycling strategy, respectively.

In Table 7, we perform the same previous experiment where we choose the right-hand side to be the left singular vector associated with the smallest singular value. We note that the SVDbased recycling subspace extracted a good approximation of the right singular vector associated with the smallest singular value.

Parallel numerical experiments

In this section we present the numerical experiments of our implementation of the method GMRES-MDR. All experiments are performed on the supercomputer PANGEA at TOTAL. We had access to a maximum of 1024 nodes. Each node has 16 cores. Our set of test matrices consists of the matrices SPE10 and Grid2D-2000x2000. We solve two sequences, each of nine linear systems in which the matrix is fixed. The right-hand sides are different for each linear system in the sequence. All variants have the same right-hand side for the i th linear system, i = 1, . . . , 9. The matrix is SPE10 ∈ R n×n where n is the number of unknowns n = 1094421. The number of the non-zero elements nnz = 7478141. SPE10 corresponds to a simulation of a black oil with a heterogeneous permeability field [START_REF] Christie | Tenth spe comparative solution project : A comparison of upscaling techniques[END_REF]. The grid has 1094421 active cells. The matrix Grid2D-2000x2000 ∈ R n×n , where n is the number of unknowns n = 4 × 10 6 . The number of non-zero elements nnz = 2 × 10 7 . This matrix corresponds to a discretization of a pressure system in which the permeability field follows a log-normal distribution with standard deviation σ = 3 and a mean µ = log(500). The discretized domain is the unit square with 2000 × 2000 grid points. In all experiments, a linear system is solved to warm up the machine before solving each sequence. Figure 9 shows the scalability of the three variants discussed in Section 5. The sequence of linear systems is related to the test matrix Grid2D-2000x2000. The number of nodes (processors) varies from 16 up to 1024 with 2 cores per processor. A block Jacobi preconditioner with LU factorization in each block is used. The number of blocks is equal to the number of processors. In Listing 1 and Listing 2, we present two post-analysis of GMRES-MDR with SVD deflation variant that is associated with Figure 9. Figure 10 shows the scalability of the three variants discussed in Section 5. The sequence of linear systems is related to the test matrix SPE10. The number of nodes (processors) varies from 16 up to 256, 1 core per MPI processor. A block Jacobi preconditioner with ILU (0) factorization in each block is used. The number of blocks is equal to the number of processors. We notice that the three variants scale well with a smaller runtime for the SVD variant. We note that the runtime for the three variants of deflation does not scale on 512 nodes in the test case SPE10. See Figure 11 and the discussion below for the explanation. Figure 11 shows the percentage of the average time (over all processors) spent in the operations during the solution of a representative linear system of the sequence by using the SVD deflation strategy. This runtime corresponds to a solution of the system on 256 nodes. We remark that the time spent in the computation of the generalized eigenvalue problem becomes much more important than all other operations. This explains why the scaling of the runtime stopped for this sequence after 256 nodes.

Conclusion

In this work, we reviewed deflation strategies in Krylov iterative methods. Usual choice of deflation subspaces is related to the approximation of eigenvectors. Based on the study in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] we introduced the deflation of singular vectors. Since we are interested in solving a sequence of linear systems, we reviewed the method GCRO-DR [START_REF] Michael | Recycling krylov subspaces for sequences of linear systems[END_REF] with an abstract deflation subspace. We briefly presented two deflation subspaces which can be found in the literature. We extended part of the theory presented in [START_REF] Simoncini | On the convergence of restarted krylov subspace methods[END_REF] on the relation between the residual and the first vector in the Krylov basis sequence to GCRO-DR. Afterwards, we proposed a deflation subspace related to an approximation of singular vectors. To validate the proposed strategy, a comparison between the three variants was made on different sequences arising from reservoir simulation. This validation was performed on a sequential code as well as on a parallel code. The numerical experiments demonstrated the impact of choosing the SVD variant. The gain in the global iteration count is up to 40%. Concerning the parallel implementation, the three variants showed a similar behavior on a short sequence of linear systems. More experiments on larger test cases and longer sequences will be considered in future work.

Inria

Algorithm 1

 1 Let k A,b = min P ∈P P (0)=1 {deg(P ) | P (A)b = 0}, where P represents the set of polynomials with complex coefficients. The vectors b, Ab, . . . , A k-1 b for k ≤ k A,b form a basis of the Krylov subspace. Due to round-off errors, an orthogonalization procedure is necessary to construct an orthogonal basis of the Krylov subspace. To this aim, the Arnoldi procedure is widely used, see Algorithm 1. Arnoldi (A, v 1 , m) Require: normal vector v 1 ∈ C n , matrix A ∈ C n×n , number of iterations m Ensure: orthonormal basis vectors V m , Hessenberg matrix H m ∈ C (m+1)×m 1: for j = 1 : m do 2:w = Av j 3:

1 19: end if 20 :i = i + 1 22 :

 12022 while r i 2 > ε do 21: perform m -k Arnoldi steps with the linear operator (I -C k C H k )A, letting v 1 = r i-1 / r i-1 2 and generating V m-k+1 , H m-k , and B m-k = C H k AV m-k 23: let D k be a diagonal scaling matrix such that Ũk = U k D k , where the columns of Ũk have a unit norm 24:

3 :

 3 call Algorithms 3 and get Ỹk , where k is the number of columns in Ỹk 4: let [Q, R] be the QR factorization of A Ỹk 5:

1 20: end if 21 :i = i + 1 23 :

 12123 while r i 2 > ε do 22: perform m -k Arnoldi steps with the linear operator (I -C k C H k )A, letting v 1 = r i-1 / r i-1 2 and generating V m-k+1 , H m-k , and B m-k = C H k AV m-k 24: let D k be a diagonal scaling matrix such that Ũk = U k D k , where the columns of Ũk have a unit norm 25:

Figure 1 -

 1 Figure 1 -Deflation and search subspace basis vectors storage for GMRES-MDR method

Inria Algorithm 5 3 :

 53 RCI loop 1: ido = 0 2: while ido ! = 99 do call the reverse communication interface routine 6 and get the value of ido 4:

18 :

 18 set ido = 1, status = 2 and return 19: end if 20: if status = 2 then 21: orthogonalization of new basis vector and Hessenberg operations and stopping criterion 22: if convergence and max interior iteration are not achieved then 23: ido = 1, status = 2 24:

  if convergence or the max outer iteration is achieved then 29: ido = 2, status = 99 30: return to compute the unpreconditioned solution 31: else if max interior iteration is achieved then 32: ido = 1, status = 3 33:

Inria

Figure 2 -

 2 Figure 2 -Convergence history, singular vectors deflated GMRES against eigenvectors deflated GMRES. The convergence history of non-deflated GMRES is also plotted, results for the matrix BIGP1 in MATLAB 2017. The 20 computed approximated vectors (singular vectors and eigenvectors) are approximated up to a tolerance of 10 -6 .

Figure 3 -

 3 Figure3-Comparison between three variants of deflation subspaces. SVD, REV, and HEV stand for the deflation methods related to the matrices A H A, A and A -1 , respectively. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Fifteen linear systems in the sequence BIGP1 are solved, with three matrices, each solved for five different right-hand sides (for index i the three methods share the same right-hand side). In the figure on the left, the iteration counts are obtained by using a solver that does not truncate the deflation subspace. In the figure on the right, the iteration counts are obtained by using a solver that truncates the deflation subspace. The truncation threshold is set to 10 -3 .

Inria

Figure 4 -

 4 Figure4-Comparison between three variants of deflation subspaces. SVD, REV, and HEV stand for the deflation methods related to the matrices A H A, A and A -1 , respectively. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Thirty linear systems in the sequence BIGCO24 are solved, with three matrices, each solved for ten different right-hand sides (for index i the three methods share the same righthand side). In the figure on the left, the iteration counts are obtained by using a solver that does not truncate the deflation subspace. In the figure on the right, the iteration counts are obtained by using a solver that truncates the deflation subspace. The truncation threshold is set to 10 -3 .

Figure 5 -

 5 Figure5-Comparison between four variants of deflation subspaces. SVD, REV, and HEV stand for the deflation method related to the matrices A H A, A, and A -1 , respectively. SVDadap stands for an adaptive strategy that alternates the deflation between SVD and REV based on the relative residual norm during the cycle. Fifteen linear systems in the sequence BIGP1 are solved, with three matrices, each solved for five different right-hand sides (for index i the three methods share the same right-hand side).

Figure 6 -

 6 Figure6-Comparison between four variants of deflation subspaces. SVD, REV, and HEV stand for the deflation method related to the matrices A H A, A, and A -1 , respectively. SVDadap stands for an adaptive strategy that alternates the deflation between SVD and REV based on the relative residual norm during the cycle. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Thirty linear systems in the sequence BIGCO24 are solved, with three matrices, each solved for ten different right-hand sides (for index i the three methods share the same right-hand side).

  i

Figure 7 - 2 SVDFigure 8 -

 728 Figure7-History of residual norm in a sequence of linear systems (a peak corresponds to the beginning of a new linear system). Comparison between three variants of deflation subspaces. SVD, REV, and HEV stand for the deflation method related to the matrices A H A, A, and A -1 , respectively. The SVD variant uses an adaptive strategy that alternates the deflation between SVD and REV deflation methods based on the relative residual norm during the cycle. The maximum dimension of the search subspace is 30. The maximum dimension of the deflation subspace is 5. Thirty linear systems in the sequence BIGCO24 are solved, with three matrices, each solved for ten different right-hand sides (for index i the three methods share the same right-hand side).

Figure 11 -

 11 Figure 11 -Average time (over processors) of different operations in GMRES-MDR as a percentage of the total runtime. The communication and the computation part of the deflation and the orthogonalization are separated. This runtime corresponds to the solution of a representative linear system during the sequence. Number of nodes = 256. Deflation strategy is the SVD.

  Prepare matrices and solve GEVP (Ritz A) n(mk + k 2 + k) + 68m 3 + 2m 2 Prepare matrices and solve GEVP (Ritz A -1 ) nmk + 68m 3 + 2m 2 Prepare matrices and solve GEVP (Ritz A H A) nk 2 + 15m 3 + m 2 Prepare C k and Ũk nk(2m

	6:	end if
	7:	if ido == 2 then
	8:	apply the preconditioner
	9:	end if
	10: end while
	Operations	cost (flops)
	Preparation of the residual, the solution, and the starting vector n(4k + 3)
	Arnoldi orthonormalization	n(m 2 + 4(m -k) -k 2 )
	Hessenberg and residual updates	2m 2 -4mk + 2k 2 + 5(m -k)
	Solve the least squares problem and expand the solution	nm + 1 2 (m 2 + m)

Table 2 -

 2

Detailed computation cost of one cycle of GMRES-MDR supposing that k deflation vector exist. The maximum dimension of the search subspace is m. GEVP stands for generalized eigenvalue problem. Ritz A, Ritz A -1 , and Ritz A H A refer to the deflation strategy.

Table 6 -

 6 Approximation measurement of recycling vectors to eigenspaces associated with the operator associated with the recycling strategy. Experiment is run on system with initial BIGCO24. The number of iterations required to reach convergence is 204, 204, and 197 by using the SVD, harmonic Ritz, and Ritz recycling strategy, respectively.

RR n°9206

Note that G j+k is not the matrix obtained from G j+k by suppressing the first line as for H j . Inria

rue Simone Iff -CS 42112

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

 ------------------------------------------------------------------------Post a n a l y s i s : ----------------------------------------------------------------------- ---------------------------------------------------------------------- ---------------------------------------------------------------------- ------------------------------------------------------------------------Post a n a l y s i s : ----------------------------------------------------------------------- ---------------------------------------------------------------------- ---------------------------------------------------------------------


Pierre Lewden
email: pierre.lewden@ims-bordeaux.fr

Adrien F Vincent
email: adrien.vincent@ims-bordeaux.fr

Charly Meyer
email: charly.meyer@ims-bordeaux.fr

Jean Tomas
email: jean.tomas@ims-bordeaux.fr

Shidoush Siami
email: ssiami@univ-pau.fr

Sylvain Saïghi
email: sylvain.saighi@ims-bordeaux.fr

Hardware Spiking Neural Networks: Slow Tasks Resilient Learning with Longer Term-Memory Bits

Keywords: neuromorphic systems, spiking neural networks, memristors, spike timing-dependent plasticity

When aiming at efficient and low-power processing of event-based data, hardware implementations of spiking neural networks that co-integrate analog silicon neurons with memristive synaptic crossbar arrays are a promising framework. Fully analog systems however commonly make it difficult to learn patterns with real-world timescales, which are typically beyond the millisecond, due to intrinsic constraints of the underlying technologies. In this work, we propose to alleviate this issue by supplementing each presynaptic unit with a single memory bit, which allows to implement a hardware-friendly Spike Timing-Dependent Plasticity. By simulation means on the N-MNIST dataset, we illustrate the potential of this concept and show its robustness to postsynaptic neuron variability. We also discuss how to circumvent challenges raised by initial weight distribution. These results could facilitate the emergence of embedded smart systems directly fed by event-based sensors.

I. INTRODUCTION

Brain-inspired neuromorphic systems have attracted significant attention in the recent years [START_REF] Indiveri | Memory and information processing in neuromorphic systems[END_REF], [START_REF] Thakur | Large-scale neuromorphic spiking array processors: A quest to mimic the brain[END_REF], offering better performance on cognitive memory-intensive tasks than conventional architectures limited by the von Neumann bottleneck [START_REF] Editorial | Big data needs a hardware revolution[END_REF]. In particular, leveraging Kirchhoff's laws and Ohm's law, neuromorphic architectures that rely on memristor-based synaptic arrays are candidates of choice to build highly integrated and low-power hardware implementations of Spiking Neural Networks (SNNs) with on-line learning capabilities. Among the possible learning rules, Spike Timing-Dependent Plasticity (STDP) naturally fits such hardware event-based systems and allows unsupervised learning, a major challenge for the processing of increasingly large volumes of data [START_REF] Lecun | Deep learning[END_REF].

Implementing STDP in low-power analog or mixed-signal hardware however raises several challenges. First, depending on its complexity, the chosen STDP algorithm may not be hardware-friendly, inducing a consequential circuit or power overhead. Besides, the timescales that reasonably fit into integrated analog hardware are typically orders of magnitude shorter than the time constants of many real-time tasks.

Among the promising research works on hardware spiking neural networks that investigate fully CMOS solutions, some Financial support from the French Agence Nationale de la Recherche (ANR) through MIRA project is acknowledged. This publication has received funding from the European Unions Horizon 2020 research innovation programme under grant agreement 732642 (ULPEC project). have decided to focus on time-accelerated tasks [START_REF] Friedmann | Demonstrating hybrid learning in a flexible neuromorphic hardware system[END_REF], which alleviates the need for implementing a slow paced learning rule. Subthreshold CMOS architectures may offer solutions to implement STDP-capable synapses that could handle reasonably long timescales, as suggested by the Dynamic Neuromorphic Asynchronous Processors (DYNAPs) architecture [START_REF] Moradi | A scalable multicore architecture with heterogeneous memory structures for dynamic neuromorphic asynchronous processors (DYNAPs)[END_REF], at the expense of silicon area. Proof of concept of an experimental small-scale hardware architecture with TiO 2 memristive synapses and on-line STDP-learning capabilities was recently reported by Prezioso et al. [START_REF] Prezioso | Spike-timing-dependent plasticity learning of coincidence detection with passively integrated memristive circuits[END_REF], with neurons that are however still implemented with discrete electronics and large capacitor values (∼ 100 nF). Highly integrated and lowpower memristor-based hardware implementations of spiking neural networks with on-line learning capabilities still remain a challenging goal.

The following work derives from studies made to design a demonstrator of a smart vision sensor (European project ULPEC [START_REF]Ultra-Low Power Event-Based Camera (ULPEC) H2020 project[END_REF]) that co-integrates an event-based retinomorphic vision sensor [START_REF] Posch | Retinomorphic event-based vision sensors: Bioinspired cameras with spiking output[END_REF] with a hardware spiking neural network implemented with analog CMOS neurons and a synaptic crossbar array of ferroelectric memristors. This paper makes several contributions in that context. First, we illustrate the need for a better match of the hardware time constants with the input event dynamics on a common event-based dataset (N-MNIST [START_REF] Orchard | Converting static image datasets to spiking neuromorphic datasets using saccades[END_REF]) and we introduce an accordingly modified STDP learning rule that requires little CMOS overhead (Section III). We then show the resilience of the system to two common variability issues of the postsynaptic circuits (Section IV). Finally, we study challenges posed by the initial synaptic weights on the learning and discuss two possible solutions (Section V).

II. SYSTEM ARCHITECTURE AND METHODOLOGY

In this work, we study the architecture shown in Fig. 1(a) using a Python in-house simulator that allows high-level description of the behavior of each hardware building block.

Each output Leaky Integrate-and-Fire (LIF) neuron includes a second generation current conveyor [START_REF] Lecerf | Silicon neuron dedicated to memristive spiking neural networks[END_REF], which allows to apply postsynaptic voltage pulses onto the memristive synapses (Fig. 1(b-c)). To save power the current conveyor only transmits current if i X > 0. The output current i Z scales i X by a factor K (constant in our simulations) in order to accommodate the limited value of the membrane capacitor C mem 1 pF, while allowing to overcome the discharge (∼ 100 V•s -1 ) due to the constant leakage current i discharge = 100 pA over the timescale of the targeted patterns. An output event is generated when the membrane voltage V mem reaches a threshold value (1 V in our simulations). Following a Winner-Take-All strategy, an output neuron that fires resets all the other output neurons. Besides, during training, a neuron that fired is reset until N refrac other output neurons have fired (in our simulations N refrac = 10), which prevents a single neuron from overtaking all the activity.

V pre1 V pre2 i X i Z forwarded current iz < 0 iz > 0 iz > 0 Vpre waveform during inference bipolar bipolar Vpre < 0 incomplete charge unnecessary discharge * Vmid = 1.1 V Vlow = 0.1 V Vhigh = 2.2 V V 0.
To prevent the shadowing of input events due to overlapping presynaptic pulses (Fig. 1(d)), we split the presynaptic pulse waveform into two components (Fig. 1(c)) that are applied independently. Inference pulses of duration T LTP = 10 µs are triggered by events received from the input visual sensor and can participate to synaptic long term potentiation. Shorter (0.5 µs) programming pulses can overlap with the lowamplitude part of a postsynaptic pulse for synaptic long term depression. From a circuit design viewpoint, the longer the time constants, the more chip area and energy is required.

If negative, the input voltage offset ε of the current conveyor (Fig. 1(b)) continuously charges the membrane capacitor, resulting in spurious output events. Thus, one must ensure a positive offset value, e.g., by adding a dedicated low-resolution digital-to-analog converter used to decrease the voltage offset at the network setup.

Many memristive technologies can have their conductance tuned by applying a voltage pulse beyond a threshold value [START_REF] Kuzum | Synaptic electronics: materials, devices and applications[END_REF]. Among them, ferroelectric memristors offer extremely low values of conductance [START_REF] Chanthbouala | A ferroelectric memristor[END_REF]. For simplicity sake and simulation speed, here we consider that the memristive devices are programmed accordingly to

∆G =      +A pot × (G max -G 0 ) for V syn ≤ -1.2 V -A dep × (G 0 -G min ) for V syn ≥ +1.2 V 0 otherwise , (1)
where G 0 , G min and G max are respectively the current, the minimal and the maximal conductance values, V syn is the voltage applied onto the memristor (Fig. 1(b)), and A pot and A dep are (dimensionless) learning rates. Eq. ( 1) is reminiscent of a multiplicative Spike-Timing Dependent Plasticity (STDP) learning rule implemented with rectangular programming pulses of constant amplitude and duration. The precise characteristics of the STDP should have only a limited impact on the insights provided by our study as several works have shown the high resilience of similar architectures to synaptic variability or STDP curves [START_REF] Querlioz | Immunity to device variations in a spiking neural network with memristive nanodevices[END_REF], [START_REF] Serrano-Gotarredona | STDP and STDP variations with memristors for spiking neuromorphic learning systems[END_REF]. Based on recent work on ferroelectric devices [START_REF] Boyn | Learning through ferroelectric domain dynamics in solid-state synapses[END_REF], we consider G min ranging from 1 nS to 10 nS and G max from 0.1 µS to 1 µS. We use A pot = A dep = 0.1, which allows single-epoch learning in our simulations.

III. WHEN ONE NEEDS A LONGER TERM INPUT MEMORY

In all our simulations, we replace the event-based visual sensor with the N-MNIST dataset, a conversion of the MNIST handwritten digit pictures to an event-based framework using a neuromorphic vision sensor [START_REF] Orchard | Converting static image datasets to spiking neuromorphic datasets using saccades[END_REF]. Each one of the 60,000 (10,000) training (test) samples was generated through 3 consistent 100 ms-saccades, resulting in 3 separate subclasses per digit class. Due to the limited amount of neurons in the hardware target, for each sample, we only use the data i) from the 28×28 central pixels ii) in the first saccade and iii) with ON polarity (i.e., increasing light events).

Fig. 2 illustrates a major issue of entirely implementing STDP through analog voltage overlapping. When an output neuron fires, because of the large mismatch between the analog hardware timescale (∼ T LTP = 10 µs) and the input event rate (∼ 7000 events per second on average), only a handful of synapses are potentiated, and all the other ones are depressed. In our simulations, this leads to a ∼ 0 % recognition rate.

To overcome this issue-unspecific to N-MNIST-we propose an STDP variation by introducing a digital control block (Fig. 1(a)) that allows tracking input activity and driving the programming events. To keep the CMOS overhead as low as possible, we use a 1-bit flag per presynaptic neuron that is set to ON when the neuron receives an input event. When a postsynaptic neuron fires, the synapses that are connected to it are then potentiated if their activity flag is ON and depressed otherwise. All presynaptic activity flags are then reset to OFF and we present the next digit sample. This form of STDP is thus entirely driven by postsynaptic events; each synapse connected to an output neuron that spikes is either potentiated or depressed, using separate pulses with constant duration and amplitude (translated into Eq. ( 1) first two lines in our simulations). This scheme is reminiscent of previous works in the literature [START_REF] Nessler | Bayesian computation emerges in generic cortical microcircuits through spike-timingdependent plasticity[END_REF], [START_REF] Querlioz | Bioinspired programming of memory devices for implementing an inference engine[END_REF]. However, the time window that triggers long term potentiation is no more constant but depends on the last previous postsynaptic event. This makes the hardware implementation simpler: the timer needed for each presynaptic neuron is replaced with logic blocks.

Similarly to other works in the field, we focus on the subset of digits 5, 6 and 9 to reduce simulation time [START_REF] Iyer | Unsupervised learning of event-based image recordings using spike-timing-dependent plasticity[END_REF]. One epoch thus corresponds to 17,288 (2859) training (test) samples.

After training (for one epoch) and before testing, we disable output neurons that fired less than 10 times. For each remaining output neuron, if at least 60 % of its last 10 events correspond to the same input class, we label it accordingly to the latter. Otherwise we disable the neuron. Empirically, this simple labeling strategy provides good agreement with labeling performed by a human, and the system achieves good selectivity (Fig. 3).

Without any device variability, the mean recognition rate is 92.4 % (Fig. 4(a)), close to the 93.68 % rate from the literature with a slightly different STDP and 400 output neurons [START_REF] Iyer | Unsupervised learning of event-based image recordings using spike-timing-dependent plasticity[END_REF].

IV. IMPACT OF POSTSYNAPTIC NEURON VARIABILITY

In a hardware implementation, some dispersion will affect the characteristics of the postsynaptic neurons, like the current conveyor offset voltages or the membrane capacitor values.

The wider the range of offset voltages ε, the larger the current i cancel , which makes it harder for postsynaptic neurons to spike. This results in significantly less active and selective output neurons after training. Cadence ® Virtuoso Monte Carlo simulations of our current conveyor CMOS design suggest that the offset voltages ε may mostly lie between 1 mV and 2 mV. However, simulations of an even worse scenario (uniform distribution of ε from 0 mV to 4 mV) show that the recognition rate remains at 87.3 % on average (solid black line in Fig. 4). Besides, adding a challenging uniform dispersion of the membrane capacitor C mem between ±20 % around 1 pF impacts only slightly more the performance (dashed gray line in Fig. 4). For example, the mean recognition rate drops to 86.9 % with the largest offset voltage dispersion.

The results in this section suggest that the mixed-signal architecture we propose should remain operational even with levels of postsynaptic circuit variability that are plausible in a hardware implementation, provided that the redundancy among output neurons is sufficient (Fig. 4(b)).

V. IMPACT OF INITIAL WEIGHT DISTRIBUTION

The mean conductance values of the patterns learned by the system determine a (rough) suitable range for postsynaptic parameters like the current conveyor factor K. Hence it is to be expected that the final system performance depends on the initial conductance distribution. This section illustrates how < R network >, defined as the average resistance of the N inputs memristors connected to each output neuron, evolves in various situations and which challenges arise.

We typically scale the current conveyor factor K respectively to the average of the expected final conductance maps <G final > target , in order to allow a reasonable triggering of the output neurons after a successful learning. Then, starting from a smaller <R network > (i.e., <G init >≥<G final > target ) allows output neurons to fire and to learn (thick blue line in Fig. 5(b), with a 92.3 % mean recognition rate after one epoch). Besides in such a situation Fig. 5(a) shows that a system with 10 times more resistive synapses (solid black line) behaves similarly to a less resistive learning-capable system (thick pink line), provided that K is made 10 times stronger; in both cases, the mean recognition rate is 92.3 %. If on the contrary <R network > is initially too large (i.e., <G init > < < <G final > target ), the output neurons have difficulties to spike, which results in unsuccessful one-epoch learning (dashed red line in Fig. 5(b)). A first approach to this challenge is to introduce an homeostasis process for each postsynaptic neuron, e.g., on K. For example, starting with conductance values similar to the previous failed attempt but with K=1/12.5 and decreasing it by 0.001 (down to 1/125) each time the neuron spikes allows to reach a 92.4 % mean recognition rate (solid orange line in Fig. 5(b)). Another solution is to explicitly initialize the synapses to decrease <R network >. As an example, for each output neuron, when initially arbitrarily programming half of the synapses to their maximum conductance value G max , the system achieves 92.1 % mean recognition rate (solid black line in Fig. 5(b)). Implementing the latter approach is likely to have a significantly smaller circuit-level overhead. Even the aforementioned crude memristor initialization allows to reach a good performance level; one could then simply leverage the programming blocks used for training to feed the digital control block with dedicated input events.

VI. CONCLUSION

To mitigate the short learning timescale of analog CMOS, we proposed adding a binary memory to each presynaptic neuron. By simulation means, we illustrated the performance of our concept on the N-MNIST dataset and its resilience to neuron variability. We also studied how the initial synaptic weights impact the system operation, and how crude initialization or neuron homeostasis both mitigate this issue.

These results provide insights to build hardware spiking neural networks able to learn tasks with real-world timescales.

Fig. 1 .

 1 Fig.1. Design overview. (a) System architecture: spiking neural network with N inputs = 784 and Noutputs= 100 analog CMOS pre-and postsynaptic neurons (numbered circles) fully connected through a single layer of memristive synapses with no selection device (purple arrows). An event-based vision sensor feeds the network while a digital control block supervises the analog neuron operation. (b) Circuit-level overview of one Leaky Integrate-and-Fire (LIF) postsynaptic neuron, which includes a second generation current conveyor (CCII) block that only forwards positive currents (i.e., i X > 0). In this sketch, two input pixels are active and drive their respective memristor (devices on the left) while the rest of the memristive crossbar can be modeled as a current source (i cancel ). For simplicity sake, we assume constant Vsyn= -1 V for active synapses. The control block on the right is actually part of the digital control block in (a). (c) The pre-(left) and postsynaptic (right) voltage waveforms applied by the neurons. (d) Comparing input voltage waveform strategies for inference and illustrating the challenges of overlapping presynaptic waveforms, depending on Vpre polarity and the polarity of the current forwarded by the current conveyor.

Fig. 2 .

 2 Fig. 2. For each digit class, statistics of the maximum amount of events (Nevents) that fall into a same long-term potentiation time window (T LTP ) over the first saccade of the 60,000 samples of the N-MNIST dataset. Gray large bars and black thin bars show the average value and the full range, respectively.

Fig. 3 .

 3 Fig. 3. Example of outcomes after one simulation run (Cmem = 1 pF, ε = 0 mV, K = 1/125, initial weights uniformly distributed between G min = 10 nS and Gmax = 1 µS). (a) Some conductance maps used after labeling. (b) Confusion matrix (on 2628 [218] correct [wrong] output events).

CFig. 4 .

 4 Fig. 4. Impact of the current conveyor offset voltage (ε), without (black solid lines) and with uniformly distributed (gray dashed lines) dispersion of the postsynaptic neuron membrane capacitor values (Cmem). Random offset voltages ε are uniformly drawn between 0 mV and the upper bound reported along the x-axis. Results without (with) membrane capacitor dispersion are averaged over 25 (125) one-epoch runs, combining 5 random sets of initial conductance values with 5 random sets of offset voltages (and with 5 random sets of membrane capacitors). The random sets are kept consistent between the simulations as much as possible. The other system parameters are the same as in Fig. 3. Symbols show the average values and the error bars show the full ranges. (a) The recognition rate, and (b) the fraction of output neurons (over Noutputs = 100 neurons) that are not used for classification after labeling.

GFig. 5 .

 5 Fig. 5. Evolution of < R network > the average resistance value seen by each output neuron. All results are averaged over five one-epoch simulation runs (ε = 0 mV and Cmem = 1 pF). (a) <R network > normalized by the current conveyor factor K, when the synaptic crossbar is on average more resistive after the learning than initially. The thick pink line shows results with initial conductance values G uniformly distributed between 10 nS and 1 µS and K = 1/125 (average recognition rate [avgRR]: 92.3 %), while the thin black line corresponds to a case with ten times more resistive synapses and a factor K scaled accordingly (K = 1/12.5, avgRR: 92.3 %). (b) Illustration of challenges that arise when the synapses are initially more resistive on average than the patterns that would be learned after successful training. Conductance values can evolve between G min = 1 nS and Gmax = 1 µS. If not mentioned, K=1/125. Solid thick blue line: a situation similar to (a), for reference purposes (avgRR: 92.3 %). All three other cases use initial conductance values that are uniformly distributed between G min and Gmax/10. Dashed red line: unchanged learning strategy (avgRR: 0 %). Solid orange line: applying some homeostasis onto K during the learning (starting from K = 1/12.5, avgRR: 92.4 %). Solid black line: half of the synapses are forcefully initialized to Gmax before the learning (avgRR: 92.1 %).