M Fisher

S Gratton

S Gürol

Y Trémolet

X Vasseur

Low rank updates in preconditioning the saddle point systems arising from data assimilation problems

Keywords: data assimilation, limited memory preconditioning, saddle point system, weak-constraint

The numerical solution of saddle point systems has received a lot of attention over the past few years in a wide variety of applications such as constrained optimization, computational fluid dynamics and optimal control, to name a few. In this paper, we focus on the saddle point formulation of a large-scale variational data assimilation problem, where the computations involving the constraint blocks are supposed to be much more expensive than those related to the (1, 1) block of the saddle point matrix. New low-rank limited memory preconditioners exploiting the particular structure of the problem are proposed and analysed theoretically. Numerical experiments performed within the Object-Oriented Prediction System are presented to highlight the relevance of the proposed preconditioners.

Introduction

This paper investigates the saddle point representation of a large-scale quadratic minimization problem subject to linear equality constraints which is solved within a Gauss-Newton (GN) method. Therefore, in this study we deal with the solution of a sequence of saddle point systems particularly arising from data assimilation problems such as applications in Earth-system modelling. For instance, in meteorological applications the result of the data assimilation problem is the initial state of a dynamical system, which is then integrated forward in time to produce a weather forecast.

It has been shown that the saddle point formulation of the four-dimensional variational data assimilation (4D-Var) problem allows parallelization in the time dimension [START_REF] Rabier | Data assimilation in meteorology[END_REF]. Therefore, this formulation represents a crucial step towards improved parallel scalability of the 4D-Var problem.

When the saddle point system is large, the problem is solved iteratively, usually by using a Krylov subspace method. When using Krylov methods, it is well-known that their convergence can be slow on indefinite systems unless a good preconditioner is used to accelerate the convergence. Therefore, finding an efficient preconditioner is crucial when solving the large-scale saddle point systems [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]. *Corresponding author. Email: gurol@cerfacs.fr

In recent years, a large amount of work has been devoted to find efficient preconditioners such as block diagonal, triangular and constrained preconditioners (see, e.g. [START_REF] Benzi | Numerical solution of saddle point problems[END_REF][START_REF] Bergamaschi | Inexact constraint preconditioners for linear systems arising in interior point methods[END_REF][START_REF] Bramble | A preconditioning technique for indefinite systems resulting from mixed approximations of elliptic problems[END_REF][START_REF] Keller | Constraint preconditioning for indefinite linear systems[END_REF][START_REF] Kuznetsov | Efficient iterative solvers for elliptic finite element problems on nonmatching grids[END_REF][START_REF] Murphy | A note on preconditioning for indefinite linear systems[END_REF] and references therein). However, the efficiency of most of these preconditioners depends on the property that the computations involving the (1,1) primal variable block are computationally more expensive than the calculations involving the off-diagonal (constraint) blocks. This property does not hold for the saddle point representation of the variational data assimilation problem as we shall see. Therefore, many preconditioners in the literature become computationally expensive and not well-suited to our problem.

One preconditioner of interest for this study is the inexact constraint preconditioner proposed in [START_REF] Bergamaschi | Inexact constraint preconditioners for linear systems arising in interior point methods[END_REF]. It has been shown that this preconditioner can be used as a first-level preconditioner [START_REF] Rabier | Data assimilation in meteorology[END_REF]. In this study, we investigate whether it is possible to further accelerate the convergence by improving the first-level preconditioner with low-rank updates generated from the inherited information.

The paper is structured as follows. In Section 2, we derive the variational data assimilation problem in which model errors are taken into account. We then focus on solution algorithms. The first standard solution approach includes the model constraints into the objective function and solves the resulting unconstrained optimization problem. By using this approach, the expensive matrix-vector products need to be performed sequentially in the time domain. As an alternative to the standard approach in terms of parallelization in the time domain, the second approach introduces additional constraints on the observation operator which leads to a quadratic problem with nonlinear equality constraints. This formulation is then solved by introducing a sequence of quadratic programs where second order derivatives of the constraints are ignored as in GN-like methods. The main computational work in this algorithm consists in solving saddle point linear system as discussed in [START_REF] Rabier | Data assimilation in meteorology[END_REF].

In Section 3, we present new low-rank limited memory updates to precondition the saddle point matrix. These updates are obtained by solving the two-sided secant equations. According to our knowledge, two-sided equations are considered only in [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF][START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF] where authors used twosided secant equations in order to approximate the Jacobian matrix. In this study, we use the approximation for preconditioning, which we believe is new in the literature. The performance of the derived preconditioners is shown in Section 4 by using a simple two-layer quasi-geostrophic model, and finally the conclusions are summarized in Section 5.

Problem definition and solution algorithms

We consider a set of states {x i ∈ R n } valid at times {t i }, i = 0, 1, . . . , N during an analysis window [t 0 , t N ]. We define the model error {q i }, i = 1, 2, . . . , N associated with these states such that

x i = M i (x i-1 ) + q i , i ≥ 1,
where M i represents the physical model integration from time t i-1 to time t i . We assume that a set of observations y i ∈ R m i is available at time t i , i = 0, . . . , N. Our aim is then to find the best model fit to the observed data in the sense that the sum of square errors between the observed data and the model prediction is minimized. During this minimization, the model error will also be taken into account. This leads to a nonlinear least squares problem min x f (x) with

f (x) = 1 2 x 0 -x b 2 B -1 b + 1 2 N i=0 H i (x i ) -y i 2 R -1 i + 1 2 N i=1 q i 2 Q -1 i , ( 1 )
where x = [x 0 , x 1 , . . . , x N ] T is a four-dimensional variable (i.e. time-distributed state variable), x b ∈ R n is the a priori information (background) for the state x 0 , H i is an observation operator that represents a mapping from the state space to the observation space, B b , R i and Q i are n × n, m i × m i and n × n symmetric positive definite error covariance matrices corresponding to the background, observation and model errors, respectively. The first term in this nonlinear least squares problem is the Tikhonov regularization term [5, p. 101] which is often needed when the problem is under-determined, as it is the case in practical applications when there are fewer observations than variables in the model. The formulation of the weighted nonlinear cost function ( 1) is known in the data assimilation community as the weak-constraint four-dimensional variational data assimilation (weak-4D-Var) problem [START_REF] Trémolet | Accounting for an imperfect model in 4D-Var[END_REF][START_REF] Yoshikazu | Numerical variational analysis with weak constraint and application to surface analysis of severe storm gust[END_REF]. This approach has its origin in the maximum likelihood estimation under the assumptions that all errors are unbiased, uncorrelated, and can be represented by zero mean Gaussian distributions [START_REF] Trémolet | Accounting for an imperfect model in 4D-Var[END_REF]. In the nonlinear cost function [START_REF] Benzi | Numerical solution of saddle point problems[END_REF], the model and observation errors are also assumed to be temporally uncorrelated. Time correlations can be included at the cost of using block-structured covariance matrices, provided that the required statistics is known [START_REF] Trémolet | Model-error estimation in 4D-Var, Q[END_REF]. We also note that, in practical applications, the Gaussian assumption may not be always realistic and can be partially relaxed [START_REF] Hólm | Assimilation and modelling of the hydrological cycle: ECMWF's status and plans[END_REF][START_REF] Ingleby | Improved variational analyses using a nonlinear humidity control variable[END_REF][START_REF] Tavolato | On the use of a Huber norm for observation quality control in the ECMWF 4D-Var[END_REF]. The salient properties of this 4D-Var problem in real-life applications such as ocean or atmospheric data assimilation can be summarized as follows

• The dimension of the state variable n exceeds 10 6 in general. Hence, this problem can be really considered as a large-scale nonlinear least squares problem. • The integrations of the physical model involving M i are typically very costly. For this reason, exact second order derivative information is too expensive to compute. • The error covariance matrices are not explicitly available. They are modelled or estimated using either smoothing functions or filters [START_REF] Courtier | The ECMWF implementation of three-dimensional variational assimilation (3D-Var). i: Formulation[END_REF][START_REF] Weaver | Correlation modelling on the sphere using a generalized diffusion equation[END_REF], and only their actions on a vector can be computed.

