François Mercier
email: francois.mercier@meteo.fr

Selime Gürol

Pierre Jolivet

Yann Michel

Thibaut Montmerle

Block Krylov methods for accelerating ensembles of variational data assimilations

Keywords: block Krylov techniques, dual approach, ensemble assimilation, Lanczos algorithm, parallelism, quasi-geostrophic model, variational assimilation 1 VARIATIONAL ASSIMILATION AND THE EDA: FORMULATION 2.1 Deterministic variational assimilation

We consider the problem of efficiently solving ensembles of variational data assimilations in the context of numerical weather prediction. Running several assimilations notably allows us to initialize ensemble prediction systems and to more accurately represent background-error statistics, but is computationally expensive, limiting ensemble size. We propose a new class of algorithms for speeding up the minimization of the ensemble of data assimilations. It consists of using block Krylov methods to simultaneously perform the minimization for all members of the ensemble, instead of performing each minimization separately. We develop preconditioned block Krylov versions of the Full Orthogonal Method and of the Lanczos algorithm in both primal and dual space. The latter works in observation space that is usually of smaller dimension than the state space, thus giving an advantage in terms of memory requirements and computational cost. We describe and compare several parallelization strategies for speeding up the minimization and limiting the communications. These methods have been tested on a quasi-geostrophic system, consisting of a simplified atmospheric circulation model equipped with an ensemble of 3D-Var schemes tuned to mimic some features of a limited-area numerical weather prediction system. Experimentation shows that the number of iterations needed to converge is drastically reduced by the block Krylov approaches. We indicate that, while working in primal space does not save significant computational time, working in the dual space may reduce the computational time by a factor of 2 to 5 (depending on ensemble size) compared to standard Krylov methods, making our approach attractive for operational use.

KE Y WOR DS

INTRODUCTION

In numerical weather prediction (NWP), the initialization of the state requires determination of the best initial state of the atmosphere from all available information, in particular from a background state (a previous short-range forecast) and a set of observations. Most meteorological centres rely on variational data assimilation schemes that are efficient in the context of very large state space (10 6 to 10 9 ) and numerous observations (10 3 to 10 7 ). This is because variational schemes (just like ensemble methods) do not explicitly store the full covariance matrices [START_REF] Dee | On-line estimation of error covariance parameters for atmospheric data assimilation[END_REF]. A minimization algorithm finds the solution iteratively; it only requires matrix-vector products. Variational data assimilation has allowed for the direct assimilation of radiances from numerous satellites and this has been a source of huge improvement in the quality of operational forecasts [START_REF] Rabier | Overview of global data assimilation developments in numerical weather-prediction centres[END_REF].