Taking into consideration these main properties, we review in Sections 2.1 and 2.2 practical algorithms for the solution of the large-scale nonlinear least squares problem given by (1).

Solution via the GN method

A widely used algorithm to minimize the large-scale nonlinear least squares cost function [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] is the GN method known in the data assimilation community as Incremental 4D-Var [START_REF] Courtier | A strategy for operational implementation of 4D-Var using an incremental approach[END_REF]. At each iteration (say the jth iteration) of the GN method, named as the outer loop, a step (an increment) δx (j) from a given x (j) is computed by minimizing a quadratic cost function which is the linearized least squares approximation of the nonlinear problem. This quadratic approximation in the neighbourhood of the iterate x (j) is given by

J (j) = 1 2 δx 0 -r (j) 2 B -1 b + 1 2 N i=0 H (j) i δx i -d (j) i 2 R -1 i + 1 2 N i=1 δq i -c (j) i 2 Q -1 i , (2)
where

r (j) = x b -x (j) 0 , c (j) i = -q (j) i , d (j) i = y i -H i (x (j) i ).
In the quadratic cost function (2), H

(j)

i represents the Jacobian matrix of the observation operator H i at x (j) i , and δq i is defined by

δx i = M (j) i δx i-1 + δq i , where M (j) i represents the Jacobian matrix of the physical model M i at x (j) i-1 .
The state is then updated according to x (j+1) = x (j) + δx where δx = [δx 0 , δx 1 , . . . , δx N ] T . We can rewrite the quadratic subproblem (2) in a more compact form as

J = 1 2 δp -b 2 D -1 + HFδp -d 2 R -1 , ( 3 )
where we have dropped the outer loop index j for sake of clarity. In Equation ( 3), given

= n × (N + 1) and m = N i=0 m i , the vectors δp ∈ R , d ∈ R m and b ∈ R are defined as δp = ⎛ ⎜ ⎜ ⎜ ⎝ δx 0 δq 1 . . . δq N ⎞ ⎟ ⎟ ⎟ ⎠ , d = ⎛ ⎜ ⎜ ⎜ ⎝ d 0 d 1 . . . d N ⎞ ⎟ ⎟ ⎟ ⎠ , b = ⎛ ⎜ ⎜ ⎜ ⎝ r c 1 . . . c N ⎞ ⎟ ⎟ ⎟ ⎠ ,
and the matrices

D ∈ R × , R ∈ R m×m , H ∈ R m× are given by D = diag(B b , Q 1 , . . . , Q N ), R = diag(R 0 , . . . , R N ), H = diag(H 0 , . . . , H N ).
The vectors δx and δp are related by δx = Fδp with F a by matrix defined as

F = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I n M 1,1 I n M 1,2 M 2,2 I n . . . . . . . . . . . . M 1,N M 2,N • • • M N,N I n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
, where M i, ī = M¯i • • • M i represents an integration of the linear model from time t i-1 to time t¯i. The entries of this matrix are not available explicitly and the computation of matrix-vector products is performed through an operator. We stress that these matrix-vector products are very expensive since they require expensive model integrations.

Minimizing J defined in Equation (3) using iterative methods involves the computation of matrix-vector products with F, which requires sequential model integrations. Hence, since the matrix-vector products with M i, ī are known to be expensive, the sequential integrations will be very costly in this approach. An alternative method referred to as the 'saddle point approach' in [START_REF] Rabier | Data assimilation in meteorology[END_REF] is briefly explained next, since it will be the basis for further developments.

Solution via a sequence of quadratic programming problems

The saddle point formulation of the weak-constrained 4D-Var provides a framework for allowing parallelization in the time domain. This is a key issue when tackling large-scale problems in data assimilation on modern architectures. In this section, we briefly present this algorithm and refer the reader to [START_REF] Rabier | Data assimilation in meteorology[END_REF] for a complete description.

We can simply reformulate the nonlinear problem (1) as a the minimum of a convex quadratic function under nonlinear equality constraints as

f (p, w) = 1 2 p 2 D -1 + 1 2 w -y 2 R -1 (4)
subject to

p i = x i -M i (x i-1 ), (i = 1, . . . , N), w i = H i (x i ), (i = 0, 1, . . . , N),
where p = [x 0x b , p 1 , . . . , p N ] T , w = [w 0 , w 1 , . . . , w N ] T , and y = [y 0 , y 1 , . . . , y N ] T . This constrained minimization problem can be solved by using a sequential quadratic programming that represents a mapping from the state space to the observation space, B b , R i and Q i are n × n, m i × m i and n × n symmetric positive definite error covariance matrices corresponding to the background, observation and model errors, respectively. The first term in this nonlinear least squares problem is the Tikhonov regularization term [5, p. 101] which is often needed when the problem is under-determined, as it is the case in practical applications when there are fewer observations than variables in the model. The formulation of the weighted nonlinear cost function ( 1) is known in the data assimilation community as the weak-constraint four-dimensional variational data assimilation (weak-4D-Var) problem [START_REF] Trémolet | Accounting for an imperfect model in 4D-Var[END_REF][START_REF] Yoshikazu | Numerical variational analysis with weak constraint and application to surface analysis of severe storm gust[END_REF]. This approach has its origin in the maximum likelihood estimation under the assumptions that all errors are unbiased, uncorrelated, and can be represented by zero mean Gaussian distributions [START_REF] Trémolet | Accounting for an imperfect model in 4D-Var[END_REF]. In the nonlinear cost function [START_REF] Benzi | Numerical solution of saddle point problems[END_REF], the model and observation errors are also assumed to be temporally uncorrelated. Time correlations can be included at the cost of using block-structured covariance matrices, provided that the required statistics is known [START_REF] Trémolet | Model-error estimation in 4D-Var, Q[END_REF]. We also note that, in practical applications, the Gaussian assumption may not be always realistic and can be partially relaxed [START_REF] Hólm | Assimilation and modelling of the hydrological cycle: ECMWF's status and plans[END_REF][START_REF] Ingleby | Improved variational analyses using a nonlinear humidity control variable[END_REF][START_REF] Tavolato | On the use of a Huber norm for observation quality control in the ECMWF 4D-Var[END_REF]. The salient properties of this 4D-Var problem in real-life applications such as ocean or atmospheric data assimilation can be summarized as follows

• The dimension of the state variable n exceeds 10 6 in general. Hence, this problem can be really considered as a large-scale nonlinear least squares problem. • The integrations of the physical model involving M i are typically very costly. For this reason, exact second order derivative information is too expensive to compute. • The error covariance matrices are not explicitly available. They are modelled or estimated using either smoothing functions or filters [START_REF] Courtier | The ECMWF implementation of three-dimensional variational assimilation (3D-Var). i: Formulation[END_REF][START_REF] Weaver | Correlation modelling on the sphere using a generalized diffusion equation[END_REF], and only their actions on a vector can be computed.

Taking into consideration these main properties, we review in Sections 2.1 and 2.2 practical algorithms for the solution of the large-scale nonlinear least squares problem given by (1).

Solution via the GN method

A widely used algorithm to minimize the large-scale nonlinear least squares cost function [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] is the GN method known in the data assimilation community as Incremental 4D-Var [START_REF] Courtier | A strategy for operational implementation of 4D-Var using an incremental approach[END_REF]. At each iteration (say the jth iteration) of the GN method, named as the outer loop, a step (an increment) δx (j) from a given x (j) is computed by minimizing a quadratic cost function which is the linearized least squares approximation of the nonlinear problem. This quadratic approximation in the neighbourhood of the iterate x (j) is given by

J (j) = 1 2 δx 0 -r (j) 2 B -1 b + 1 2 N i=0 H (j) i δx i -d (j) i 2 R -1 i + 1 2 N i=1 δq i -c (j) i 2 Q -1 i , (2)
where

r (j) = x b -x (j) 0 , c (j) i = -q (j) i , d (j) i = y i -H i (x (j) i ).
In the quadratic cost function ( 2), H

(j)

i represents the Jacobian matrix of the observation operator H i at x (j) i , and δq i is defined by

δx i = M (j) i δx i-1 + δq i , where M (j) i represents the Jacobian matrix of the physical model M i at x (j) i-1 .
The state is then updated according to x (j+1) = x (j) + δx where δx = [δx 0 , δx 1 , . . . , δx N ] T . We can rewrite the quadratic subproblem (2) in a more compact As mentioned in [START_REF] Fisher | Data assimilation in weather forecasting: A case study in PDE-constrained optimization[END_REF], this approach is a particular case of the simultaneous analysis and design method (SAND) [START_REF] Biegler | DAE-constrained optimization[END_REF]. We note that the SAND formulation differs from the one introduced in [START_REF] Fisher | Data assimilation in weather forecasting: A case study in PDE-constrained optimization[END_REF] by the presence of an additional constraint in the observation space. The crucial property of the SAND approach is that it already offers a first level of parallelism in the time domain. Indeed, matrix-vector products with

F -1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I n -M 1 I n -M 2 I n . . . . . . -M N I n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
can be obtained by performing the integrations of the linear model in parallel. Since the model integrations are known to be expensive, this approach represents a crucial step towards higher computational efficiency. In the SAND formulation [START_REF] Björck | Numerical Methods for Least Squares Problems[END_REF], by defining additional constraints on the observation operator, we introduce additional parallelism in H and H T operators. Nevertheless, the dimension of the global problem (3 + 2m) has significantly increased. Hence, we consider the application of Schur complement reduction techniques to find the solution in a space of reduced dimension. Using the following partitioning of unknowns in ( 6)

⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ D -1 0 0 I 0 0 R -1 0 0 I 0 0 0 -F -T -H T I 0 -F -1 0 0 0 I -H 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ δp δw λ μ δx ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ = ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ D -1 b R -1 d 0 0 0 ⎞ ⎟ ⎟ ⎟ ⎟ ⎠
, a reduced system of order (2 + m) involving the negative Schur complement is obtained as

⎛ ⎝ D 0 F -1 0 R H F -T H T 0 ⎞ ⎠ ⎛ ⎝ λ μ δx ⎞ ⎠ = ⎛ ⎝ b d 0 ⎞ ⎠ . ( 7
)
The order of this reduced system is still large. Nevertheless, this formulation preserves the parallelization in the time domain. As a summary, the solution of the minimization problem ( 5) can be obtained by solving the saddle point system [START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF], which is then used to update the state vector x. In the remainder of this paper, we focus on efficient preconditioned iterative methods for the solution of the saddle point system given by (7).

Solution of the linearized subproblem.

We consider preconditioned Krylov subspace methods for the solution of ( 7) and refer the reader to [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] for a survey of numerical methods related to the solution of indefinite linear systems. First, we point out an important feature of our problem that we need to take into consideration when choosing an appropriate preconditioner. We write the saddle point matrix in [START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF] in a compact form as

A = ⎛ ⎝ D 0 F -1 0 R H F -T H T 0 ⎞ ⎠ = A B T B 0 , ( 8
)
with A ∈ R ( +m)×( +m) and B ∈ R ×( +m) , respectively. The computational costs related to the matrix-vector products with A in Equation ( 8) differ from most of the situations studied in the literature, in the sense that performing matrix-vectors with B is expensive, while computations involving A are relatively cheap. Therefore, many of the well-established preconditioners [START_REF] Benzi | Numerical solution of saddle point problems[END_REF][START_REF] Keller | Constraint preconditioning for indefinite linear systems[END_REF][START_REF] Murphy | A note on preconditioning for indefinite linear systems[END_REF] are not appropriate in our particular setting. Hence, we focus our attention on the inexact constraint preconditioner P [START_REF] Bergamaschi | Inexact constraint preconditioners for linear systems arising in interior point methods[END_REF][START_REF] Bergamaschi | Erratum to: Inexact constraint preconditioners for linear systems arising in interior point methods[END_REF], that is,

P = A BT B 0 = ⎛ ⎝ D 0 F-1 0 R H F-T HT 0 ⎞ ⎠ , ( 9 )
where

B ∈ R ×( +m) is a full row rank approximation of B, F ∈ R × , H ∈ R m× , approximations
of F and H, respectively. The simple choice

P 1 = ⎛ ⎝ D 0 F-1 0 R 0 F-T 0 0 ⎞ ⎠ with F-1 = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I n -I n I n -I n I n . . . . . . -I n I n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ (10)
has been found in [START_REF] Fisher | Parallelisation in the time dimension of four-dimensional variational data assimilation[END_REF] to provide good results in terms of rate of convergence of preconditioned Krylov subspace method. Therefore, this inexact constraint preconditioner can be used as a preconditioner that we will call a first-level preconditioner. Our next goal will be to improve it by updating techniques, that will be called second-level preconditioning techniques in what follows.

Second-level preconditioning for the saddle point approach

In this section, we are interested in further improving an initial constraint preconditioner P 1 so that the linear system at the (j + 1)-th nonlinear iteration reads

P -1 j+1 A j+1 u j+1 = P -1 j+1 f j+1 , with A j+1 = A B T j+1 B j+1 0 , P j+1 = A BT j+1 Bj+1 0 , (11)
with

δy j+1 = λ j+1 μ j+1 , u j+1 = δy j+1 δx j+1 , f j+1 = ⎛ ⎝ b j+1 d j+1 0 ⎞ ⎠ .
P j+1 is built from pairs of vectors collected during the iterative solution of the linear systems arising in the previous nonlinear iterations. This idea has been well explored on sequences of symmetric positive definite linear systems coming from GN iterations and has proved effective on large-scale applications [START_REF] Gratton | On a class of limited memory preconditioners for large scale linear systems with multiple right-hand sides[END_REF][START_REF] Morales | Automatic preconditioning by limited memory quasi-Newton updating[END_REF]. While constructing a new preconditioner for the (j + 1)th iteration, we must take into consideration the features of the saddle point system (8) as explained above, that is, computations involving the matrix B j+1 are expensive, while computations involving A are relatively cheap. Therefore, in Sections 3.2 and 3.3, we consider the class of inexact constraint preconditioners and focus on specific updates of the form

P j+1 = P j + 0 BT B 0 , with B ∈ R ×( +m) .

The direct and adjoint secant equations

Given a positive integer k and a set of vectors (δy i , δx i ) T (i = 1, . . . , k), respectively, we define

A B T j+1 B j+1 0 δy i δx i = g i h i , (i = 1, . . . , k). ( 12
)
We impose that P j+1 satisfies the k relations

A BT j+1 Bj+1 0
δy i δx i = g i h i , (i = 1, . . . , k). ( 13
)
Relations ( 12) and ( 13) simply lead to

BT j+1 δx i = g i -Aδy i , (i = 1, . . . , k), ( 14
)
Bj+1 δy i = h i , (i = 1, . . . , k). ( 15
)
We call Equations ( 14) and ( 15) the adjoint secant equation and the direct secant equation, respectively, in analogy with the corresponding notions in quasi-Newton methods. Since Bj+1 = Bj + B, we obtain from Equations ( 14) and ( 15) the set of 2k relations

BT δx i = g i -Aδy i -Bj T δx i := pi , (16)
Bδy i = h i -Bj δy i := qi . ( 17
)
To satisfy all these relations at once, we consider the following multiple secant equations

BT X = P, ( 18
)
BY = Q, ( 19
)
where

X ∈ R ×k , P ∈ R ( +m)×k , Y ∈ R ( +m)×k and Q ∈ R ×k defined as X = [δx 1 , . . . , δx k ], P = [p 1 , . . . , pk ], Y = [δy 1 , . . . , δy k ], Q = [q 1 , . . . , qk ], (20)
are supposed to be of full column rank. The consistency of these multiple secant equations depends on a condition relating X,Y,P and Q, as stated in the next lemma. [START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF] are of full column rank. A necessary and sufficient condition for the linear systems defined in Equations ( 18) and (19) to be consistent is