Originally, variational assimilation schemes used climatological estimates of the background-error covariances [START_REF] Bannister | A review of forecast-error covariance statistics in atmospheric variational data assimilation. II: modelling the forecast-error covariance statistics[END_REF]. However, those estimates lacked the flow-dependency and the cycling of the statistics as in the Kalman filter. More recently, the covariances have been estimated from a current ensemble of forecasts, or by making a weighted (hybrid) combination of such ensemble-based covariances together with (modelled) climatological ones [START_REF] Hamill | A hybrid ensemble Kalman filter-3D variational analysis scheme[END_REF]. This is fully compatible with the variational framework [START_REF] Lorenc | The potential of the ensemble Kalman filter for NWP: a comparison with 4D-Var[END_REF][START_REF] Buehner | Ensemble-derived stationary and flow-dependent background-error covariances: evaluation in a quasi-operational NWP setting[END_REF]. In some implementations, this ensemble comes from a separate ensemble data assimilation scheme that is not variational, but rather based on some form of the Ensemble Kalman Filter (EnKF; [START_REF] Evensen | Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics[END_REF][START_REF] Houtekamer | A sequential ensemble Kalman filter for atmospheric data assimilation[END_REF] or on the Local Ensemble Transform Kalman Filter [START_REF] Hunt | Efficient data assimilation for spatiotemporal chaos: a local ensemble transform Kalman filter[END_REF].

In the stochastic EnKF, every ensemble member includes observation perturbations which are drawn from the observation-error covariance matrix, and model perturbations which are constructed to be representative of model errors. This scheme can be mimicked by running an ensemble of variational data assimilations (EDAs) which includes similar observation and model perturbations [START_REF] Fisher | Background-error covariance modelling[END_REF][START_REF] Žagar | Balanced tropical data assimilation based on a study of equatorial waves in ECMWF short-range forecast errors[END_REF][START_REF] Kucukkaraca | Use of analysis ensembles in estimating flow-dependent background error variances[END_REF]. This EDA can be used to initialize ensemble prediction systems as the EnKF, and is believed to more accurately represent the background-error statistics [START_REF] Berre | The representation of the analysis effect in three error simulation techniques[END_REF][START_REF] Bowler | The effect of improved ensemble covariances on hybrid variational data assimilation[END_REF]. For limited-area modelling (LAM) at Météo-France, a regional EDA based on the AROME model [START_REF] Seity | The AROME-France convective-scale operational model[END_REF] is being developed.

Performing EDA with many members in the ensemble is costly, as many variational assimilation systems have to be run. So EDAs are generally performed at a coarser resolution than the deterministic system and with a limited ensemble size. This implies spurious sampling noise in the estimation of the covariances, which can be mitigated with spatial filtering techniques [START_REF] Buehner | Spectral and spatial localization of background-error correlations for data assimilation[END_REF][START_REF] Berre | Filtering of background-error variances and correlations by local spatial averaging: a review[END_REF][START_REF] Ménétrier | Linear filtering of sample covariances for ensemble-based data assimilation. Part I: optimality criteria and application to variance filtering and covariance localization[END_REF]. However, larger ensemble sizes or more computationally efficient EDAs are clearly desirable. This is particularly true for the operational AROME-France which is run at high resolution (currently 1.3 km in its deterministic version and 3.8 km in its EDA version) over western Europe, implying large state space (of dimension 10 8 for the EDA). Reducing the cost of the minimization of the EDAs would certainly be very beneficial for operations.

Different methods have been proposed in order to speed up the minimizations in the EDA. [START_REF] Desroziers | Accelerating and parallelizing minimizations in ensemble and deterministic variational assimilations[END_REF] use the similarity of the minimization problems that have to be solved. They suggest building a better starting point for a new perturbed minimization using Lanczos vectors estimated through one or several previous perturbed minimizations. [START_REF] Lorenc | Computing an ensemble of variational data assimilations using its mean and perturbations[END_REF] propose the Mean-Pert method: the minimizations of the EDA are replaced by a minimization for the ensemble mean followed by minimizations for the perturbations from it that are accelerated by preconditioning and early stopping. Those two schemes do obtain gains in efficiency, but at the detriment of more sequential computations; a minimization with unperturbed data or for the ensemble mean has to be performed first, possibly degrading the total time to solution.

In this study we propose another complementary approach based on the use of block Krylov methods. Under an appropriate linearization hypothesis, the EDA results in an ensemble of linear systems sharing the same Hessian but with different (perturbed) right-hand sides. The standard implementation of EDA consists of applying a Krylov method such as the conjugate gradient algorithm [START_REF] Hestenes | Methods of conjugate gradients for solving linear systems[END_REF] to independently solve these systems. However, block Krylov methods [START_REF] O'leary | The block conjugate gradient algorithm and related methods[END_REF] are especially designed to solve this problem jointly, in one pass. They broaden the subspace in which an approximate solution is iteratively searched, and so can be able to speed up the convergence rate.

The purpose of this paper is to discuss, implement and test block Krylov approaches in the context of EDAs with large systems as encountered in NWP. More precisely, we will use in this paper a simplified but realistic numerical model, the quasi-geostrophic (QG) model of [START_REF] Fandry | A two-layer quasi-geostrophic model of summer trough formation in the Australian subtropical easterlies[END_REF], associated with an ensemble of 3D-Var that mimics some properties of the AROME EDA run at Météo-France. The QG EDA system used here will be run at a resolution implying a control vector of dimension 4 × 10 5 . The implementation will be made in the context of the Oriented Object Prediction System (OOPS; [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF][START_REF] Bonavita | A strategy for data assimilation[END_REF], a common C++ framework developed by ECMWF, Météo-France and their partners. This software allows an easy use of new assimilation schemes with different numerical models, including the QG model and AROME. [START_REF] Gratton | An observation-space formulation of variational assimilation using a restricted preconditioned conjugate gradient algorithm[END_REF] have introduced range-space variants of Krylov methods that work in observation space. This allows us to speed up the minimization when the number of observations is small compared to the dimension of the problem. Indeed, the inner products are less expensive to compute, and the algorithm also uses less memory to store the vectors. We propose in this paper to adapt this dual formulation to block Krylov methods. A particular focus will be given to the development of efficient parallelization strategies in the framework of OOPS.

Section 2 briefly introduces the classical variational formulation and the algorithms commonly used to minimize the cost function, as well as the EDA problem. Section 3 presents the block Krylov methods used to solve the EDA problem, both in primal and dual spaces, and details the algorithms implemented in this work. Then section 4 focuses on the implementation of these algorithms in the OOPS framework and on the associated parallelization strategies. Section 5 presents the results obtained with the QG model, both in terms of convergence and in terms of time consumption. Finally section 6 recalls the main conclusions of this paper and introduces future ongoing developments.

The following derivation directly applies to 3D state estimation by 3D-Var, 3DEnVar or hybrid schemes, which differ only with respect to the way the background-error covariance matrix is formed and which we do not detail. The 4D extension of those formulations, e.g. 4D-Var (Le Dimet and Talagrand, 1986;[START_REF] Courtier | A strategy of operational implementation of 4D-Var using an incremental approach[END_REF] and 4DEnVar [START_REF] Liu | An ensemble-based four-dimensional variational data assimilation scheme. Part I: technical formulation and preliminary test[END_REF][START_REF] Buehner | Four-dimensional ensemble-variational data assimilation for global deterministic weather prediction[END_REF][START_REF] Desroziers | 4DEnVar: link with 4D state formulation of variational assimilation and different possible implementations[END_REF] are also of great interest. Following [START_REF] Lorenc | Comparison of hybrid-4DEnVar and hybrid-4DVar data assimilation methods for global NWP[END_REF], the same derivation can be applied with states that become four-dimensional. However, in the 4D-Var case the propagation by the nonlinear model has to be accounted for within the four-dimensional observation operator. We will return to this specific case in section 2.4. The derivations of sections 2.1, 2.2, and 2.3 are valid only for 3D-Var, 3DEnVar and 4DEnVar (including hybrid versions).

Following the notations of [START_REF] Ide | Unified notation for data assimilation: operational, sequential and variational[END_REF], x denotes a model state, x b the background state, and x a the best estimate of the model state -the analysis. y o refers to the observation vector. In this paper, we denote by N the dimension of the problem, i.e. x ∈ R N , and N obs the number of observations, i.e. y o ∈ R N obs . B ∈ R N×N is the error covariance matrix of the background and R ∈ R N obs ×N obs is the error covariance matrix of the observations.  is the observation operator, allowing us to estimate an observation vector from a model state. Under these notations, the analysis x a is defined as the model state minimizing the cost function:

J(x) = 1 2 x -x b T B -1 x -x b + 1 2 {y o - (x)} T R -1 {y o - (x)}. (1)
This cost function is in general non-quadratic, due to the presence of the nonlinear observation operator . All schemes use an incremental approach, where the minimization is performed with respect to the deviations from the background state. Two different approaches can be taken to perform the (approximate) minimization of Equation ( 1). The first one is a quasi-Newton method [START_REF] Nocedal | Numerical Optimization[END_REF], as in the operational version of AROME or in the Met Office 4D-Var [START_REF] Lorenc | The Met Office global three-dimensional variational data assimilation scheme[END_REF]which both use the M1QN3 minimizer [START_REF] Gilbert | The Module M1QN3-Version 3.1[END_REF]. Another option, used for instance in the global 4D-Var of Météo-France and ECMWF, consists of applying an approximate Gauss-Newton method [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Gratton | Approximate Gauss-Newton methods for nonlinear least squares problems[END_REF]. This amounts to solving a sequence of quadratic problems (referred as the outer loop) where the operators are relinearized against the current estimate (the guess). Each quadratic problem is then solved with a Conjugate Gradient (CG) or Lanczos method (referred to as the inner loop).

For the sake of simplicity, and because we are only discussing the inner loop problem, we write the incremental linearized cost function as:

J( x) = 1 2 x T B -1 x + 1 2 (d -H x) T R -1 (d -H x) , (2)
in which x ∈ R N is an increment related to x and x b by x = x b + x, and H is a linearized version of the observation operator. The innovation vector d measures the gap between the observations and the corresponding values predicted by the model:

d = y o -(x b ).
In the simpler case described by Equation 2, the analyzed increment x a is determined by solving the linear system which allows the gradient to vanish:

B -1 + H T R -1 H x = H T R -1 d.
(3)

Krylov subspace methods can be used for numerically solving this linear system. They include the Full Orthogonalization Method (FOM), the Lanczos method, or the Conjugate Gradient (CG) method. All these algorithms produce mathematically equivalent iterates in exact arithmetic if the linear system matrix in Equation 3 is symmetric and positive definite (cf. Chapter 6 of [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF].

In order to improve the convergence of Krylov subspace methods, a preconditioner is desirable. An efficient preconditioner accelerates the convergence sufficiently such that the gain in convergence speed overcomes the extra cost of constructing and applying the preconditioner (Chapter 10 of [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF]. For variational data assimilation, a common choice to precondition Equation 3 is the background-error covariance matrix [START_REF] Bannister | A review of forecast-error covariance statistics in atmospheric variational data assimilation. II: modelling the forecast-error covariance statistics[END_REF]. In order to apply B as a preconditioner, there are two equivalent variants of the algorithms. The first one applies a split preconditioner through the square root of B and its adjoint to a Krylov method using the canonical inner product. The second performs a right B preconditioning to a Krylov method using the B inner product1 [START_REF] Derber | A global oceanic data assimilation system[END_REF][START_REF] Gürol | B-preconditioned minimization algorithms for variational data assimilation with the dual formulation[END_REF] and does not require the factorization of B.

Alternatively, an approximate solution of Equation 3 can be found by solving the dual problem [START_REF] Courtier | Dual formulation of four-dimensional variational assimilation[END_REF][START_REF] Gratton | Preconditioning and globalizing conjugate gradients in dual space for quadratically penalized nonlinear-least squares problems[END_REF] which performs the minimization in the observation space. The dual formulation may be preferable when the number of observations is smaller than the dimension of the control variable x since it reduces the computational cost and requires less memory. In the dual space, the Krylov subspace method is applied to the following linear system:

R + HBH T = d, (4)
with ∈ R N obs a control variable in the observation space related to x by x = BH T .

The Physical Space Assimilation System (PSAS; [START_REF] Courtier | Dual formulation of four-dimensional variational assimilation[END_REF][START_REF] Cohn | Assessing the effects of data selection with the DAO physical-space statistical analysis system[END_REF]) solves this linear system using split R -1 preconditioning and CG with canonical inner product. Unfortunately, PSAS has a non-monotonic behaviour in the reduction of the primal quadratic cost function Equation 2 [START_REF] El Akkraoui | Intercomparison of the primal and dual formulations of variational data assimilation[END_REF]. Therefore, a better alternative is to use the Restricted Preconditioned CG (RPCG) algorithm from [START_REF] Gratton | An observation-space formulation of variational assimilation using a restricted preconditioned conjugate gradient algorithm[END_REF]. Following [START_REF] Gürol | B-preconditioned minimization algorithms for variational data assimilation with the dual formulation[END_REF], RPCG is solving the left R -1 preconditioned dual system:

Id N obs + R -1 HBH T = R -1 d, (5)
using CG equipped with the (semi-definite) HBH T inner product, and where Id N obs is the identity matrix of R N obs ×N obs . RPCG generates the same iterates as those generated when solving the primal linear system Equation 3. Equivalent restricted Lanczos or FOM algorithms can be used as well [START_REF] Gratton | Range-space variants and inexact matrix-vector products in Krylov solvers for linear systems arising from inverse problems[END_REF]2012;[START_REF] Gürol | B-preconditioned minimization algorithms for variational data assimilation with the dual formulation[END_REF].

Ensemble of variational data assimilations (EDA)

Unlike the EnKF or LETKF, variational schemes do not provide an estimate of the background-and analysis-error covariances [START_REF] Fisher | Estimating the covariance matrices of analysis and forecast error in variational data assimilation[END_REF]. A consistent way of estimating the errors in the system is to perform an EDA with perturbed statistics, following a Monte Carlo approach. We will denote the ensemble size by m and use the subscript k to design the kth member. In the EDA, observation errors are simulated by adding random perturbations y o k drawn from R, and some model error effects are simulated by inflation of forecast perturbations [START_REF] Raynaud | Accounting for model error in the Météo-France ensemble data assimilation system[END_REF] or by perturbing physical tendencies of the model. These perturbations are propagated through the data assimilation cycling, which leads to analysis and background perturbations (denoted by x b k ) whose covariances are consistent with those of associated errors [START_REF] Belo-Pereira | The use of an ensemble approach to study the background-error covariances in a global NWP[END_REF]. The EDA, in most implementations, is thus the (parallel) minimization of m perturbed cost functions based on the deterministic cost function Equation 1 [START_REF] Desroziers | Accelerating and parallelizing minimizations in ensemble and deterministic variational assimilations[END_REF][START_REF] Lorenc | Computing an ensemble of variational data assimilations using its mean and perturbations[END_REF]:

J k (x k ) = 1 2 x k -x b -x b k T B -1 k x k -x b -x b k + 1 2 y o k + y o k - k (x k ) T × R -1 k y o k + y o k - k (x k ) . ( 6
)
Note that in Equation 6, we have so far retained nonlinear effects of variational data assimilation through the dependence of B, R and  on the ensemble member k. The quality control applied to the perturbed observations can potentially lead to their rejection for some members. Then the observation vector (and thus the observation operator and the observation-error covariance matrix) can be different between the members. Variational quality control also changes the weights in R, possibly in a different way between the members. Finally, the dependence of B with k may arise from the use of nonlinear balance equations or nonlinear humidity change of variable.

Linearized EDA

The EDA as formulated by Equation 6 includes nonlinear effects. As discussed by [START_REF] Lorenc | Computing an ensemble of variational data assimilations using its mean and perturbations[END_REF], the EDA describes the background-error statistics under a linear approximation that is common to the EnKF. There are several ways to formulate a linearized EDA that differ on the linearization state (or trajectory). [START_REF] Lorenc | Computing an ensemble of variational data assimilations using its mean and perturbations[END_REF] propose to first solve for the ensemble mean analysis, retaining all nonlinear effects, then solve for the ensemble perturbations in a linear context. We will see that the block Krylov approach is directly applicable in such a linearized EDA context. Another approach would be to linearize against the ensemble mean for the background, which is simpler and more efficient but possibly detrimental as the EDA is not recentred. We first plan to implement this method with the AROME EDA, which is linearized except for the quality control (through the first-guess check) and the nonlinearity of some observation operators (in particular due to the dependence on the vertical coordinate). We suppose that in the linearized EDA we have a common observation vector y o and statistics. To this end we suppose that there is no quality control applied to the observations after their perturbation (although one can remain before the perturbation). The observation operator H is linearized against a common state, the (background or analysis) ensemble mean. This linearization around a common state, as well as the absence of quality control on perturbed observations, are obviously simplifications as H is supposed to depend on the member index k. However in the 3D-Var application considered in this paper (QG model, section 5) there is no quality control in any case and  is independant of k, so that we will not further explore the implications of these assumptions. Nevertheless for real 3D-Var systems, we would verify the consistency of our solution by re-evaluating the nonlinear cost function value (using  k ) at the end of the minimization process. Then, for each member k, the cost functions have the same expression as in Equation 2. The corresponding linearized subproblems in incremental assimilation are then

J k ( x k ) = 1 2 x T k B -1 x k + 1 2 (d k -H x k ) T R -1 (d k -H x k ), (7) where x k = x k -x b -x b
k are the deviations from the perturbed backgrounds for each member and

d k = y o + y o k -  k (x b + x b k ) are the innovations with perturbed observations, using  k (x b + x b k + x k ) ≈  k (x b + x b k ) + H x k ).
So, in the primal space, the minimization of these m linearized subproblems results in m linear systems with the same Hessian 3):

B -1 + H T R -1 H (cf. Equation
B -1 + H T R -1 H x k = b k , (8) with right-hand sides b k = H T R -1 d k for k ∈ ⟦1; m⟧.
Similarly, in the dual space, we have m linear systems with the same Hessian R + HBH T (cf. Equation 4) for k ∈ ⟦1; m⟧:

R + HBH T k = d k . (9)
The standard EDA implementation, in primal or dual space, consists of running m independent minimizations in parallel. This does not take full advantage of the fact that, in the linearized context defined above, all the linear systems have the same Hessian (with perturbed right-hand sides). The purpose of this paper is to introduce block Krylov methods which exploit the mentioned structure of the problem to speed up the numerical minimization process.

Comments for 4D-Var application

The propagation of the state by a nonlinear model in the 4D-Var formulation makes the previous derivation more complex. First, we note that in this case, the nonlinear operator  of Equation 1 incorporates the propagation by the nonlinear model. Then, when using the Gauss-Newton approach to solve the incremental approach, the problem is relinearized around the guess through the outer loop process, which is not taken into account in Equation 2. There are additional nonlinear refinements to some 4D-Var schemes which we do not discuss here. They include the quality control of observations with non-Gaussian errors (Andersson and Jarvinen, 1999;[START_REF] Tavolato | On the use of a Huber norm for observation quality control in the ECMWF 4D-Var[END_REF], the variational bias correction [START_REF] Derber | The use of TOVS cloud-cleared radiances in the NCEP SSI analysis system[END_REF][START_REF] Auligné | Adaptive bias correction for satellite data in a numerical weather prediction system[END_REF], the dependence of background-error statistics on the first guess through nonlinear balance equations [START_REF] Fisher | Background-error covariance modelling[END_REF][START_REF] Barker | A three-dimensional variational data assimilation system for MM5: implementation and initial results[END_REF] or the nonlinear humidity change of variable [START_REF] Žagar | Balanced tropical data assimilation based on a study of equatorial waves in ECMWF short-range forecast errors[END_REF][START_REF] Ingleby | Improved variational analyses using a nonlinear humidity control variable[END_REF].

Moreover, for performing an EDA written in its linearized incremental form as in Equation 7, we would have to linearize the nonlinear operators  k of Equation 6(which include the nonlinear model trajectories) against a reference trajectory given by the ensemble mean background or analysis, and to perform a single outer loop. However, the influence of nonlinear effects in the EDA is outside the scope of this paper and deserves further research.

B L O CK KR Y L OV M E T H O DS FOR LINEARIZED EDA

Block Krylov methods in primal space

In this section, we introduce block Krylov methods for solving simultaneously the linearized EDA in primal space as described in section 2.2. General introductions to block Krylov techniques can be found in O' Leary (1980) and [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: an introduction[END_REF].

Notation / right-B preconditioning

The numerical experiments carried out in this paper apply a block Krylov subspace method equipped with B inner product to the right-B preconditioned linear system. For an arbitrary member k in the ensemble, it consists of rewriting Equation 8in the form

Id N + H T R -1 HB u k = b k , with Id N the identity matrix of R N×N , still b k = H T R -1 d k ,
and with the new transformed control variable u k ∈ R N linked to x through

x k = Bu k .
Then, we introduce the notation A = Id N + H T R -1 HB for the Hessian of this problem, so that it can be written

∀k ∈ ⟦1; m⟧, Au k = b k . ( 10
)
The block Krylov methods we describe here are iterative processes. So in this section, subscript p stands for the final iteration, while subscript i denotes the quantities corresponding to an arbitrary iteration. Moreover, as previously, subscript m is the total number of members while subscript k denotes an arbitrary member.

At final iteration p, an iterative process produces an approximate solution for u k which will be written u k,p . We note that we suppose here that the starting point of the minimization is u k,0 = 0. Then, the initial residuals for member k, denoted

r k,0 , is r k,0 = b k -Au k,0 = b k .
Under these assumptions, Equation 10 can be written in a block form:

Au = b, (11)
with the block control variable u ∈ R N×m :

u = u 1 • • • u m ,
and the block right-hand side b ∈ R N×m :

b = b 1 • • • b m .
The block initial residuals are

r 0 = r 1,0 • • • r m,0 = b.
At the end of an iterative process with p iterations, the approximate solution is u p ∈ R N×m , while kth member approximate solution is still u k,p .

Presentation

A block Krylov method is an iterative projection method which consists, at iteration p, in looking for an approximate

u p of u with u p ∈  p (A, r 0 ) , in which  p (A, r 0 ) is the block Krylov subspace defined by  p (A, r 0 ) = Blockspan r 0 , Ar 0 , … , A p-1 r 0 .
Here, "Blockspan" (Gutknecht, 2006, p. 9) is defined by

 p (A, r 0 ) = r 0 1 +Ar 0 2 +• • • + A p-1 r 0 p , with 1 , … , p ∈ R m×m . ( 12
)
We can write this as well for any single member u k,p :

u k,p ∈ Span A i r j,0 | | (i, j) ∈ ⟦0; p -1⟧ × ⟦1; m⟧ . ( 13
)
This highlights the fact that all members use the same research space, which is spanned from the entire ensemble.

u p ∈  p (A, r 0 ) means u p = V p s , ( 14
)
where V p is an orthonormal basis of  p (A, r 0 ), i.e. m⋅p vectors and V p ∈ R N×m⋅p if all the vectors implied in Equation 13are independent, and s ∈ R m⋅p×m the coefficients of a linear combination. We call the process allowing us to estimate V p the Arnoldi process, and V p is called the Arnoldi basis.

To define uniquely the approximate u p , we impose a Galerkin condition. Namely, we impose that r p , the block of residuals at iteration p, is orthogonal to the pth block Krylov subspace. This can be used to define the block FOM. If A is symmetric, it leads to the block Lanczos method which is mathematically equivalent to the block CG and the block FOM. The condition of orthogonality of the residuals

r p ⟂  p (A, r 0 ) can be written V T p r p = 0.
Using the definition of the block of residuals r p = b -Au p and the fact that b = r 0 ,

V T p r 0 = V T p Au p . ( 15
)
Combining Equations 14 and 15 leads to an expression for u p :

u p = V p T -1 V T p r 0 , (16)
with T = V T p AV p ∈ R m⋅p×m⋅p . T is constructed in the iterative process during the orthonormalization of the basis of the block Krylov subspace (section 3.1.3). Generally, in the NWP context, we have small p and m compared to the dimensions of the problem (m ⋅ p < 2 × 10 3 versus N = 10 8 and N obs = 10 5 ), so that the inversion of T is not challenging.

We can note that the search subspace at iteration 1 for any member is the space spanned by r 0 , so that we can write

r 0 = v 1 0 ,
where v 1 ∈ R N×m are the m first base vectors and 0 ∈ R m×m is a projection matrix resulting from a QR decomposition of r 0 .

Then, Equation 16 can be written

u p = V p T -1 ⎡ ⎢ ⎢ ⎣ v T 1 ⋮ v T p ⎤ ⎥ ⎥ ⎦ v 1 0 .
Using the orthogonality of the base vectors,

u p = V p T -1 e 1 0 , (17)
with e 1 = Id m 0 m • • • 0 m , where Id m is the identity matrix of R m×m and 0 m the null matrix of the same dimension. We remark that the search subspace of the block Krylov methods is the same for all the members. It combines information coming from the Hessian and from all the right-hand sides of the problem. Then its dimension is lower than or equal to p ⋅ m, with equality if all the vectors implied in Equation 13are independent. At each iteration, the dimension of the block Krylov subspace increases by (at most) m, the number of members.

If at some point of the iterative process we generate vectors which appear to be linked, the dimension of the search subspaces decreases; we call this an inexact breakdown. Such a breakdown may appear for example when an exact solution is found for one of the right-hand sides in the block (also known as happy breakdown in the context of non-block Krylov methods), or when there is a loss of orthogonality between residuals. Dealing with this phenomenon involves making an explicit reduction of the number of right-hand sides, for instance using deflation [START_REF] Gutknecht | Block Krylov space methods for linear systems with multiple right-hand sides: an introduction[END_REF]. Deflation may be used not only in the case of exact linear dependence of the generated vectors; techniques have been developed to force deflation by defining a tolerance threshold below which vectors are considered as linked, so that we can reduce the size of the Arnoldi basis.

However, in this paper, we deal with vectors of large dimensions (> 10 5 ) which are randomly perturbed (both backgrounds and observations); right-hand sides are unlikely to be dependent. Also, we perform a limited number of iterations (typically less than 50), such that the minimization is stopped quite early, i.e. before full convergence to machine precision. This may limit the occurrence of finding dependencies between the vectors through the algorithm. In all experiments presented in this paper, we have not encountered such cases of inexact breakdowns. Therefore, deflation has not been found to be necessary and we will not look into this issue in this paper. This point may be considered again with a more realistic system, although preliminary experiments confirm that we do not need deflation.

Block B-FOM

In our numerical experiments we use a block Krylov subspace method with the right-B preconditioner based on the Arnoldi process, the block B-FOM. First, we notice that the right-B preconditioned Hessian, i.e.

A = Id N + H T R -1 HB, is not symmetric with respect to the canonical inner product. Although block B-FOM can be used for unsymmetric systems, it is known that the iterative solution of a symmetric system is more stable. For this reason, we use the B-inner product with respect to which the preconditioned matrix is symmetric. The corresponding block B-FOM is given in Algorithm 1.

In Algorithm 1, V p and Z p , both in R N×m⋅p at iteration p, contain the orthonormal basis of the block Krylov subspace respectively in the physical space and in the B-preconditioned space, i.e. Z p = BV p . During an iteration (lines 5-17), we firstly define the new m search directions (w i ∈ R N×m , line 6). Then, lines 7-11 orthogonalize these w i against all the base vectors computed at previous iterations. The projection matrices j are stored as block components of matrix T (line 10). In line 14, matrix i represents the lower-diagonal block component of matrix T. Lines 18-19 allow us to calculate

Algorithm 1 Block B-FOM 1: z 0 = Br 0 2: v 1 , z 1 , 0 = qr (r 0 , z 0 ) 3: V 1 = [v 1 ] 4: Z 1 = [z 1 ] 5: for i = 1 ∶ p -1 do 6: w i = H T R -1 Hz i + v i 7: for j = 1 ∶ i do 8: j = z T j w i 9:
w i ∶= w iv j j 10:

T jm∶(j+1)m,im∶(i+1)m = j 11:
end for 12:

z i = Bw i 13: v i+1 , z i+1 , i = qr (w i , z i ) 14: T (i+1)m∶(i+2)m,im∶(i+1)m = i 15: V i+1 ∶= [V i ; v i+1 ] 16: Z i+1 ∶= [Z i ; z i+1
] 17: end for 18: solve T 1∶p⋅m,1∶p⋅m s p = e 1 0 19: x p = Z p s p the approximate solution at iteration p, corresponding to Equation 17, using matrix T 1∶m⋅p,1∶m⋅p ∈ R m⋅p×m⋅p . (We notice that x p = Bu p = BV p s p = Z p s p .)

We note that the quantities in the B-preconditioned space (z i and Z i ) are used as intermediates of calculations to facilitate the application of the B inner product and of the Hessian, so that any operator (B, H, H T , R -1 ) is applied only once per iteration.

We note as well that there is no criterion allowing us to stop the minimization in the case of satisfying convergence in Algorithm 1. (We suppose that we perform a predefined number of iterations.) Obviously we could add a criterion calculated at the end of each iteration to monitor the convergence, for instance on the residuals or on the background part of the cost function. However because many criteria exist, we leave this discussion for section 5.3.1 on our results.

The QR decomposition (lines 2 and 13) has to be made with respect to the B inner product as well (oblique QR decomposition). It consists of writing w i in the form

w i = v i+1 i , with v i+1 ∈ R N×m , i ∈ R m×m , v T
i+1 Bv i+1 = Id m , and i upper triangular. Several techniques exist to do so. The first method considered in this work is the Modified Gram-Schmidt QR algorithm (MGSQR). It is based on the fact that the columns of v i+1 form an orthonormal base with respect to the B inner product of the space spanned by the columns of w i . The modified orthogonalization process of the Gram-Schmidt method allows to estimate these base vectors, and the coefficients of the matrix i are estimated as the projection coefficients of the process; the Appendix gives details. Another algorithm is the Cholesky QR method (CholQR; [START_REF] Stathopoulos | A block orthogonalization procedure with constant synchronization requirements[END_REF]. It is based on the fact that if

w i = v i+1 i , we can write w T i Bw i = T i v T i+1 Bv i+1 i = T i i
, thanks to the orthogonality of v i+1 . Thus, i appears to be the Cholesky factor of w T i Bw i . Finally, solving the triangular system v i+1 i = w i gives v i+1 .

The last method, Pre-Cholesky QR (PRE-CHOLQR), algorithm 2 in [START_REF] Lowery | Stability analysis of QR factorization in an oblique inner product[END_REF], is a variant of the previous one in which we perform in a first step a standard QR decomposition with respect to the canonical inner product of w i . Then, CholQR is used to perform an oblique QR decomposition of the Q factor resulting from the first step. All these algorithms are described and evaluated in terms of stability in [START_REF] Lowery | Stability analysis of QR factorization in an oblique inner product[END_REF]: overall, these experiments on small systems and various conditioning both for B and w i show that PRE-CHOLQR is more stable than the other two, notably in terms of loss of orthogonality, but it is more costly as well. We have never found any serious orthogonality issue in our numerical experiments, nor have we found significant numerical differences between the three algorithms. So we finally decided not to use the more costly PRE-CHOLQR algorithm, and to work with MGSQR instead, detailed in the Appendix.

Block B-Lanczos algorithm

In the case in which the Hessian A is symmetric, the matrix T of the block B-FOM is symmetric (Saad, 2003, p. 185).

Using the notations of Algorithm 1, this leads to j = 0 for j ∉ {i -1, i} and i-1 = T i . Finally, this leads to a block tridiagonal matrix T with blocks of size m ⋅ m.

This results in a considerable simplification of the block B-FOM (the for loop at line 7 is replaced by the two cases j = i -1; j = i), and defines the block B-Lanczos algorithm.

However, it appears that, in finite precision because of the round-off errors, the orthogonality may not be ensured [START_REF] Gürol | B-preconditioned minimization algorithms for variational data assimilation with the dual formulation[END_REF]. It is therefore often necessary to add reorthogonalization in the block B-Lanczos procedure to ensure convergence (lines 8-9 of Algorithm 1 ∀j). Nevertheless, for systems in which p and m are small compared to the dimension N of the problem, the inversion of T in Algorithm 1 is not costly compared to the inner products. In this case, it is better to use the block B-FOM algorithm which does not restrict T to a block tridiagonal matrix.

3.2

Case m = 1

The purpose of this paper is to evaluate the efficiency of block Krylov methods for minimizing the cost function of an EDA. To do so, we will compare the results of such a m-member block minimization with results given by m independent non-block minimizations (m minimizations with a single member). As briefly described in section 2.1, Krylov methods, such as the FOM, can be used to perform a single member minimization. Actually, the block B-FOM in Algorithm 1 applied with m = 1 results in the B-FOM. In the numerical sections of this paper, we will compare the results given by a m-member block B-FOM with results given by m independent applications of the B-FOM.

Algorithm 2 B-FOM

1: z 0 = Br 0 2: 0 = √ (r 0 , z 0 ) 3: v 1 = r 0 ∕ 0 4: z 1 = z 0 ∕ 0 5: V 1 = [v 1 ] 6: Z 1 = [z 1 ] 7: for i = 1 ∶ p -1 do 8: w i = H T R -1 Hz i + v i 9: for j = 1 ∶ i do 10: T j,i = w i , z j 11: w i ∶= w i -T j,i v j 12:
end for 13:

z i = Bw i 14: T i+1,i = √ (v i , z i ) 15: v i+1 = w i ∕T i+1,i
16:

z i+1 = z i ∕T i+1,i 17: V i+1 ∶= [V i ; v i+1 ] 18: Z i+1 ∶= [Z i ; z i+1 ] 19: end for 20: solve T 1∶p,1∶p s p = 0 e 1 21: x p = V p s p
For readability purposes, we rewrite in Algorithm 2 the B-FOM and recall that it can be compared to Algorithm 1 when it is applied m times, potentially concurrently, with a single member.

In Algorithm 2, we now have w i , v i , and z i in R N and the bases V p and Z p in R N⋅p at iteration p. Matrices j and i of Algorithm 1 become scalars, directly stored in T. The QR decompositions result in simple normalizations.

First remarks on efficiency

It is known that the rate of convergence of Krylov and block Krylov methods is linked to the condition number of the Hessian and to the clustering of its eigenvalues. The condition number of a symmetric matrix is the ratio of the largest and the smallest eigenvalues. The greater the condition number, the slower the convergence. In NWP models, the unpreconditioned Hessian B -1 + H T R -1 H is generally very poorly conditioned (cf. table 1 in [START_REF] Desroziers | Accelerating and parallelizing minimizations in ensemble and deterministic variational assimilations[END_REF]. Using B as a preconditioner reduces the condition number, with the lowest eigenvalue equal to 1, though the condition number will be equal to the largest eigenvalue. This could allow the minimization to be speeded up. Additionally, it is shown on a simplified model in [START_REF] Haben | Conditioning of incremental variational data assimilation, with application to the Met Office system[END_REF] that the condition number of A is linearly linked to the number of observations available and to 2 b ∕ 2 o , with 2 b the variance of the background errors and 2 o the variance of the observations errors.

Concerning the differences between block B-FOM (Algorithm 1) and m times B-FOM (Algorithm 2), it can be observed that block B-FOM looks, for any member, for approximates of x in a search subspace of size m ⋅ p (number of members × number of iterations) whereas B-FOM looks for approximate solutions in search subspaces of size p. So by using a larger search subspace, we expect a faster convergence with block methods. Moreover, both algorithms imply one application of each operator (B, H, H T , and R -1 ) for each member and for each iteration (m applications by iteration). However, there are many more inner products and vector summations in block B-FOM than in m B-FOM. For instance, in block B-FOM, in the j loop at line 8 of Algorithm 1, we need m 2 inner products to calculate j , compared to m inner products for m B-FOM (line 10 of Algorithm 2 applied m times). Because we expect fewer iterations to converge with block B-FOM, we can already guess that the block approach will be much more efficient if the main part of the computational cost lies in the application of the operators (typically B or H) and less efficient if inner products and vector summations are very costly. Although this question will be addressed in detail in section 5, we can already note that this is a very good reason for using dual space algorithms. In the dual space, vectors of lower dimensions (N obs ≪ N) are used so that the inner products and vector summations will be less costly while the number of operators to be used will be the same. In the next section, block Krylov methods in dual space are introduced.

Block Krylov methods in dual space

Similarly to Equation 11 in primal space, we can write the ensemble problem in dual space defined by Equation 9 under the block form:

 = b, with  the Hessian in dual space with left R -1 preconditioning  = Id N obs + R -1 HBH T , = 1 • • • m , and b = b1 • • • bm , bk = R -1 d k , k ∈ ⟦1; m⟧.
Then, if we use a Krylov method equipped with the HBH T -inner product to solve this linear system, it will yield the same iterations as in the primal space. The Full Orthogonalization Method applied under these conditions gives the RB-FOM algorithm detailed in Algorithm 3. R stands for "Restricted" [START_REF] Gratton | Range-space variants and inexact matrix-vector products in Krylov solvers for linear systems arising from inverse problems[END_REF][START_REF] Gürol | B-preconditioned minimization algorithms for variational data assimilation with the dual formulation[END_REF].

IMPLEMENTATION

Oriented object prediction system

OOPS is a generic driving layer for NWP systems initiated by ECMWF and developed both at ECMWF and Météo-France [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF][START_REF] Arbogast | A parallel implementation of a 4DEnVar ensemble[END_REF][START_REF] Bonavita | A strategy for data assimilation[END_REF]. The purpose of OOPS is to develop a flexible code, written in the object-oriented language C++, to deal with the complex codes of NWP and oceanography.

It was notably developed to make easier the development Algorithm 3 Block RB-FOM 1: ẑ0 = HBH T r0 2: v1 , ẑ1 , 0 = qr r0 , ẑ0

3: V1 = [v 1 ] 4: Ẑ1 = [ẑ 1 ] 5: for i = 1 ∶ p -1 do 6: ŵi = R -1 ẑi + vi 7: for j = 1 ∶ i do 8: j = ẑT j ŵi 9:
ŵi ∶= ŵivj j 10:

T jm∶(j+1)m,im∶(i+1)m = j 11:

end for 12:

ẑi = HBH T ŵi 13:

vi+1 , ẑi+1 , i = qr ŵi , ẑi 14: [START_REF] Bonavita | A strategy for data assimilation[END_REF], and both NWP models of Météo France, ARPEGE [START_REF] Arbogast | A parallel implementation of a 4DEnVar ensemble[END_REF] and AROME, in their deterministic and EDA versions.

T (i+1)m∶(i+2)m,im∶ (
In this work, we have developed our block Krylov algorithms in the abstract layer of OOPS and have interfaced these algorithms with a FORTRAN model: the quasi-geostrophic (QG) numerical model [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF], cf. section 5, a simplified model of atmospheric circulation in midlatitudes. Distributed memory parallelism through the Message Passing Interface (MPI) is possible in the OOPS C++ layer. However, there can also be another MPI parallelization layer in the underlying FORTRAN code of the models (e.g. in AROME), even if it will not be the case for the QG model.

In this section, we introduce the three versions (mathematically equivalent) of the block Krylov algorithms that we have implemented in OOPS. The first one in section 4.2 does not include any parallelism in OOPS, while the two others in sections 4.3 and 4.4 do. For introducing these codes, we use as an illustration the block B-FOM described in Algorithm 1 (in primal space). However, the conclusions which will be drawn are also applicable for dual space block RB-FOM.

Sequential version

The first version of block B-FOM implemented in OOPS will be referred to as the "SEQ" version. This is a sequential implementation in the sense that there is no parallelization included in the OOPS layer. It can nevertheless include a MPI parallelization but only over the geographical domain (if such a parallelization is available). Figure 1a illustrates this situation for m = 2 and for four MPI processes on the Western Europe domain of AROME-France. Each process deals with small geographical areas but for all members in the ensemble. This implies for each process to perform more applications of operators and inner products/vector summations but with smaller vectors.

Then in the SEQ version, all the computations implied by Algorithm 1 are defined as loops on members. For instance, we detail the algorithm used to calculate line 6 of block B-FOM: application of the Hessian at iteration i to compute the new search directions.

for k = 1 ∶ m do (w i ) k = H T R -1 H (z i ) k + (v i ) k end for
And for line 8 of Algorithm 1, at iteration i:

for k = 1 ∶ m do for l = 1 ∶ m do j l,k = (w i ) k , z j l end

for end for

The next paragraph explains to what extent this SEQ approach could be competitive with a NWP model including an underlying geographical parallelization.

Suppose that we run EDA experiments with m members using N MPI MPI processes. Our interest is to compare the speed of convergence between concurrent independent B-FOM and block B-FOM. First, in the non-block case, each minimization for each member will be run on N MPI ∕m processes. So the geographical domain of the NWP model will be divided into N MPI ∕m parts. Figure 1b illustrates this situation with m = 2 and N MPI = 4. Processes 1 and 2 perform the minimization for the first member, while processes 3 and 4 perform the minimization for the second member. On the other hand, the SEQ block B-FOM algorithm has no OOPS parallelization, so that when it will be run on N MPI processes, the geographical domain will be divided into N MPI parts, giving a finer decomposition as is illustrated on Figure 1a with m = 2, N MPI = 4: the geographical domain is divided into four parts, and each process deals with both members.

Table 1 summarizes this analysis and compares the number of operator applications, inner products, and vector summations for non-block and SEQ block versions (columns 2 The SEQ method leads to m times more operator applications, and to roughly m 2 more inner products and vector summations. But the geographical decomposition is m times finer than in the non-block version, so that we can expect the same computational time for the operator applications and m more computational time for the inner products and vector summations. This is true if we neglect the communication time between geographical domains, which is not obvious. Knowing that we expect fewer iterations to converge with the block version, this approach could be competitive, especially if the operators are costly compared to the inner products.

If there is no parallelization in the underlying FORTRAN model, this SEQ version has no interest. All the MPI processes would run the same entire problem, and it should be faster to perform m non-block B-FOM concurrently. This is the case for the QG model described in section 5, so no results using the SEQ version will be studied in this paper. It has nevertheless been described here both as a simpler framework before further MPI developments and for future use in NWP models including geographical parallelization.

MPI version

Another version of block B-FOM, referred to as the "MPI" version, implies a workload distribution by member. This distribution can be combined with an underlying geographical decomposition. Figure 1b illustrates this situation, still with m = 2 and four MPI processes. Now each process deals only with one member, but on larger geographical areas. This reduces the number of operator applications and inner products/vector summations but they have to be performed with larger vectors. Under these conditions, line 6 of block B-FOM becomes very simple, without any loop. For example, for the process tied to member k:

(w i ) k = H T R -1 H (z i ) k + (v i ) k
This implies of course communication between processes at the OOPS level, for instance to compute the projections. Line 8 of Algorithm 1 becomes:

for l = 1 ∶ m do l ′ = (m -k + l) mod m if (l ′ ≠ k) then send z j k to proc l ′ get z j l ′ from proc l ′ j l ′ ,k = (z j ) l ′ , (w i ) k else j k,k = (z j ) k , (w i ) k end if end for
Here, for member k at iteration i, we want to compute the kth column of j , namely the projection of (w i ) k on all the z j l . The permutation implied by the definition of l ′ is the one used by [START_REF] Arbogast | A parallel implementation of a 4DEnVar ensemble[END_REF] For the number of increments exchanged in OOPS, we note that it is proportional to m 2 and j, the index of the current iteration. Indeed, we need to exchange all the previous base vectors to all processes at each iteration for the reorthogonalization part of block B-FOM (lines 7-11 of Algorithm 1). An idea for avoiding these communications would be to store the full Krylov bases corresponding to all members on each MPI process. This is the purpose of the next section.

MPI stored version

The third version of block B-FOM, referred to as the "MPI stored" version, is a variant of the MPI version which reduces the number of communications at the expense of increased memory. In this version, the full block Krylov base is stored redundantly on all processes. Computationally, it is nevertheless very similar but, at the end of each iteration, all processes send their new base vectors to all the others. Consequently, we do not need any communication in the reorthogonalization part of Algorithm 1 (lines 7-11), given that the full w i and z j are stored on each process. Line 8 of Algorithm 1 becomes:

for l = 1 ∶ m do j l,k = (z j ) l , (w i ) k end for
We can see in Table 1 that we avoid the dependence in j of the number of increments exchanged in OOPS by iteration.

The disadvantage of this approach is that we store all the Krylov bases m times. This implies storing m times objects of dimension 2 ⋅ j ⋅ m ⋅ N (two base vectors of dimension N per member and per iteration, though divided in reality in N MPI ∕m pools by the geographical decomposition). This can be intractable if m, j, and/or the dimension of the problem are large. So we can guess that this method will be mainly usable in the dual space approach in which smaller dimensions are in play. In the case that is considered in this paper, working in dual space with N obs = 2 × 10 4 , 40 iterations and 25 members, allows us to handle vectors of dimension (2 ⋅ j ⋅ m ⋅ N) upper bounded by 4 × 10 7 . Given that these experiments will be run on m supercomputer nodes, each with 64GB of memory, memory will not be an issue. For future applications with the AROME France NWP system, we expect to run ensembles of 75 members with around 10 5 assimilated observations. With 20 iterations, this implies storing vectors of dimension around 3 × 10 8 , which is still a tractable value given that 2-5 nodes by member are generally available to perform such experiments.

R E SU L T S O N T H E QG M OD E L

The two-layers QG model

The QG model [START_REF] Fandry | A two-layer quasi-geostrophic model of summer trough formation in the Australian subtropical easterlies[END_REF] allows to mimic QG flow in a cyclic channel with two layers. It is widely used as a simple way to realistically represent a large part of large-scale dynamics in the atmosphere in midlatitudes, notably to understand and to model cyclogenesis and storms (e.g. [START_REF] Plu | A cyclogenesis evolving into two distinct scenarios and its implications for short-term ensemble forecasting[END_REF]. A version of the QG model has been implemented in OOPS at ECMWF in order to test different assimilation strategies in a simple framework [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF]. We choose as well to test these new block Krylov developments on this model because it is available in OOPS and the developments made for the QG model should be easily portable to the OOPS version of AROME. Moreover, it is easy to tune for having realistic characteristics (see below).

It also allows us to test our code in a fully controlled and "perfect" framework: we have a well-posed problem with, for instance, the same operator H for all members. In this section, we briefly give the equations of the QG model. Readers can refer to [START_REF] Snyder | Linear evolution of error covariances in a quasigeostrophic model[END_REF] or [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF] for more details about this framework and its use in data assimilation.

The values of the parameters used in this study as well as the discretization strategy are detailed later in section 5.2. The QG model is based on the conservation of the potential vorticity. In the equations below, subscripts 1 and 2 denote respectively the top and bottom layers of the model. q is the potential vorticity, defined as:

q 1 = ΔΨ 1 -F 1 (Ψ 1 -Ψ 2 ) + y, q 2 = ΔΨ 2 -F 2 (Ψ 2 -Ψ 1 ) + y + R s , (18)
with Δ the Laplacian, Ψ the two-dimensional stream function, the northward derivative of the Coriolis parameter, R rms a coefficient representing the ground effects in the bottom layer (orography, heating), and F 1 and F 2 constant coefficients allowing the non-dimensionalization [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF] gives their detailed expressions), depending notably on the depths of the two layers, D 1 and D 2 .

The Lagrangian conservation of q can be written:

⎧ ⎪ ⎨ ⎪ ⎩ Dq 1 Dt = q 1 t + u 1 q 1 x + v 1 q 1 y = 0, Dq 2 Dt = q 2 t + u 2 q 2 x + v 2 q 2 y = 0, (19)
with u and v the zonal and meridional 2D wind fields, which can be derived from the stream function:

u i = Ψ i x , v i = Ψ i y . ( 20
)
The numerical resolution of this model is done with a pseudo-spectral method which involves the following sequence at each time step:

• Propagate q 1 and q 2 to the next time step using their conservation Equation 19, using a semi-Lagrangian scheme;

• Invert Equation 18to derive the stream function;

• Apply Equation 20 to get the horizontal wind fields in both layers.

This model is relevant for the large-scale dynamics of the atmosphere in terms of error growth, as it represents baroclinic instability [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF].

Experimental framework

Parametrization of the QG model

The domain on which these experiments are performed covers an horizontal area of 12,000 × 6,300 km 2 , discretized on a grid with 640 × 320 points. The spatial resolution is around 20 km in both horizontal directions. Then, the dimension of the problem will be N = 2 × 640 × 320 ≈ 4 × 10 5 . The true state and the true observations are generated using this model with depths of D 1 = 6 and D 2 = 4 km for the vertical layers. There are N obs = 1.2 × 10 4 observations of the stream function randomly distributed over the domain. Both the background and the observations are perturbed offline by appropriate codes in order to generate an ensemble. These perturbations are made according to the chosen statistics of errors. For the observation, R is supposed diagonal, with a constant standard deviation of o = 0.4. Concerning the background, the horizontal correlations are supposed to be Gaussian, isotropic, and homogeneous, with a length-scale of 1000 km, as in [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF]. The standard deviation is set to 1.6. The vertical correlation coefficient between the two layers is set to 0.2. We also add a model error by using in the assimilation step D 1 = 5.5 km and D 2 = 4.5 km for the depths of the layers.

The dimensions of this system have been chosen in order to be tractable but realistic to allowing us to be confident in the portability of the results to the AROME model. First, we note that, as has been mentioned in section 4, there is no underlying geographical workload distribution in the QG OOPS code. Consequently it is not possible to extent the QG model dimensions to values comparable to the AROME ones. That is why we limit ourselves to N = 4×10 5 , yielding a large but tractable system. Then, we remark that in dual space, the performance of the system will be driven both by N via the operators acting in primal space, e.g. B, and N obs . So we choose a value for N obs close to the one acting in the operational EDA of AROME (1.2×10 4 versus 10 4 -10 5 , depending on the weather conditions). It will nevertheless be inferior by an order of magnitude to the N obs used in the operational deterministic system or expected in a few years even for the EDA.

We have seen in section 3.1 that the speed of convergence of the block Krylov methods is mainly driven by the condition number of the problem and by the spectrum of the Hessian in general. The condition number is closely related to the covariances of the background and observation errors, as well as to the number of observations available. So it is quite important to use realistic values for all these parameters. In AROME, the number of assimilated observationsdepends highly on the weather situation. Notably, radars provide usable information when it is raining over France, so that heavy rain tends to significantly increase the condition number.

In AROME experiments carried on both for rainy and dry days, we have found that the condition number ranged from 2 × 10 2 up to 6 × 10 3 , corresponding to 1.4 × 10 4 to 4.9 × 10 4 assimilated observations. In the QG experiment described here, we have a condition number around 9 × 10 3 . We can reasonably believe that this matches the worst possible case of AROME. However, the relative dimensions of both problems produce eigenvalues decreasing much faster in our QG experiment than in any case of AROME. Eventually, what we can expect in the block Krylov methods applied to this QG problem is a convergence similar to the one expected with AROME in the first iterations of the minimization processes, but ultimately a faster convergence after a few iterations. So we have created a reasonably realistic model given the dimension constraints, but we still need to be careful before extrapolating the results to real systems.

Data assimilation experiments carried out

We finally describe the EDA experiments carried out with this QG model. First, it has to be noted that, in the QG experiments, we have found no numerical differences between all the algorithms described previously, namely the SEQ, MPI and MPIstored versions, in primal and dual space. Second, we recall that the QG model has no underlying geographical workload distribution, so the SEQ version has no computational interest. Thus we will compare the time performances of four algorithms: the MPI version of the block B-FOM (B-MPI) and of the dual block RB-FOM (RB-MPI), and the corresponding primal (B-MPIstored) and dual (RB-MPIstored) versions of the MPIstored algorithm. Regarding the number of members m considered, this will be set to 1, 5, 10, 20, or 40. All the experiments will use m MPI processes. This strategy allows to estimate the time needed to perform EDA minimization: for example, a one-member experiment run on one MPI process would need the same time to complete as 40 independent one-member experiments launched concurrently on 40 MPI processes.

Numerical results

In this section, the results for the QG EDA experiments are presented, first in terms of number of iterations, and then with respect to wall-clock time.

Convergence

We firstly study the number of iterations needed to reach convergence depending on the number of members in the ensemble. Figure 2 shows the convergence in terms of iterations according to three different criteria. We note again that the numerical results are the same for the four algorithms used here, so that there is no distinction between them. Figure 2a shows the B-norm of the residuals for the first member, i.e. at iteration p,

r p 1 = ‖ ‖ ‖ A x p 1 -b 1 ‖ ‖ ‖B
. We can see that the convergence speed, which is quite slow for the one-member case, increases on a regular basis when the number of members increases. Figure 2b shows the background part of the cost function J b which converges to a value of approximately 2.1 × 10 2 at the end of the minimization process. This final value is reached sooner and sooner when the number of members increases. We note that J b = ‖ x‖ 2 B ∕2. It shows that this criterion can be estimated through the minimization process, with the notations of Algorithm 1, as: Vs p T Zs p ∕2. Figure 2c shows the distance to its final value of the linearized cost function along the minimization processes: J lin -J * lin , with J * lin being the value of J lin at the end of the minimization. The same conclusion can be drawn: convergence is reached faster as the number of members increases. Now we choose a criterion for stopping the minimization. Several criteria can be used. The ECMWF system uses a criterion based on the relative variation of J b from an iteration to the next one: the minimization is stopped when it drops below 5%, but with three outer loops to ensure convergence to the nonlinear problem [START_REF] Trémolet | Incremental 4D-Var convergence study[END_REF]. In this paper, we choose a criterion based on the one currently acting operationally in deterministic AROME. In this system, the minimization is stopped after 40 iterations. So we set the stopping criterion at 40 iterations for the non-block one-member minimization, and we note the corresponding value of the B-norm of residuals (dotted horizontal grey line in Figure 2a). Then, block methods with various number of members are considered as converged when they reach the same value for the B-norm of residuals. This is obviously an a posteriori criterion, but it allows us to easily compare the non-block and block methods. Column 2 of Table 2 gives the corresponding number of iterations needed to converge for different numbers of members. We notice that the number of iterations is divided by more than four compared to the non-block method when using more than 20 members. We also note in Figure 2b

Wall-clock time

Increasing the number of members thus allows us to decrease the number of iterations needed to reach convergence, at the cost of more expensive iterations. In this section, the final cost of the minimization is studied in order to quantify the potential efficiency of block methods for speeding up the minimization process overall. We firstly explore the time needed to perform two individual iterations of the minimization, according to the algorithm used and the number of members in the ensemble. This is the purpose of Figure 3.

For the B-MPI algorithm in primal space, Figure 3a, we see that the computational time dramatically increases when the 3a) is to decrease the number of communications. While the time spent in the operators remains obviously the same, the time spent in other tasks significantly decreases. For instance, it goes from 12.5 s to 3 s for the 15th iteration of the minimization with 40 members. This proves that the communications are quite costly and that it is important to reduce their number. Nevertheless, the time needed to perform the 15th iteration is still eight times greater than the one needed to perform the same iteration with the non-block one-member version. Even if we perform far fewer iterations with forty members, it is not sure that this B-MPIstored algorithm will be able to improve the overall results. In order to decrease the cost of the inner products and vector summations, and to further decrease the cost of the communications, we have to work in the dual space.

The results for the dual space algorithm without base storage (RB-MPI) are presented in Figure 3b. Now, no matter the number of members, most of time is spent in the operators. Then, the time needed to perform iteration 15 with 40 members is no more than twice the time needed to perform the same iteration with one member. This is due to the fact that dual vectors have smaller dimension than primal ones: N obs = 1.2 × 10 4 versus N = 4 × 10 5 . This decreases dramatically the cost of the inner products, vector summations, and communications, since vectors of smaller dimensions are used. The cost of the operators remains the same.

The RB-MPIstored algorithm (RB-MPIst in Figure 3b) further reduces the number of communications compared to RB-MPI. With this approach, the time needed to compute an iteration is almost independent of the number of members and of the iteration index. For instance, 0.4 s are needed to compute iteration 15 with one member, and less than 0.6 s with 40 members. Since increasing the number of members allows us to decrease the number of iterations, we guess that this method will allow us to decrease the computational time needed to perform the full minimization.

The time spent in tasks other than operator applications for the best described algorithm, namely the RB-MPIstored, cf. grey part of Figure 3b, is decomposed in Figure 4 for iteration 15. These other tasks are divided in three groups: orthogonalization (Algorithm 3, lines 7-11 with j = i -1, i) and reorthogonalization (Algorithm 3, lines 7-11 with j < i -1) phases of RB-FOM (inner products and vector summations), and QR decomposition, which includes all the MPI communications in the MPIstored version of RB-FOM. We can see in Figure 4 that the reorthogonalization and the QR decomposition roughly account for the same part of the total extra-operator time consumption, whatever the number of members. This time consumption is proportional to the number of members for the reorthogonalization, which was expected (inner products and vector summations; Table 1). Now, we can look at the total time needed to perform the minimization process, given that this process ends when we reach the criterion defined in section 5.3.1, which corresponds to the number of iterations listed in Table 2. These results for the four algorithms are plotted in Figure 5, and summarized in columns 3-6 of Table 2. In Figure 5a corresponds to the primal algorithms while Figure 5b shows the dual ones. On each plot, bars show the total time needed to perform the minimization for different number of members. We still distinguish the time spent in the operators (dark) and in the other tasks (grey). This figure confirms the previous results. With the B-MPI algorithm, the best results are obtained with the non-block version (m = 1). The inner products and vector summations are too costly to have any advantage in using block methods, even if they require far fewer iterations. In block methods, the number of inner products and vector summations scales as m, and the number of communications as (m -1) 2 . This explains the increase with m in the time spent in other tasks in Figure 5a despite the decrease of the number of iterations that is performed. The reduction of the control vector dimension and of the communication number is thus crucial. Then if we store the bases with B-MPIstored (B-MPIst in Figure 5a) we obtain intermediate results. The performance of the algorithm is best with twenty members, but the gains are not really convincing. All methods imply a total computational time ranging from 10 to 20 s. The use of block methods starts to be attractive when using the dual space, as in algorithm RB-MPI (Figure 5b). We clearly see that increasing the number of members speeds up the minimization. We have even better results with RB-MPIstored (RB-MPIst in Figure 5b). Because most time is spent in operators, the time needed by the iteration is almost independent of the number of members. Then, the total computational time is around 3 s for 40 members compared to the 17 s of the non-block one-member version (Table 2). We also note in this table that the non-block version has performance almost independent of the algorithm used: 17 to 18 s in any case, almost completely spent in the operators.

CONCLUSION

In this paper, block Krylov methods have been considered for solving ensembles of variational data assimilations. In a linearized context, the EDA consists of an ensemble of linear systems sharing the same Hessian but with different right-hand sides. In this study, this implies linearization of the observation operator around the background ensemble mean and performing no quality control after the observation perturbation. Then Block Krylov methods are a natural solution to this problem. We propose to use the block Full Orthogonal Method, and derive variants both in primal space and in dual space with suitable preconditioning and convergence properties. We note that we limit ourselves to an assimilation system without outer loops. The implementation of these algorithms has been made under OOPS, in order to be easily applicable to different numerical weather and ocean prediction systems in the future. Several parallelization strategies have been developed to make the minimization tractable in terms of memory and efficient in terms of time consumption, relying on the infrastructure developed by [START_REF] Arbogast | A parallel implementation of a 4DEnVar ensemble[END_REF]. Finally, all these implementations have been tested on a QG numerical system with an ensemble of 3D-Var, a simplified but nevertheless realistic model. The dimension of its domain (4 × 10 5 ) was much smaller than is commonly used in NWP systems, but it has been run with a number of observations (1.2 × 10 4 ), consistent with the configuration of the EDA of the limited-area model of Météo-France, AROME.

On the one hand, these experiments on the QG system have shown that the block Krylov methods allow us to notably reduce the number of iterations needed to reach convergence (from 40 with one member to 14 with 10 members and 6 with 40 members). But on the other hand, block Krylov methods applied in primal space do not appear really convincing in terms of time consumption: the wall-clock times are greater or similar (depending on the parallelization strategy applied) to the ones obtained with independent non-block classical Krylov approaches. Nevertheless, performing block Krylov methods in dual space gives much more interesting results and considerably speeds up the minimization, so that the time needed to perform a single iteration is almost independent of the number of members in the ensemble. This is because we gain in the inner products and in the cost of communications. This approach allows us to gain in wall-clock time almost what is gained in iterations, reducing the time consumption from 17 s with one member down to 6.3 s with 10 members and 3.1 s with 40 members.

Moreover, because the dimension of the observation vectors is similar to the one in AROME, we can expect close results in dual space with the real system, even if these results need to be confirmed because of the existing differences in the Hessian eigenvalue spectra. We may need to further investigate the consequences of using a linear approximation in the EDA, notably the implication of the observation operator linearization around the background ensemble mean. Also, we note that the FOM algorithm may handle inexact matrix-vector products [START_REF] Gratton | Range-space variants and inexact matrix-vector products in Krylov solvers for linear systems arising from inverse problems[END_REF]. Therefore, we expect our block FOM algorithm also to perform with limited (single) precision, which may further reduce the computational cost. So exploring these block Krylov methods in a real NWP context is now the purpose of our ongoing work and already gives encouraging results which will be reported in the future.

FIGURE 1

 1 FIGURE 1 Illustration of the two workload distributions implemented in this work on the AROME-France domain for four MPI processes and two members in the ensemble. (a) Workload geographical distribution (no distribution by member), as used in the SEQ version of block B-FOM. (b) Workload distribution by member, combined with an underlying geographical distribution, as used in the MPI and MPIstored versions of block B-FOM, and in B-FOM

FIGURE 2

 2 FIGURE 2For different numbers of members in the ensemble (m from 1 to 40), convergence of the minimization for the first member in terms of iterations for different criteria: (a) B-norm of the residuals, (b) background part of the cost function, J b , and (c) distance to the final value of the linearized cost function, J lin

FIGURE 3

 3 FIGURE 3 Computational time for the second iteration (six left bars of each plot) and the 15th iteration (six right bars) of the minimization, according to the number of members (m) in the ensemble, for four different algorithms: (a) primal versions B-MPI and B-MPIstored, and (b) dual versions RB-MPI and RB-MPIstored. Dark shading denotes time spent in the operators, and light shading in other tasks (axpys, inner products, communications) [AUTHOR: The abbreviation 'axpys' has not been used elsewhere.]

FIGURE 4

 4 FIGURE 4For the RB-MPIstored algorithm at iteration 15, time spent in tasks which are not operator applications for different numbers of members (corresponding to the light shading part of the three right bars on Figure3d). These tasks are orthogonalization and reorthogonalization parts of RB-FOM (lines 7-11 of Algorithm 3) and QR decomposition and MPI communications (lines 13-16 of Algorithm 3)

FIGURE 5

 5 FIGURE 5As Figure3, but showing total computational time for performing the minimization process, according to the number of members in the ensemble.

B

  inner product used for performing the QR decomposition implied by block B-FOM Algorithm 1 and block RB-FOM Algorithm 3. MGSQR is given in Algorithm 4. Its inputs are two matrices w and z, both in R N×m , with z = Bw. The outputs are matrices v and z, again in R N×m , and in R m×m , which have to verify z = Bv, w = v (QR factorization), upper triangular and v T z = Id m (orthonormality with respect to the B inner product). In Algorithm 4, w p denotes the pth column of matrix w (same for all matrices of the same size). , z 0 )2: v 0 = w 0 ∕ 0,0 3: z 0 ∶= z 0 ∕ 0,0 4: for k = 1 ∶ mk , z k ) 12: v k ∶= v k ∕ k,k 13:z k ∶= z k ∕ k,k 14: end for

TABLE 1

 1 Characteristics of the different algorithms compared in this paper

	Non-block	SEQ	MPI	MPIstored

j =current iteration of the minimization process, m =number of members in the EDA, and k =current member

  to speed up the exchanges in a similar context. So this version implies many communi-Compared to the non-block version, even if there are more inner products and vector summations, the number of operator applications is the same, so that we can expect benefits if fewer iterations are needed to reach convergence and if it is costlier to apply the operators than computing inner products, which is expected at least in dual space. Compared to the SEQ version, this MPI version could be competitive only when dealing with extremely fine geographical decompositions.

cations and the OOPS decomposition by member implies a rougher geographical decomposition. While the geographical domain was divided into N MPI parts with the SEQ version, it is divided in N MPI ∕m parts in this MPI version, as in the non-block version (Figure

1b

). Table

1

presents the number of operator applications, inner products, and vector summations implied by this method, as well as the total number of communications per iteration.

TABLE 2

 2 For different number of members, m, the number of iterations needed to reach convergence (details in the text) and corresponding computation times (s) needed to perform the minimization process (columns 3-6) for the four algorithms tested here. For each algorithm, the best time performance is in bold MPI algorithm is not scalable and cannot compete against non-block methods. On the one hand, the time spent in the operators is quite constant with respect to the number of members and to the iteration index. This was expected since we perform in any case one application of each operator by member and by iteration. On the other hand, the time spent in inner products, vector summations, and communications soar with the number of members and the iteration index. The block minimization is thus intractable.The goal of the B-MPIstored algorithm (B-MPIst in Figure

		No. of	Primal			Dual
	m	iterations	B-MPI	B-MPIst	RB-MPI	RB-MPIst
	1	40	18.3	18.1	16.60	16.90
	5	22	22.8	13.6	9.89	9.60
	10	14	21.6	11.0	6.45	6.32
	20	9	20.7	10.4	4.36	4.27
	40	6	25.6	12.5	3.23	3.10
	number of members increases. It takes less than 0.5 s to per-
	form iteration 2 with one member, but around 3 s with forty
	members. With forty members, iteration 15 takes 12.5 s. The
	B-

∀(u, v) ∈ R N × R N , (u, v) B = v T Bu

ACKNOWLEDGEMENTS

This research has been supported by the AVENUE project funded by the RTRA-STAE foundation, Toulouse, France. The authors would also like to thank Etienne Arbogast for his help with OOPS.