Lemma 3.1 Suppose that X ∈ R ×k , P ∈ R ( +m)×k , Y ∈ R ( +m)×k and Q ∈ R ×k defined in
P T Y = X T Q. ( 21
)
Proof For sake of readability, the proof is postponed to the appendix.

Finding the update B, that is, solving the secant equations given by ( 18) and ( 19) simultaneously, is not a very well covered area in the literature. Most of the studies focus on the case of square matrices and use only one secant equation at a time. To the best of our knowledge, two-sided equations are only considered in [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF][START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF], to obtain an approximation of the Jacobian matrix. In this study, we rather use the approximation for preconditioning, which we believe is new. Throughout the paper, we aim at finding certain solutions of the two-sided equations given by ( 18) and [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF]. We especially focus on the rank-k and rank-2k cases, giving a strong emphasis on the following three properties:

• Hereditary: When the new update is repeatedly used in a sequence of consistent secant equations, the new update still satisfies the previous secant equations.

• Transformation invariance:

The update is invariant with respect to linear transformations. For instance, in data assimilation, this can be a transformation of model variables to a set of control variables that are assumed to be uncorrelated [START_REF] Katz | Correlations of control variables in variational data assimilation[END_REF]. When such a transformation occurs, this property ensures that the update maintains all the relations of the transformed problem that are initially satisfied. • Least-change: The update has a least change characterization in terms of a particular matrix norm. In this study, we focus on the least change characterization in terms of a (possibly weighted) Frobenius norm.

A rank-k update

The proof in the appendix has derived a general form of the update B. Now, as in [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF][START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF], we specifically analyse the case of a rank-k update. Indeed, provided that P T Y is non-singular and P T Y = X T Q , a solution of both the adjoint and direct secant equations can be written as a rank-k update B of the form B = QZP T , [START_REF] Ingleby | Improved variational analyses using a nonlinear humidity control variable[END_REF] with Z ∈ R k×k being non-singular. This is stated in the next theorem.

Theorem 3.2 Given X ∈ R ×k , P ∈ R ( +m)×k , Y ∈ R ( +m)×k and Q ∈ R ×k of full column rank defined in [START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF], suppose that relationship [START_REF] Hólm | Assimilation and modelling of the hydrological cycle: ECMWF's status and plans[END_REF] holds and that P T Y is non-singular . Then the rank-k matrix of the form [START_REF] Ingleby | Improved variational analyses using a nonlinear humidity control variable[END_REF] given by

B = Q(P T Y ) -1 P T ( 23
) or B = Q(X T Q) -1 P T ( 24
)
is a solution of the multiple secant equations [START_REF] Gratton | On a class of limited memory preconditioners for large scale linear systems with multiple right-hand sides[END_REF] and [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF].

Proof The necessary and sufficient condition for the direct secant equation in Z

QZP T Y = Q (25)
to have a solution [26, Ch. 2, Theorem 13] reads

QQ † Q(P T Y ) † (P T Y ) = Q,
where † denotes the Moore-Penrose pseudo-inverse. This relation effectively holds, since Q is supposed to be of full column rank (the relation Q † Q = I k is then satisfied) and P T Y is nonsingular . Hence a possible solution of ( 25) can be obtained as [26, Ch. 2, Theorem 13]

Z = Q † Q(P T Y ) -1 , that is, Z = (P T Y ) -1
, which completes the proof. We note that the adjoint secant equation BT X = P holds due to Equation [START_REF] Hólm | Assimilation and modelling of the hydrological cycle: ECMWF's status and plans[END_REF].

Remark 1 When k = 1, we recover the rank-one update derived in [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF][START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF].

Remark 2 If the equality P T Y = X T Q does not hold, the update (23) obviously solves only the direct multiple secant equations [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF] exactly, whereas the update (24) solves only the adjoint multiple secant equations ( 18) exactly.

We have derived the rank-k block formula, when k pairs satisfy the multiple secant equations. Block formulas ensure that the secant equations do hold at once. Therefore, by construction, these updates satisfy the hereditary property. As an alternative to the block formula, one can recursively use the rank-1 update 1 as

Bi+1 = Bi + (h i -Bi δy i )(g i -Aδy i -BT i δx i ) T (h i -Bi δy i ) T δx i , B1 = B(j=0) (26)
to obtain Bj+1 = Bk+1 , where we have assumed that (h i -Bi δy i ) T δx i = 0 for i = 1, . . . , k. This rank-1 update is actually used in [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF] to update the Jacobian matrix in the framework of the full quasi-Newton method. It is referred as to the two-sided rank one update (TR1) update. A similar rank-1 update is also used in [START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF]. The properties of the TR1 update can be summarized as follows:

• The TR1 update generalizes the classical symmetric rank one (SR1) update [START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF][START_REF] Davidon | Variable metric method for minimization[END_REF], in the sense that this update reduces to the SR1 update, when Bj+1 is symmetric. • When the secant equations are consistent, it is shown in [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF] that the TR1 recursive formula satisfies the hereditary property, (i.e. it maintains the validity of all previous secant equations), and is invariant under linear transformations. • The TR1 update has no obvious least change characterization.

We note that, in general, the updates obtained from the recursive formula and the block formula are not necessarily equivalent [START_REF] Byrd | Representations of quasi-Newton matrices and their use in limited memory methods[END_REF]. However, concerning the TR1 update, when the secant equations are consistent, both block and recursive formulas become equivalent. In practice, since it is easier to monitor the values of the denominator (to avoid blow up of the update) from the recursive formula involving vectors rather than matrices as in the block formula, recursive formula are often preferred to the block ones.

A rank-2k update

Least-Frobenius norm update.

In this section, as an alternative to the TR1 update proposed in Section 3.2, we present a new update B, which satisfies a least change characterization in terms of the Frobenius norm. [START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF] , suppose that relationship [START_REF] Hólm | Assimilation and modelling of the hydrological cycle: ECMWF's status and plans[END_REF] holds. Then the solution of the minimization problem

Theorem 3.3 Given X ∈ R ×k , P ∈ R ( +m)×k , Y ∈ R ( +m)×k and Q ∈ R ×k of full column rank defined in
min B∈R ×( +m) B F subject to BT X = P, BY = Q, ( 27
)
is given by B * = X T † P T + (I -XX † )QY † , [START_REF] Murphy | A note on preconditioning for indefinite linear systems[END_REF] or equivalently

B * = (Y T † Q T + (I -YY † )PX † ) T . ( 29
)
Proof Any rank-2k solution B satisfying Equation ( 18) (the first constraint of the minimization problem [START_REF] Morales | Automatic preconditioning by limited memory quasi-Newton updating[END_REF]) can be written as [38, Lemma 2.1]

B = X T † P T + (I -XX † )C, ( 30
)
where C ∈ R ×( +m) . Inserting the relation [START_REF] Nowak | Efficient computation of linearized cross-covariance and auto-covariance matrices of interdependent quantities[END_REF] into the second constraint (Equation ( 19)) yields

X T † P T Y + (I -XX † )CY = Q.
The least Frobenius norm solution of the minimization problem min

C∈R ×( +m) (Q -X T † P T Y ) -(I -XX † )CY F is then given by [38, Lemma 2.3] C * = (I -XX † ) † (Q -X T † P T Y )Y † ,
or equivalently

C * = (I -XX † ) † QY † . ( 31
)
Since (I -XX † ) defines an orthogonal projection, substituting Equation (31) into Equation [START_REF] Nowak | Efficient computation of linearized cross-covariance and auto-covariance matrices of interdependent quantities[END_REF] yields the solution of the minimization problem as

B * = X T † P T + (I -XX † )QY † .
This rank-2k update can be also written as

B * = X T † P T + QY † -X (X T X ) -1 X T QY † .
Using Equation ( 21), we deduce

B * = X T † P T + QY † -X (X T X ) -1 P T YY † , = X T † P T + QY † -X T † P T YY † , = QY † + X T † P T (I -YY † ), = (Y T † Q T + (I -YY † )PX † ) T ,
which completes the proof.

Remark 3 In the case of a single secant equation, Equation [START_REF] Nowak | Efficient computation of linearized cross-covariance and auto-covariance matrices of interdependent quantities[END_REF] simplifies into the following rank-2 update

B * = δx T † p T + (I -δxδx † )C, ( 32
)
where C is a rank-2 matrix. Let us now investigate the role of C in the geometrical properties of the update B. To do so, we denote u = δx/ δx 2 and choose U ⊥ ∈ R ×( -1) such that U = [u, U ⊥ ] is an orthogonal matrix. Then Equation ( 32) can be written as

B * = up T δx 2 + U ⊥ U ⊥ T C.
Therefore, a matrix-vector multiplication with B * leads to a vector in the column space of U. When C is chosen as the zero matrix, the resulting vector has no components in the direction of the column vectors of U ⊥ , that is, BT * v = 0 for any vector v orthogonal to u and the rank-one Broyden update is recovered. When C = (I -δxδx † )αqp T , this yields the TR1 update. As shown in Theorem 3.3, choosing C = (I -δxδx † ) † qδy † results in a rank-two least-change update. Therefore, different choices of C define the information to be included from the directions orthogonal to the vector u.

We name the new update given in Theorem 3.3 as the least-Frobenius two-sided rank-2k update (FTR2). As an alternative to the block formulas, the FTR2 update can be recursively obtained from Equation (28) (with k = 1 and B1 = B(j=0) ) as

Bi+1 = Bi + δx i (g i -Aδy i -BT i δy i ) T δx T i δx i + (h i -Bi δy i )δy T i δy T i δy i - δx i δx T i (h i -Bi δy i )δy T i δx T i δx i δy T i δy i . ( 33
)
The properties of the FTR2 update can be summarized as follows:

• The FTR2 update generalizes the Powell-symmetric-Broyden (PSB) update [START_REF] Powell | A new algorithm for unconstrained optimization[END_REF], in the sense that the PSB update is recovered, when Bj+1 is symmetric. • As outlined before, the block FTR2 update satisfies the hereditary property. However, the recursive FTR2 (similar to the PSB update) does not satisfy the hereditary property. Hence, it satisfies only the last direct and adjoint secant equations. When using new pairs, the previous secant equations do not hold any longer. • The FTR2 update is affected by linear transformations on the variable domain.

• It is shown in Theorem 3.3 that the FTR2 update satisfies a least-change property with respect to the Frobenius norm.

The weighted least-Frobenius norm update.

In this section, we are interested in solving the minimization problem [START_REF] Morales | Automatic preconditioning by limited memory quasi-Newton updating[END_REF] in a weighted Frobenius norm, such that the update enjoys a transformation invariance property. The next theorem provides such an update. [START_REF] Haber | Quasi-Newton methods for large-scale electromagnetic inverse problems[END_REF] suppose that relationship (21) holds. Let S ∈ R ×k and T ∈ R ( +m)×k be such that X T S and Y T T are symmetric and positive definite matrices, respectively. Let W 1 ∈ R × and W 2 ∈ R ( +m)×( +m) be non-singular matrices such that the relations W 1 W T 1 X = S and W T 2 W 2 Y = T hold. Then the solution of the minimization problem

Theorem 3.4 Given X ∈ R ×k , P ∈ R ( +m)×k , Y ∈ R ( +m)×k and Q ∈ R ×k of full column rank defined in
min B∈R ×( +m) W -1 1 BW -1 2 F subject to BT X = P, BY = Q, ( 34
)
is given by B * = S(X T S)

-1 P T + (Q -S(X T S) -1 X T Q)(Y T T) -1 T T , ( 35
)
or equivalently

B * = (T(Y T T) -1 Q T + (P -T(Y T T) -1 Y T P)(X T S) -1 S T ) T . ( 36
) Proof We define ˆ B ∈ R ×( +m) as ˆ B = W -1 1 BW -1 2 .
Then the minimization problem (34) can be written as min

ˆ B∈R ×( +m) ˆ B F subjectto ˆ BT W 1 T X = W -T 2 P, ˆ BW 2 Y = W -1 1 Q.
This problem is structurally identical to the minimization problem [START_REF] Morales | Automatic preconditioning by limited memory quasi-Newton updating[END_REF]. Hence, using the update formula [START_REF] Murphy | A note on preconditioning for indefinite linear systems[END_REF] in this context yields

ˆ B * = (X T W 1 ) † (W -T 2 P) T + (I -W 1 T X (W 1 T X ) † )W -1 1 Q(W 2 Y ) † , = W 1 T X (X T W 1 W 1 T X ) -1 P T W -1 2 + (I -W 1 T X (X T W 1 W 1 T X ) -1 X T W 1 )W -1 1 Q(Y T W T 2 W 2 Y ) -1 Y T W T 2 . ( 37
)
Substituting the relations

W 1 T X = W -1 1 S, (38)
W 2 Y = W -T 2 T, ( 39
)
into Equation [START_REF] Shanno | Matrix conditioning and nonlinear optimization[END_REF], we obtain

ˆ B * = W -1 1 S(X T S) -1 P T W -1 2 + (W -1 1 Q -W -1 1 S(X T S) -1 X T Q)(Y T T) -1 T T W -1 2 , = W -1 1 [S(X T S) -1 P T + (Q -S(X T S) -1 X T Q)(Y T T) -1 T T ]W -1 2 .
From Equation ( 29), the solution of the minimization problem (34) can be written as

ˆ B * = ((Y T W T 2 ) † (W -1 1 Q) T + (I -W 2 Y (W 2 Y ) † )W -T 2 P(W 1 T X ) † ) T , = (W 2 Y (Y T W T 2 W 2 Y ) -1 Q T W -T 1 + (I -W 2 Y (Y T W 2 T W 2 Y ) -1 Y T W 2 T )W -T 2 P(X T W 1 W 1 T X ) -1 X T W 1 ) T . ( 40
)
Substituting Equations ( 38) and (39) into Equation ( 40), we obtain

ˆ B * = (W -T 2 T(Y T T) -1 Q T W -T 1 , + (W -T 2 P -W -T 2 T(Y T T) -1 Y T P)(X T S) -1 S T W -T 1 ) T , = (W -T 2 [T(Y T T) -1 Q T + (P -T(Y T T) -1 Y T P)(X T S) -1 S T ]W -T 1 ) T ,
which completes the proof. Due to Lemma 3.1, the formulas ( 35) and (36) are equivalent.

Remark 4

We note that X T S and Y T T may not be necessarily symmetric and positive definite for some given pairs S and T. However, the strategy proposed in [START_REF] Schnabel | Quasi-Newton methods using multiple secant equations[END_REF] can be applied in order to find least-change perturbations S and T such that X T (S + S) and Y T (T + T) become symmetric and positive definite. Theorem 3.4 can then be applied straightforwardly.

We name the update given by formulas [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] or [START_REF] Schnabel | Quasi-Newton methods using multiple secant equations[END_REF] as the weighted least Frobenius norm twosided rank-(2k) (WFTR2) update. We note that different choices of S and T (meaning different weighted Frobenius norms) lead to different updates.

As a special case (with k = 1 and B1 = B(j=0) ) of the block WFTR2 update, the recursive update can be obtained as

Bi+1 = Bi + s i (g i -Aδy i -BT i δx i ) T δx T i s i + (h i -Bi δy i )t T i δy T i t i - s i δx T i (h i -Bi δy i )t i T δx T i s i δy T i t i . ( 41
)
The properties of the WFTR2 update can be summarized as follows

• The block WFTR2 update satisfies the hereditary property, whereas the recursive formula does not. • Let us assume that we have applied the following change of variables

X = G -1 X and Ȳ = V -T Y
where G and V are non-singular matrices of order and ( + m), respectively. Provided that the relations S = G T S and T = VT hold, it can be simply shown by induction that the WFTR2 update is transformation-invariant.

• It is shown in Theorem 3.4 that the WFTR2 update satisfies a least-change property.

In this section, we are left with the choice of S and of T in the WFTR2 update. We note that when

S = [h 1 , . . . , h k ], (42)
T = [g 1 -Aδy 1 , . . . , g k -Aδy k ], (43)
with h i and g i (i = 1, . . . , k) defined in [START_REF] Fandry | A two-layer quasi-geostrophic model of summer through formation in the Australian subtropical easterlies[END_REF] and Bj+1 is a symmetric and positive definite matrix, the WFTR2 update reduces to the Davidon-Fletcher-Powell update [START_REF] Davidon | Variable metric method for minimization[END_REF][START_REF] Fletcher | A rapidly convergent descent method for minimization[END_REF]. Choosing S = Q and T = P leads to the TR1 update.

In our experiments, since the column vector matrices given by ( 42) are available as byproducts of the Krylov subspace method, we consider to choose these pairs for the WFTR2 update.

So far, we have introduced different update formulas detailed in Section 3.2 and 3.3 whose properties for the recursive formulas are summarized in Table 1. Concerning the block formulas, the hereditary property holds whatever the update, whereas the other two properties remain similar to those given in Table 1.

The inverse of the second-level preconditioner

In this section, we derive a practical formula concerning the inverse of the updated preconditioner to be used when solving the saddle point system [START_REF] Davidon | Variable metric method for minimization[END_REF]. The block update formulas of Theorems 3.2-3.4 can be written in the following generic form

B * = VZU T , ( 44
)
with V ∈ R ×p , Z ∈ R p×p and U ∈ R ( +m)×p , where p is equal to either k or 2k, respectively. For instance, the update (35) of Theorem 3.4 reads

B * = S Q -(X T S) -1 (X T Q)(Y T T) -1 (X T S) -1 (Y T T) -1 0 T T P T . ( 45
)
Using Equation ( 44), the inexact constraint preconditioner P 2 at the second nonlinear iteration can be thus updated as

P 2 = P 1 + 0 UZ T V T VZU T 0 , = P 1 + 0 U V 0 T ZU T 0 0 Z T V T G , ( 46
)
where T ∈ R (2 +m)×(2p) and G ∈ R (2p)× (2 +m) , respectively. Using the Sherman-Morrison-Woodbury formula, we obtain the inverse of the second-level preconditioner as

P -1 2 = P -1 1 -P -1 1 T (I + GP -1 1 T ) -1 GP -1 1 . ( 47
)
We note that this formula involves the application of the inverse of the first-level preconditioner P -1 1 , which is assumed to be easily tractable. This is indeed the case when considering the inexact constraint preconditioner as detailed next.

Construction and application of the second-level preconditioner

Let us assume that the following set of vectors (δx i , δy i , g i , h i ), (i = 1, . . . , k) satisfying Equation ( 12) have been saved during the previous outer loop of the approximate sequence of quadratic programming problems, and that a first-level preconditioner P 1 in the form of ( 9) is used. We present in Algorithm 1 both the setup phase and the application of the inverse of the second-level preconditioner on a given vector of appropriate dimension.

We want to discuss next a few questions related to the computational cost of the application of the first-and second-level preconditioners, respectively. The inverse of the first-level preconditioner [START_REF] Courtier | A strategy for operational implementation of 4D-Var using an incremental approach[END_REF] is given by

P -1 1 = ⎛ ⎝ 0 0 FT 0 R -1 0 F 0 -FD FT ⎞ ⎠ with F = ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ I n I n I n I n I n I n I n . . . . . . . . . I n . . . I n I n I n ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ . ( 48
)
Thus the action of P -1 1 requires one application of D, R -1 , F and FT in total. We note that applying F or FT is really cheap, since it only involves summations of vectors. Concerning the second-level preconditioner, we have decided to compute and store G, P -1 1 T Algorithm 1 Setup and application of the inverse of the second-level preconditioner # Setup phase

(1) Obtain P, Q, S and T from relations ( 16), ( 17), ( 42) and ( 43), respectively.

(2) Deduce T ∈ R (2 +m)×(2p) and G ∈ R (2p)×(2 +m) defined in (46) depending on the update formula.

(3) Compute P -1 1 T and (I 2p + GP -1 1 T ) -1 . # Application of the inverse of the second-level preconditioner on r ∈ R (2 +m) (1) Apply the first-level preconditioner

s = P -1 1 r. (2) Compute P -1 2 r = s -(P -1 1 T )(I 2p + GP -1 1 T ) -1 (G s).
and (I + GP -1 1 T ) -1 in memory, respectively. Consequently, the action of the inverse of the second-level preconditioner in Equation (47) involves only a single application of the first-level preconditioner P -1

1 and a matrix-vector product with G, (I + GP -1 1 T ) -1 and P -1 1 T , respectively (see Algorithm 1). By exploiting the structure of G and T , the total additional storage corresponds to 3p vectors of length (2 + m) and to 2p vectors of length 2p. This amount can be considered as affordable, since p is assumed to be much smaller than in this study. Computing P -1

1 T requires a single application of P -1 1 on a set of 2p vectors, while the computation of GP -1 1 T involves 8p 2 (2 + m) additional floating point operations. Both computations have to be performed at each nonlinear iteration. Applying the second-level preconditioner is thus slightly more expensive than the first-level preconditioner. Nevertheless, in our setting, we note that the application of A is the most expensive task in the Krylov subspace method since it involves multiple calculations of Jacobians. This will be illustrated in Section 4.

Spectrum of the preconditioned saddle point system

We briefly comment on the spectrum analysis of the saddle point system preconditioned by the inexact constraint preconditioner [START_REF] Courtier | A strategy for operational implementation of 4D-Var using an incremental approach[END_REF]. We have

P -1 1 A = ⎛ ⎝ FT F -T FT H T 0 0 I R -1 H FD(I -FT F -T ) -FD FT H T FF -1 ⎞ ⎠ .
A direct application of Corollary 2.2 of [START_REF] Bergamaschi | Erratum to: Inexact constraint preconditioners for linear systems arising in interior point methods[END_REF] reveals that the eigenvalues of P -1

1 A (denoted by λ(P -1 1 A)) are either equal to one or bounded by

|λ(P -1 1 A) -1| ≤ (F -T -F-T )D -1/2 H T R -1/2 2 σ 1 F-T D -1/2 0 ,
where σ 1 (A) denotes the smallest singular value of A. The last inequality is then equivalent to

|λ(P -1 1 A) -1| 2 ≤ λ max ((F -T -F-T )D -1 (F -1 -F-1 ) + H T R -1 H) λ min ( F-T D -1 F-1 ) . ( 49
)
When F = F, the eigenvalues of P -1 1 A are constrained to lie on the line x = 1 in the complex plane [START_REF] Fisher | Weak-constraint and long-window 4D-Var[END_REF]. When F is approximated, numerical computations reveal that the eigenvalues do lie on this line or belong to a disc of centre (1, 0) in the complex plane. The bound (49) gives us an insight on how the approximation F influences the radius of this disc.

We note that the inequality (49) is mainly of theoretical interest since the exact computation of the upper bound is out of reach for large-scale applications. In practice, we usually rely on the computation of Ritz or harmonic Ritz values to cheaply investigate the convergence of GMRES [START_REF] Van Der | Iterative Krylov Methods for Large Linear Systems[END_REF]. Numerical experiments in our setting have revealed that most of these values lie in the right-half plane, which is often found to be beneficial to the convergence of GMRES. Finally, we leave to a future study the spectral analysis of P -1 2 A, which we believe is outside the scope of the paper.

Numerical experiments

Our numerical experiments are performed by using a simple two-layer quasi-geostrophic model (QG-model) in the Object-Oriented Prediction System platform [START_REF] Rabier | Data assimilation in meteorology[END_REF]. In this section, we first describe the two layer QG-model, then explain the numerical set-up and finally present the numerical results.

A two-layer quasi-geostrophic model

This section describes the simple two-layer quasi-geostrophic 2 (nearly geostrophic) channel model which is widely used in theoretical atmospheric studies, since it is simple enough for numerical calculations and adequately captures the most relevant large-scale dynamics in the atmosphere.

The two-layer quasi-geostrophic model equations are based on the non-dimensional quasigeostrophic potential vorticity, whose evolution represents large-scale circulations of the atmosphere. The quasi-geostrophic potential vorticity on the first (upper) and second (lower) layers can be written, respectively, as

q 1 = ∇ 2 ψ 1 - f 2 0 L 2 g H 1 (ψ 1 -ψ 2 ) + βy, ( 50
)
q 2 = ∇ 2 ψ 2 - f 2 0 L 2 g H 2 (ψ 2 -ψ 1 ) + βy + R s , ( 51
)
where ψ is the stream function, ∇ 2 is the two-dimensional Laplacian operator, R s represents orography or heating, β is the (non-dimensionalized) northward variation of the Coriolis parameter at the fixed latitude y, f 0 is the Coriolis parameter at the southern boundary of the domain. L is the typical length scale of the motion we wish to describe, H 1 and H 2 are the depths of the two layers, g = g θ/ θ is the reduced gravity where θ is the mean potential temperature, and θ is the difference in potential temperature across the layer interface. Details of the derivation of these non-dimensional equations can be found in [START_REF] Fandry | A two-layer quasi-geostrophic model of summer through formation in the Australian subtropical easterlies[END_REF][START_REF] Pedlosky | Geophysical Fluid Dynamics[END_REF]. The conservation of potential vorticity in each layer is thus described by

D i q i Dt = 0, i = 1, 2. ( 52
)
This evolution equation is given in terms of the advective (total) derivative D i /Dt defined by

D i Dt = ∂ ∂t + u i ∂ ∂x + v i ∂ ∂y ,
where

u i = - ∂ψ i ∂y and v i = ∂ψ i ∂x , (53)
are the horizontal velocity components at each layer. Therefore, the potential vorticity at each time step is determined by using the conservation of potential vorticity given by Equation (52). In this process, time stepping consists of a simple first order semi-Lagrangian advection of potential vorticity. Given the potential vorticity at fixed time, Equations ( 50) and (51) can be solved for the stream function at each grid point and then the velocity fields obtained through Equation ( 53). Periodic boundary conditions in the west-east direction and Dirichlet boundary conditions in the north-south direction are employed. For the experiments in this paper, we choose L = 10 6 m, H 1 = 6000 m, H 2 = 4000 m, f 0 = 10 -4 s -1 and β = 1.5. For more details on the model and its solution, we refer the reader to [START_REF] Rabier | Data assimilation in meteorology[END_REF]. The domain for the experiments is 12,000 km by 6300 km for both layers. The horizontal discretization consists of 40 × 20 points, so that the east-west and the north-south resolution is approximately 300 km. The dimension of the state vector of the model is 1600.

Experimental set-up

Our numerical experiments are performed by using a two-layer model as described above. A reference stream function is generated from a model with layer depths of H 1 = 6000 m and H 2 = 4000 m, and the time step is set to 300 s, whereas the assimilating model has layer depths of H 1 = 5500 m and H 2 = 4500 m, and the time step is set to 3600 s. These differences in the layer depths and the time step provide a source of model error.

For all the experiments presented here, observations of the non-dimensional stream function, vector wind and wind speed were taken from the reference of the model at 100 points randomly distributed over both levels. Observations were taken every 12 hours. Observation errors were assumed to be independent from each others and uncorrelated in time, the standard deviations were chosen to be 0.4 for the stream function observation error, 0.6 for the vector wind and 1.2 for the wind speed. The observation operator is the bi-linear interpolation of the model fields to horizontal observation locations.

The background error covariance matrix (B b matrix) and the model error covariances (matrices Q i ) used in these experiments correspond to vertical and horizontal correlations. The vertical and horizontal structures are assumed to be separable. In the horizontal plane covariance matrices correspond to isotropic, homogeneous correlations of stream function with Gaussian spatial structure obtained from a Fast Fourier Transform approach [START_REF] Dietrich | Fast and exact simulation of stationary Gaussian processes through circulant embedding of the covariance matrix[END_REF][START_REF] Nowak | Efficient computation of linearized cross-covariance and auto-covariance matrices of interdependent quantities[END_REF]. For the background error covariance matrix B b , the standard deviation and the horizontal correlation length scale in this experiments are set to 0.8 and 10 6 m, respectively. For the model error covariance matrices Q i , the standard deviation and the horizontal correlation length scale are set to 0.2 and 2 × 10 6 m, respectively. The vertical correlation is assumed to be constant over the horizontal grid and the correlation coefficient value between the two layers was taken as 0.5.

The length of the assimilation window is set to 24 hours. The assimilation window is divided into three sub-windows (sub-intervals). Thus, the length of the each subwindow is equal to 8 hours and the model is considered as sufficiently accurate within each subwindow. The control variable is defined as the initial state of each subwindow. We have performed 3 GN iterations, and 10 inner loops at each GN iteration. Inner loops are solved by using full GMRES [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF].

Numerical results

We aim at solving the equality constrained least squares problem (4) by using our two-level preconditioning technique in which a sequence of saddle point systems (Equations ( 7)) is solved. The least squares problem solves for the stream function of the two layer quasi-geostraphic model as explained above by using a priori information x b and observations y i over the assimilation window (24 hours). Figure 1(a) shows the reference stream function (the truth) used to generate a priori information and observations, while Figure 1(b) shows a priori fields (forecast) for the stream function and potential vorticity that we have used in our experiments.

The computational times (in milliseconds) needed to apply the operators appearing in the saddle point system [START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF] are given in Table 2. These values have been obtained by running a numerical simulation on a single core of a desktop computer.

As outlined in Table 2, applying the Jacobian matrix of the QG-model M i (at time t i ) is more expensive by at least a factor of 70 than applying the covariance matrices D and R. Hence, this shows that the computational cost related to the application of either the first-or second-level preconditioner is significantly much lower than the cost of a matrix-vector product with the saddle point matrix of (7).

We have performed numerical experiments by using both the two levels of preconditioning, that is, the inexact preconditioner given by [START_REF] Courtier | A strategy for operational implementation of 4D-Var using an incremental approach[END_REF], and the second-level preconditioners obtained from the TR1 and WFTR2 updates, respectively. When constructing the updates, the last 8 pairs coming from the previous outer loop are used. Figure 2 shows the objective function values along the iterations when using these preconditioners. As seen especially from Figure 2(b), the effect of the second-level preconditioners can be remarked only at the third outer loop.

The improvement on large problems when using updates of Broyden type in nonlinear optimization is strongly affected by an improved numerical stability. It is indeed shown in [START_REF] Shanno | Matrix conditioning and nonlinear optimization[END_REF] that if the initial update matrix is unscaled, or scaled without knowledge of the actual magnitude of the elements of the true inverse Hessian, a severe loss of accuracy due to round-off errors can occur. Therefore, it is very common to use an initial scaling based on the inverse approximation of the Hessian to possibly improve the convergence [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Oren | Optimal conditioning of self-scaling variable metric algorithms[END_REF][START_REF] Shanno | Matrix conditioning and nonlinear optimization[END_REF][START_REF] Sun | Optimization Theory and Methods[END_REF]. One of the most successful choice for the initial scaling is shown to be an approximation of one of the eigenvalues of the inverse Hessian [START_REF] Nocedal | Numerical Optimization[END_REF][START_REF] Oren | Optimal conditioning of self-scaling variable metric algorithms[END_REF]. This choice actually attempts to make the range of the entries of the initial matrix similar to those of the inverse Hessian. With the same motivation in mind now for the initial preconditioner P 1 , we choose a scaling factor as

α = δy k δx k T g k h k g k h k T g k h k
, with pairs satisfying Equation ( 12) (k being the index of the last pair obtained from the previous minimization). By using Equation ( 12) and defining δw k = g k h k , α can be simply rewritten as a Rayleigh quotient

α = δw T k A -1 j+1 δw k δw k T δw k .
Therefore, this scaling factor attempts to make the magnitude of entries of P 1 closer to those of A -1 j+1 . Figure 3 shows the numerical results when the scaled preconditioner is used in our experiments. In this setting, the second-level preconditioners constructed from the scaled initial preconditioner accelerate the convergence also at the second outer loop. We stress that, for a large-scale data assimilation system, this gain is very crucial (here approximately 30%), since the cost of each iteration is very large (see Table 2). Thus, the role of the scaling factor for Table 3. RSE values for the forecast field and solutions obtained by using different preconditioners. the initial preconditioner may be quite significant in practical applications, since the operational 4D-Var is typically performed only with two or three outer loops. We want to mention that the performance of the second-level preconditioners is subject to small changes arising in the sequence of saddle point systems. Taking larger steps in an approximate sequence of quadrating problems changes the information crucially in the next saddle point system which prevents from using inherited information for the new system. Since in practical applications, only a limited number of inner iterations can be performed, we usually face a slowly varying sequence of saddle point systems which explains why these preconditioners are efficient in our experiments.

Examining the root square error (RSE) from the truth at the initial time given by

RSE = 1 n (truth -x * 0 ) T (truth -x * 0 ),
we notice that all preconditioners are able to reduce the error for a given number of iterations, the second-level preconditioners leading to a further reduction in the error (see Table 3). Finally, the difference fields are shown in Figure 4. We observe that both second-level preconditioners lead to very similar results for these experiments.

Conclusions

Geophysical inverse problems where a state vector is estimated using physical observations are very common in the simulation of complex systems such as those arising from meteorology or oceanography. In this study, we focus on the estimation of the state vectors at different times over an assimilation window in which the observations are assimilated. This process is known as a data assimilation problem. When perfect models are used for the dynamics, the estimation process leads to a large-scale unconstrained problem where the constraints are directly inserted into the objective function. The solution method is then based on a truncated GN technique. When accounting for model errors, this approach is known to be difficult to parallelize in the time dimension which makes it unsuitable for modern parallel computer architectures. It has been recently shown that a saddle point formulation of the estimation process can introduce highly parallel algorithms in the time dimension. In this paper, we present this formulation that involves a sequence of QP problems, in which the constraints represent both the system dynamics and the observation operator. Each QP results in a linear system where the matrix has the usual saddle point structure.

Efficient iterative solution algorithms for saddle point systems strongly rely on the availability of efficient preconditioners. The structure of the saddle point system arising in our application is markedly different from those usually considered in the literature, in the sense that the computational costs of matrix-vector products with the constraint blocks is much more important than those related to the diagonal blocks. For our preconditioning strategy we proposed in this paper to start with a low-cost first-level preconditioner for the first QP, and to further improve this preconditioner along the minimization process by using inherited information form earlier QP's. Due to the particular nature of our problem, we focus on the limited memory updates for such blocks. One interesting property of these blocks is that they are rectangular matrices, which is again a case that is less covered in the literature than the standard square case.

It is shown that these updates can be derived as a solution of the two-sided secant equations, namely adjoint and direct secant equations. We have also shown on numerical simulations that these second-level preconditioners do improve the convergence and that they are promising when solving either the sequence of saddle point systems with multiple right-hand sides or the sequence of slowly varying systems.

We are aware that open issues still need to be addressed. An important complementary aspect is related to the analysis of the properties of the resulting preconditioned operator to characterize the convergence of the Krylov subspace method. Another aspect is to exploit globalization strategies such as trust region methods or line searches strategies. Both issues are left as future lines of research.

Figure 1 .

 1 Figure 1. Stream function (top) and potential vorticity (bottom) values of (a) truth and (b) forecast ( a priori ) fields.

Figure 2 .

 2 Figure 2. (a) Objective cost function values along iterations for the first-and second-level preconditioners updated by using the TR1 or WFTR2 approaches. The first-level preconditioner is taken as the inexact constraint preconditioner P 1 . (b) Ratios of the cost function values associated with the second-level preconditioners to the cost function value associated with P 1 .

Figure 3 .

 3 Figure 3. (a) Objective cost function values along iterations for the first-and second-level preconditioners updated by using the TR1 or WFTR2 approaches. The first-level preconditioner is taken as the scaled P 1 . (b) Ratios of the cost function values associated with the second-level preconditioners shown in (a) to the cost function value associated with P 1 .

Figure 4 .

 4 Figure 4. Difference between the true solution and a priori (a), solution obtained with P 1 (b), solution obtained with TR1 update (c), solution obtained with WFTR2 update (d), respectively.

Table 1 .

 1 Properties of the different recursive updates of Sections 3.2 and 3.3.

	Update type	Rank	Hereditary	Trans. invariance	Least-change
	TR1 (Equation (26))	rank-1	yes	yes	no
	FTR2 (Equation (33))	rank-2	no	no	yes
	WFTR2 (Equation (41))	rank-2	no	yes	yes
	with pairs defined as (42)

Table 2 .

 2 Computational time (ms) associated with the application of the different operators appearing in Equation[START_REF] Broyden | Quasi-Newton methods and their application to function minimization[END_REF].

B b and Q i R M i and M T i H i and H T i Time (ms) 0.11 0.017 7.9 0.19

Acknowledgments

The authors wish to thank the two anonymous referees for their constructive comments which helped to improve the manuscript.

Notes

1. Since the derivations related to the update formula [START_REF] Katz | Correlations of control variables in variational data assimilation[END_REF] are quite similar, we omit the corresponding formula. 2. Quasi-geostrophic motion means that, in the horizontal direction of the atmospheric flow, the Coriolis force caused by the rotation of the Earth, and the pressure gradient force are in approximate balance.

Disclosure statement

No potential conflict of interest was reported by the authors.

Appendix . Proof of Lemma 3.1

Proof We note that Equations [START_REF] Gratton | On a class of limited memory preconditioners for large scale linear systems with multiple right-hand sides[END_REF] and [START_REF] Griewank | On constrained optimization by adjoint based quasi-Newton methods[END_REF] are consistent since both X and Y are supposed to be of full column rank [START_REF] Magnus | Matrix Differential Calculus with Applications in Statistics and Econometrics[END_REF]Chapter 2,Theorem 13]. We write the QR decompositions of the by k matrix X and of the ( + m) by k matrix Y as (1) , (A.1)

We note that both R (1) and R (2) are invertible k by k matrices, since X and Y are supposed to be of full column rank. Equations ( 18) and ( 19) can then be recast as

or equivalently as BT Q (1) = P(R (1) ) -1 , (A.3)

We refer to the byk matrix Q (1) ⊥ as the orthogonal complement of Q (1) in R , whereas the ( + m) by ( + mk) matrix Q (2) ⊥ denotes the orthogonal complement of Q (2) in R ( +m) , respectively. B can then be written as

where B is a by ( + m) matrix of the form

.

Postmultiplying Equation (A.5) by Q (2) leads to

Using Equation (A.4), we obtain

Similarly, by premultiplying Equation (A.5) by (Q (1) ) T , we obtain

Using Equation (A.3) then leads to

Finally, postmultiplying this last relation by Q (2) gives 11 = (R ( 1) ) -T P T Q (2) . (A.9)

From Equations (A.7) and (A.9), we conclude that

which completes the proof. Finally, we note that the proof is based on the equivalence of Equations (A.8) and (A.6) to Equations (A.3) and (A.4), respectively. Equations (A.8) and (A.6) helped us to determine specific blocks of rows or columns of B in the given decomposition basis. More precisely, we have determined a possible solution of the direct and adjoint secant equations as

⊥ 21 (Q (2) ) T + Q (1) ⊥ 22 (Q (2) ⊥ ) T , with 11 given in Equation (A.9), 12 = (R (1) ) -T P T Q (2) ⊥ , 21 = (Q (1) ⊥ ) T Q(R (2) ) -1 , and 22 a (k) by ( + mk) matrix that can be freely chosen.