Stanley C Eisenstat
email: stanley.eisenstat@yale.edu

Serge Gratton
email: serge.gratton@enseeiht.fr

David Titley-Peloquin
email: david.titley-peloquin@mcgill.ca

ON THE SYMMETRIC COMPONEN TWISE RELATIVE BACKWARD ERROR FOR LINEAR SYST EMS OF EQUATIONS *

Keywords:

We derive an upper bound on the symmetric componentwise relative backward error for symmetric linear systems of equations. Since the bound can be computed efficiently and, except for some artificial examples, seems to be of the same order of magnitude as the true symmetric componentwise backward error, we believe that it is suitable for practical use. Our results also provide new insight into the relationship between the symmetric and unsymmetric backward errors.

Introduction.

During an "open problems" session at the Householder Symposium XIX on Numerical Linear Algebra, held in Spa, Belgium, James Demmel presented the following challenge: Solve the symmetric componentwise relative backward error problem for linear systems of equations. In this note we present some progress that we have made toward its solution.

Suppose that we are given an approximate solution x ∈ R n to the linear system Au = b, where A ∈ R n×n is nonsingular, b ∈ R n , and u ∈ R n . The backward error problem is to find minimal (in some sense) perturbations ∆A ∈ R n×n and ∆b ∈ R n such that x is the exact solution of the perturbed problem (A + ∆A)x = b + ∆b.

The size of the smallest perturbation is called the backward error and can be measured normwise or componentwise. Here we consider a componentwise relative measure that is often used in conjunction with sparse direct solvers (see, e.g., [START_REF] Arioli | Solving sparse linear systems with sparse backward error[END_REF]).

When the matrix A has some structure, it is natural to require that the perturbation ∆A also have this structure. Here we consider the case where A is symmetric, and we require that ∆A be symmetric as well.

The unstructured componentwise and normwise backward error problems have been solved for nearly fifty years [START_REF] Oettli | Compatibility of approximate solution of linear equations with given error bounds for coefficients and right-hand sides[END_REF][START_REF] Rigal | On the compatibility of a given solution with the data of a linear system[END_REF]. The symmetric normwise backward error problem has also long been settled [START_REF] Bunch | The strong stability of algorithms for solving symmetric linear systems[END_REF]. The symmetric componentwise variant is more difficult: Although some results exist [START_REF] Smoktunowicz | A note on the strong componentwise stability of algorithms for solving symmetric linear systems[END_REF][START_REF] Higham | Backward error and condition of structured linear systems[END_REF][START_REF] Rump | The componentwise structure and unstructured backward errors can be arbitrarily far apart[END_REF], it remains open.

In this work we derive an upper bound on the symmetric componentwise relative backward error, whose precise definition is given in [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] below. The bound can be computed efficiently and, except for some artificial examples, seems to be of the same order of magnitude as the true symmetric backward error. Therefore, we believe that it is suitable for practical use. Our results also provide new insight into the relationship between the symmetric and unsymmetric backward errors.

Error bounds for the approximate solution x can be obtained by combining the structured backward error and a structured perturbation analysis. For an overview of structured perturbation results for linear systems of equations we refer to [START_REF] Higham | Backward error and condition of structured linear systems[END_REF][START_REF] Higham | A survey of componentwise perturbation theory in numerical linear algebra[END_REF][START_REF] Rump | Structured perturbations Part I: Normwise distances[END_REF][START_REF] Rump | Structured perturbations Part II: Componentwise distances[END_REF] or [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]Chapter 7] and the references therein.

1.1. Notation. We use subscripts to denote elements of vectors or subvectors of partitioned vectors, and similarly for matrices, while superscripts denote terms in a sequence. For example, A (k) ij denotes the (i, j)th element of the kth term in the sequence {A (k) }. For u ∈ R n and v ∈ R we define

diag (u) =    u 1 . . . u n    , sign (v) =      +1 if v > 0, -1 if v < 0, 0 if v = 0.
Absolute values and inequalities are meant componentwise; that is, |A| is the matrix whose entries are |A ij |, and |A| ≤ |B| means that |A ij | ≤ |B ij | for all pairs (i, j).

1.2. The componentwise backward error. Given A ∈ R n×n , b ∈ R n , and x ∈ R n , the componentwise relative backward error for the linear system Au = b is defined as

(1) ǫ * = min ǫ,∆A,∆b ǫ : (A + ∆A)x = b + ∆b, |∆A| ≤ ǫ|A|, |∆b| ≤ ǫ|b| .
Oettli and Prager [START_REF] Oettli | Compatibility of approximate solution of linear equations with given error bounds for coefficients and right-hand sides[END_REF] showed that

ǫ * = max i (|b -Ax|) i (|A| |x| + |b|) i
with the convention that 0/0 = 0. Since (|b -Ax|) i ≤ (|b| + |A||x|) i , it must be the case that 0 ≤ ǫ * ≤ 1 with ǫ * = 0 if and only if Ax = b.

Let r = b -Ax, let d be the vector whose elements are

(2)

d i = (|A||x| + |b|) i if (|A||x| + |b|) i = 0, 1 otherwise, and let (3)
D = diag (d) , z = D -1 r, Z = diag (z) , S = diag (sign (x)) . Then (4) ǫ * = z ∞ ,
and the minimum in (1) is attained by the perturbations 6) is identical to (1) except for the additional symmetry constraint, it must be the case that ǫ * sym ≥ ǫ * . Rump [START_REF] Rump | The componentwise structure and unstructured backward errors can be arbitrarily far apart[END_REF] has recently given an example (see ( 20)) where ǫ * sym can be arbitrarily larger than ǫ * . However, this result cannot be used to compute or bound ǫ * sym given a specific triplet (A, b, x). Variants of [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] have appeared in the literature. For example, Smoktunowicz [START_REF] Smoktunowicz | A note on the strong componentwise stability of algorithms for solving symmetric linear systems[END_REF] (see also [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]Problem 7.12]) considers perturbations only in the matrix A,

min ǫ,∆A ǫ : (A + ∆A)x = b, |∆A| ≤ ǫ|A| , (7a) min ǫ,∆A ǫ : (A + ∆A)x = b, |∆A| ≤ ǫ|A|, ∆A = (∆A) T , (7b)
and shows that if A is symmetric and diagonally dominant, the minimum in (7b) is at most a factor of three larger than that in (7a). This is also true for symmetric positive definite matrices [START_REF] Smoktunowicz | A note on the strong componentwise stability of algorithms for solving symmetric linear systems[END_REF]Theorem 3.1]. On the other hand, Higham and Higham [START_REF] Higham | Backward error and condition of structured linear systems[END_REF] give an example of a general symmetric A (see ( 19)) where the minimum in (7b) can be arbitrarily larger than that in (7a). Since perturbations to b are readily accounted for in structured perturbation analyses, we include them in [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF].

More generally, one can also replace |A| and |b| in (1) and ( 6) with arbitrary nonnegative tolerances E and f . For example, in the symmetric case Higham and Higham [START_REF] Higham | Backward error and condition of structured linear systems[END_REF] consider [START_REF] Oettli | Compatibility of approximate solution of linear equations with given error bounds for coefficients and right-hand sides[END_REF] min However, E = |A| and f = |b| is a natural and common choice, and for simplicity we have chosen to work with the less general (but still useful) formulation [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF]. Each of these backward error problems is a linear program (LP) and can in principle be solved using any general LP solver. Alternatively, Higham and Higham [START_REF] Higham | Backward error and condition of structured linear systems[END_REF] show how the solution can be computed as the minimum ∞-norm solution of an underdetermined linear system of equations. However, when the problem dimensions are large, these approaches are too expensive to be of practical use. It would be useful to have an expression for (or a tight upper bound on) the symmetric backward error ǫ * sym that can be computed cheaply, analogously to (4) in the unsymmetric case. This was the goal of the present work.

1.3. Our contributions. We make some progress toward solving (6) by deriving an upper bound ǫ ˜sym on ǫ s * ym . Combined with a struc tured perturbation analysis, it leads to upper bounds on the error in x. It can also b e used to compute a bound on the ratio ǫ s * ym /ǫ * for a given triplet (A, b, x).

In the following section we derive our upper bou nd. In section 3 we investigate whether it is tight. In section 4 we show how it can be computed efficiently, even when the problem dimensions are very large. In sect ion 5 we give some illustrative examples and a summary of our numerical experiment s.

2.

A bound on the symmetric componentwise backward error. We start by presenting our main result.

Theorem 2.1. Given a symmetric matrix A ∈ R n×n and vectors b ∈ R n and x ∈ R n , define r = b -Ax, the vector d as in (2), and the vector z and the matrices D and S as in [START_REF] Davis | The University of Florida sparse matrix collection[END_REF]. Then the matrix

(9) N = D -1 diag 1 2 |A||x| + |b| + 1 2 S|A|S diag (|x|)
is diagonally dominant with nonnegative diagonal elements, and the linear system N z = z is consistent. Moreover, for any solution z, the perturbations

(10) ∆A = 1 2 Z|A|S + S|A| Z , ∆b = -Z|b|, where Z = diag (z), satisfy (11)
(A + ∆A)x = b + ∆b, ∆A = ( ∆A) T ;
and the symmetric backward error ǫ * sym in (6) satisfies

(12) ǫ * sym ≤ max ǫA , ǫb ≤ z ∞ ≡ ǫsym , where (13) ǫA
= max i,j | ∆A ij | |A ij | , ǫb = max j | ∆b j | |b j |
with the convention that 0/0 = 0.

Note the similarity to the Oettli and Prager [START_REF] Oettli | Compatibility of approximate solution of linear equations with given error bounds for coefficients and right-hand sides[END_REF] result for the unsymmetric backward error. The final bound z ∞ in [START_REF] Rump | The componentwise structure and unstructured backward errors can be arbitrarily far apart[END_REF] has the same form as (4) and the perturbations ∆A and ∆b in [START_REF] Rump | Structured perturbations Part I: Normwise distances[END_REF] are the same as those in [START_REF] Higham | A survey of componentwise perturbation theory in numerical linear algebra[END_REF], except that z and Z have been replaced by z and Z, respectively, and ∆A has been symmetrized.

Proof. Consider the perturbations in [START_REF] Rump | Structured perturbations Part I: Normwise distances[END_REF]. Clearly ∆A is symmetric. In order to satisfy the other equality constraint in [START_REF] Rump | Structured perturbations Part II: Componentwise distances[END_REF] we must have

0 = (A + ∆A)x -(b + ∆b) = 1 2 Z|A|S + S|A| Z x + Z|b| -(b -Ax) = Z 1 2 |A|Sx + |b| + 1 2 S|A| Zx -r = diag 1 2 |A||x| + |b| z + 1 2 S|A|S diag (|x|) z -Dz = D N z -z .
Thus, if N z = z, then the perturbations ∆A and ∆b satisfy both equality constraints in [START_REF] Rump | Structured perturbations Part II: Componentwise distances[END_REF], and by [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] we have ǫ * sym ≤ max ǫA , ǫb .

From [START_REF] Rump | Structured perturbations Part I: Normwise distances[END_REF] it follows that ∆A ij = 0 if A ij = 0, and

| ∆A ij | |A ij | ≤ 1 2 |z i | + 1 2 |z j | ≤ z ∞
otherwise, while ∆b j = 0 if b j = 0, and

| ∆b j | |b j | = |z j | ≤ z ∞
otherwise. Therefore, ǫ * sym ≤ z ∞ = ǫsym . We now show that N is diagonally dominant with nonnegative diagonal entries. From (9) it follows that

N ii = 1 2 |A||x| + |b| i d i + |A ii ||x i | 2d i = |A ii ||x i | + |b i | d i + j =i |A ij ||x j | 2d i ≥ 0, while, for j = i, |N ij | = 0 if x i = 0, 1 2di |A ij ||x j | otherwise. Therefore (14) |N ii | - j =i |N ij | ≥ |A ii ||x i | + |b i | d i ≥ 0.
Finally we show by contradiction that the linear system N z = z is consistent. Suppose it is not. Then the matrix N = N (A, b, x) must be singular.

If any entries in x are zero, we can symmetrically permute the rows and columns so that

A = A 11 A T 21 A 21 A 22 , b = b 1 b 2 , x = x 1 0 ,
where |x 1 | > 0, and partition the other variables conformally. Then N z = z is equivalent to

N (A 11 , b 1 , x 1 ) 0 0 D -1 2 diag 1 2 |A 21 ||x 1 | + |b 2 | ) z1 z2 = z 1 z 2 .
Because

|z 2 | = D -1 2 r 2 = D -1 2 |b 2 -A 21 x 1 | ≤ 2D -1 2 1 2 |A 21 ||x 1 | + |b 2 | ,
if a diagonal entry in the (2, 2) block of the coefficient matrix above is zero, then the corresponding entry on the right-hand side must also be zero and that equation is consistent. If the diagonal entry is nonzero, that equation has a unique solution. That is, the second block of equations is consistent, so the system N (A 11 , b 1 , x 1 )z 1 = z 1 must be inconsistent. Thus without loss of generality we may assume that |x| > 0.

If A is reducible we can symmetrically permute the rows and columns so that

A = A 11 0 0 A 22
and partition the other variables conformally. Then N z = z is equivalent to

N (A 11 , b 1 , x 1 ) 0 0 N (A 22 , b 2 , x 2 ) z1 z2 = z 1 z 2 ,
and at least one of the systems N (A ii , b i , x i )z i = z i must be inconsistent. By applying the same reasoning recursively, we see that without loss of generality we may assume that A is irreducible. Since |x| > 0, the matrix N has the same off-diagonal nonzero structure as A and must also be irreducible; and by (14) it is diagonally dominant. If any entry in b were nonzero, then the corresponding row of N would be strictly diagonally dominant, so that N would be irreducibly diagonally dominant and nonsingular. Thus without loss of generality we may assume that b = 0.

Because N is singular, there exists a vector u = 0 such that

0 = N u = D -1 diag 1 2 |A||x| u + 1 2 D -1 S|A|S diag (|x|) u.
Let u i = 0 be an element of u of maximum magnitude. Then

0 = (N u) i = u i 2d i n j=1 |A ij ||x j | 1 + sign (x i ) sign (x j ) u j u i .
Since |u j | ≤ |u i |, all terms in the above sum are nonnegative. For the sum to be zero we must have A ii = 0 and u j = -u i sign (x i x j ) for every j for which |A ij | > 0; that is, u j must also be an element of maximum magnitude. Applying the same argument repeatedly, we must have |A kk | = 0 and u k = (-1) ℓ u i sign (x i x k ) whenever there is a path of length ℓ from node i to node k in the graph of A. Since A is irreducible, every node k is reachable and all entries in u are so defined. Thus u is unique up to a scaling factor, and the null-space of N is one-dimensional. Since the sign of u k depends only on whether the lengths of all paths from node i to node k in the graph of A are all even or all odd, the matrix A must have Young's property A: there exist two disjoint subsets S 1 and S 2 of W = {1, . . . , n} such that S 1 ∪ S 2 = W; and if A ij = 0, then i ∈ S 1 and j ∈ S 2 or i ∈ S 2 and j ∈ S 1 . It follows (see, e.g., [15, section 2.6]) that after symmetrically permuting the rows and columns of A we can write

A = 0 A T 21 A 21 0 .
Partitioning all other variables conformally, we have

N = 1 2 D -1 diag |A T 21 ||x 2 | S 1 |A T 21 |S 2 diag (|x 2 |) S 2 |A 21 |S 1 diag (|x 1 |) diag (|A 21 |x 1 |) . Let v = D [ x1 -x2 ]. Then N T v = 1 2 diag |A T 21 ||x 2 | x 1 -diag (|x 1 |) S 1 |A T 21 |S 2 x 2 diag (|x 2 |) S 2 |A 21 |S 1 x 1 -diag (|A 21 |x 1 |) x 2 = 0;
that is, v spans the null space of N T . On the other hand, since r = -Ax and z = D -1 r, we have

v T z = - x 1 -x 2 T A T 21 x 2 A 21 x 1 = x T 1 A T 21 x 2 -x T 2 A 21 x 1 = 0;
that is, z is orthogonal to the null space of N T . Thus z lies in the range of N and N z = z must have a solution, contrary to our original assumption.

3. Improvements to the bound. While Theorem 2.1 gives a simple expression for an upper bound, it is possible to improve on this result. We mention some of these improvements in this section.

First, we show how to eliminate zero elements on the diagonal of N , since such elements would complicate the algorithms for computing ǫsym that are presented in section 4. Suppose that N ii = 0 for some i. Since N is diagonally dominant, we must have N ij = 0 for all j, so that (N y) i = 0 for any y ∈ R n . Since the nonzero structure of N is symmetric, we must have N ji = 0 for all j; and since N z = z has a solution, we also have z i = (N z) i = 0. Thus the linear system N z = z does not place any constraints on zi . Changing N ii from 0 to 1 leaves N diagonally dominant with nonnegative diagonal entries but imposes the added constraint zi = 0. However, setting z i = 0 in any solution to the original system yields a solution to both the original and the new systems that does not increase z ∞ . After eliminating all zero diagonal entries in this manner, we get the following result.

Corollary 3.1. In the notation of Theorem 2.1, let the matrix N be defined by

Nij = 1 if i = j and N ii = 0, N ij otherwise.
Then N is diagonally dominant with positive diagonal entries, N z = z is consistent, and ǫ * sym ≤ z ∞ for any solution z. It is also possible to improve our bound by taking advantage of zeros in x, b, and the diagonal of A.

If any entries in x are zero, then without loss of generality we can write A, b, and x as

A = A 11 A T 21 A 21 A 22 , b = b 1 b 2 , x = x 1 0 ,
where |x 1 | > 0, and any symmetric perturbation ∆A and ∆b as

∆A = ∆A 11 ∆A T 21 ∆A 21 ∆A 22 , ∆b = ∆b 1 ∆b 2 .
To satisfy (A + ∆A)x = b + ∆b this perturbation must satisfy

(A 11 + ∆A 11 )x 1 = b 1 + ∆b 1 , (A 21 + ∆A 21 )x 1 = b 2 + ∆b 2 .
But replacing ∆A 21 and ∆b 2 by the Oettli-Prager perturbation for A 21 x 1 = b 2 and replacing ∆A 22 by the zero matrix gives another perturbation that satisfies these requirements without increasing ǫ. Indeed, the new minimum ǫ can be nearly a factor of 2 smaller (see ( 21)).

If corresponding entries on the diagonal of A and in b are 0, then the second term max{ǫ A , ǫb } in ( 12) can be as much as 2n -1 times smaller than ǫsym (see ( 22)).

If a diagonal entry of A is 0 and the corresponding entry of b is not, then ǫA can be smaller than ǫb . More generally, whenever ǫA < ǫb there is a gap between ǫ * sym and max{ǫ A , ǫb } that can be reduced by the following procedure. Suppose that we only allow perturbations in A as in problem (7b). Analogously to Theorem 2.1, let

Z = diag (ẑ) , ∆A = 1 2 Z|A|S + S|A| Z , ǫA = max i,j | ∆A ij | |A ij | ,
where ẑ is any vector such that

1 2 D -1 diag (|A||x|) + S|A|S diag (|x|) ẑ = z.
(Such a ẑ will exist if |b i | = 0 whenever (|A||x|) i = 0; the proof is similar to that of Theorem 2.1.) Then (A + ∆A)x = b and ∆A is symmetric, so that ǫA is also an upper bound on ǫ * sym . Now consider the blended perturbations

∆A (β) = β ∆A + (1 -β) ∆A, ∆b (β) = β ∆b,
for 0 ≤ β < 1. We still have

(A + ∆A (β) )x = b + ∆b (β) , ∆A (β) = (∆A (β) ) T , while max i,j |∆A (β) ij | |A ij | ≤ βǫ A + (1 -β)ǫ A , max j |∆b (β) j | |b j | ≤ βǫ b . Since ǫA < ǫb , ( 15
) ǫ * sym ≤ min 0≤β<1 max {βǫ A + (1 -β)ǫ A , βǫ b } = ǫA if ǫA < ǫA , β min ǫb otherwise,
where

β min = ǫA /(ǫ A + (ǫ b -ǫA )) < 1.
This bound is smaller than ǫb = max{ǫ A , ǫb } and ǫsym the bounds from Theorem 2.1 (see ( 23)).

In the tests performed in section 5.3, the bound ǫsym almost always gives a good order of magnitude estimate of ǫ * sym , so that these improvements are usually unneeded. However, the blending procedure may provide a much tighter bound in some particular examples (see (23)).

4. Efficient computation of the bound. The perturbations ∆A and ∆b in [START_REF] Rump | Structured perturbations Part I: Normwise distances[END_REF] and the upper bound ǫsym in [START_REF] Rump | The componentwise structure and unstructured backward errors can be arbitrarily far apart[END_REF] for the symmetric componentwise backward error are based on a solution z of the linear system N z = z with N given in [START_REF] Rigal | On the compatibility of a given solution with the data of a linear system[END_REF] and z in (3). In principle one can explicitly compute a z by solving the linear system directly, taking advantage of the fact that no pivoting is required for stability since N is diagonally dominant. However, this is not always practical. Furthermore, since N is based on |A| (as opposed to A), even if exact or inexact triangular factors of A are available, it is not clear how to reuse them to obtain a cheap estimate of z.

For large problems it is typically much more efficient to compute z using an iterative procedure. In practice one does not usually require the backward error and corresponding perturbations with a great deal of accuracy; often only an order-ofmagnitude estimate is sufficient. This makes the use of iterative strategies particularly attractive. Below we outline two such strategies that we have found to be quite efficient in our numerical experiments. To simplify the presentation we assume that N ii > 0 for all i. If this is not the case, one can simply replace N by N (see Corollary 3.1).

4.1. Gauss-Seidel. Write N = E + L + U , where E is diagonal with positive diagonal entries and L and U are strictly lower and upper triangular, respectively. To solve N z = z the Gauss-Seidel method with initial guess 0 computes the sequence of iterates

(16) z(k) = (E + L) -1 z -U z(k-1) , k > 0; z(0) = 0. Let ǫ(k) sym = z(k) ∞ . Since z(k) -z = -(E + L) -1 U z(k-1) -z , k > 0; z(0) -z = -z, we have ǫ(k) sym -ǫsym = z(k) ∞ -z ∞ ≤ z(k) -z ∞ ≤ (E + L) -1 U k ∞ ǫsym .
Since E is diagonal with positive diagonal entries,

(E + L) -1 U k ∞ ≤ (E -|L|) -1 |U | k ∞ = (E -|L|) -1 |U | k e ∞ ≡ α (k) ,
where e ∈ R n is the vector whose entries are all 1. Thus

(17) 1 + α (k) -1 ǫ(k) sym ≤ ǫsym ≤ 1 -α (k) -1 ǫ(k) sym .
The iteration can be stopped when the upper and lower bounds are within a small factor, say, 2, of each other, corresponding to α (k) ≤ 1/3. Note that α (k) can be computed as α (k) = q(k) ∞ , where

q(k) = (E -|L|) -1 |U |q (k-1) , k > 0; q(0) = e
are the Gauss-Seidel iterates for the linear system (E -|L| -|U |)q = 0 with initial guess e. This doubles the cost per iteration.

Convergence of Gauss-Seidel.

To show that Gauss-Seidel converges even when N is singular, we need a special case of Keller [7, Corollary 2.1].

Theorem 4.1. Let T be a symmetric nonnegative definite matrix with positive diagonal entries. Then the Gauss-Seidel iterates u (k) for the linear system T u = v converge for any initial guess u (0) . Moreover, if u satisfies T u = v and the initial error u (0) -u is orthogonal to the null space of T , then u (k) converges to u.

Since N need not be symmetric, we cannot apply this result directly. However, let P = diag (p), where

p i = |x i |d i if x i = 0, 1 otherwise.
Then PNP -1 is symmetric and has the same positiv e diagonal entries as N. Moreover, it has the same eigenvalues as N (by similarity), and these eigenvalues are real (by symmetry) and lie in the right half plane (by Gersh gorin's theorem applied to N).

That is, PNP -1 is symmetric nonnegative definite wi th positive diagonal entries. Let T = P 2 N. Since T = P (PNP -1 )P and P is diagonal, T is symmetric with positive diagonal entries, and by Sylvester's theorem i t is nonnegative definite. Since Gauss-Seidel is invariant under row scaling, the itera tes for T z ˜ = P 2 z with initial guess 0 are the same as those for Nz ˜ = z. Finally, s ince P is nonsingular, the null spaces of T and N are the same.

By Theorem 2.1 the linear system Nz ˜ = z is consi stent. Let z ¯ denote the projection of some solution z ˜ onto the orthogonal compleme nt of the null space of N. Then Nz ¯ = z so that T z ¯ = P 2 z. Since the initial error 0 -z ¯ is orthogonal to the null space of T , by Theorem 4.1 the Gauss-Seidel iterates must converge to z ¯.

A similar argument applied to (E -|L| -|U |)q = 0 shows that the Gauss-Seidel iterates q(k) with the initial guess e must converge. If E -|L| -|U | is nonsingular, its null space is {0}, and the initial error e -0 corresponding to the solution q = 0 is orthogonal to it. Thus the iterates converge to 0, as does α (k) . However, if E-|L|-|U | is singular, then it can be shown that neither sequence does, so that α (k) does not provide a useful bound.

GMRES.

While the Gauss-Seidel iterates always converge to a solution, they may converge slowly and α (k) need not converge to 0. An alternative is to use GMRES with initial guess 0 to solve the preconditioned system

G -1 L N G -1 R y = G -1 L z,
where G L and G R are nonsingular, and compute z =

G -1 R y. Let ỹ(k) denote the kth iterate, let z(k) = G -1 R ỹ(k) , and let ǫ(k) sym = z(k) ∞ .
If N is nonsingular, the error can be bounded as follows:

ǫ(k) sym -ǫsym ≤ z(k) -z ∞ ≤ G -1 R G R N -1 G L G -1 L N z(k) -z 2 ≤ G -1 R 2 G -1 L N G -1 R -1 2 G -1 L (N z(k) -z) 2 ≈ G -1 R 2 σ min (H k ) -1 G -1 L N z(k) -z 2 ,
where H k is the (k + 1) × k upper Hessenberg matrix generated by GMRES and the last factor is the norm of the residual of the preconditioned system. Letting

α (k) = G -1 R 2 σ min (H k ) -1 G -1 L N z(k) -z 2 z(k) ∞ , we get (18) 1 -α (k) ǫ(k) sym ǫsym 1 + α (k) ǫ(k) sym .
Again the upper bound and lower bound are within a factor of 2 of each other when

α (k) ≤ 1/3. Since σ min (H k ) -1 is a lower bound on (G -1 L N G -1 R ) -1
2 , this bound is only approximate, unlike that for the Gauss-Seidel iteration. However, the cost of computing it is negligible.

In our tests the best preconditioner was Gauss-Seidel (G L = E + L and G R = I, while taking advantage of the fact that

G -1 L N G -1 R = I + (E + L) -1 U
to reduce the cost of multiplying by the preconditione d matrix to that of one Gauss-Seidel iteration). Since Gauss-Seidel with initial guess z ˜(k) = 0 always converges and GMRES with the Gauss-Seidel preconditioner genera tes the same Krylov subspaces but chooses iterates that minimize the norm of the r esidual, the latter combination must also converge to a solution even when N is singu lar. As noted in section 4.2, the matrix PNP -1 is sy mmetric nonnegative definite. Thus GMRES with the preconditioner G L = P -1 and G R = P reduces to the conjugate residual method, which unlike GMRES requires a constant amount of work per iteration. However, in our experiments this advantage was not enough to overcome a slower rate of convergence.

Numerical experiments.

5.1. Some illustrative examples. We start by giving a few small examples that illustrate some important points.

Higham and Higham [START_REF] Higham | Backward error and condition of structured linear systems[END_REF] show that allowing perturbations in the right-hand side vector b can reduce the symmetric componentwise backward error. Let

(19) A = 1 1 1 0 , b = 1 δ , x = δ 1 ,
where δ > 0. If only perturbations in A are allowed, then that minimum is δ/(1 + δ) in the unsymmetric case (7a) and 1 in the symmetric case (7b). On the other hand, if perturbations in b are allowed, we have ǫ * = ǫ * sym = δ/(2 + δ). Furthermore, ǫsym = 3δ/(4 + 3δ), which is of the same order of magnitude.

In contrast, Rump [START_REF] Rump | The componentwise structure and unstructured backward errors can be arbitrarily far apart[END_REF] shows that the symmetric backward error can be arbitrarily larger than the unsymmetric one, even when perturbations to b are allowed. Let (20)

A =       δ 1 1 0 -1 1 0 1 -1 0 1 1 0 1 1 0 -1 1 0 1 -1 0 1 1 0       , x =       1 1 0 1 1       , b =       0 0 4 0 0       .
Then ǫ * = δ/(2 + δ), while ǫ * sym = 1. Furthermore, ǫsym = 1. As shown in section 3, we can take advantage of zeros in x. Let (21)

A =   2 -2 1 -2 2 -2 1 -2 2   , x =   1 1 0   , b =   1 1 δ   .
Then ǫsym = 2(1 + δ)/(3 + 2δ) corresponding to the perturbation

∆A =    2 5 2 5 1+δ 3+2δ 2 5 2 5 2(1+δ) 3+2δ 1+δ 3+2δ 2(1+δ) 3+2δ 0    , ∆b = -    1 5 1 5 δ(1+δ) 3+2δ    .
However, if we replace the third row/column of ∆A and the third row of ∆b by the Oettli-Prager perturbation for 1 -2 2 x = δ , then we get the perturbation

∆A =    2 5 2 5 1+δ 3+δ 2 5 2 5 2(1+δ) 3+δ 1+δ 3+δ 2(1+δ) 3+δ 0    , ∆b = -    1 5 1 5 δ(1+δ) 3+δ    , which corresponds to ǫ * sym = (1 + δ)/(3 + δ).
The latter is smaller than ǫsym by a factor of 2 asymptotically as δ → 0.

While z ∞ is easier to compute, max{ǫ A , ǫb } can be as much as 2n -1 times smaller. Let

(22) A =       1 -1 -1 0 1 1 . . . . . . . . . 0 (-1) n-1 (-1) n-1 0       , x =      1 1 . . . 1      , b =      0 0 . . . 0      .
Then ǫ * = ǫ * sym = max{ǫ A , ǫb } = 1, but ǫsym = 2n -1. As shown in section 3, blending can lead to tighter bounds. Let

(23) A =       1 -1 -1 0 1 1 . . . . . . . . . 0 (-1) n-1 (-1) n-1 (-1) n       , x =      1 1 . . . 1      , b = δ      -1 1 . . . (-1) n      .
Then ǫ * = δ/(2 + δ), while ǫ * sym = nδ/(2 + nδ). Thus the ratio ǫ * sym /ǫ * tends to n as δ → 0. Furthermore, for sufficiently small δ we have ǫsym ≈ n 2 δ/4 if n is even and ǫsym ≈ (n 2 + 1)δ/4 if n is odd. Thus the ratio ǫsym /ǫ * sym tends to roughly n/2 as δ → 0. On the other hand, by using the blending procedure and taking β = 0 in [START_REF] Young | Iterative Solution of Large Linear Systems[END_REF], we have ǫA = nδ/2. In other words, in this example,

ǫ * ≪ ǫ * sym ≈ ǫA ≪ ǫsym . 5
|A ij | |A ii | ,
and has also suggested that [14] (25)

ǫ * sym ≤ 2βǫ * , where β = max i (|A||x| + |b|) i |b i |
with the convention that α/0 = ∞ if α = 0 and 0/0 = 0. The bound (24) is ∞ when any diagonal entry of A is 0 (unless the entire row/column and right-hand side are also 0). As the tests in the next subsection show, it is not very tight.

The bound (25) is at least as large as z ∞ for any solution z of the linear system N z = z defined in Corollary 3.1. If (|A||x| + |b|) i = 0, then d i = 1 and

N ii = 1 2 |A||x| + |b| i + 1 2 |A ii ||x i | = 0 so that Nii = 1. Since N is diagonally dominant, Nij = N ij = 0 for all j = i. Thus we have (|I -N |e) i = 0. If (|A||x| + |b|) i > 0, then d i = (|A||x| + |b|) i , (I -N ) ii = 1 d i d i -1 2 |A||x| + |b| i -1 2 |A ii ||x i | = 1 2d i ((|A||x|) i -|A ii ||x i |) ,
and

j =i |(I -N ) ij | = j =i |N ij | ≤ 1 2d i j =i |A ij ||x j | = 1 2d i ((|A||x|) i -|A ii ||x i |) .
Thus we have

(|I -N |e) i ≤ 1 d i ((|A||x|) i -|A ii ||x i |) = 1 - |A ii ||x i | + |b i | (|A||x| + |b|) i .
Putting these together,

I -N ∞ = max i (|I -N |e) i ≤ max i : (|A||x|+|b|) i >0 1 - |A ii ||x i | + |b i | (|A||x| + |b|) i
with the convention that the maximum is 0 if

|A||x| + |b| = 0. If I -N ∞ < 1, then N -1 exists, N -1 ∞ ≤ 1 1 -I -N ∞ ≤ max i : (|A||x|+|b|) i >0 (|A||x| + |b|) i |A ii ||x i | + |b i | ≤ β,
and we have

ǫ * sym ≤ ǫsym = z ∞ ≤ N -1 ∞ z ∞ ≤ βǫ * . If I -N ∞ = 1, then |A ii ||x i | + |b i | = 0
for some i with (|A||x| + |b|) i > 0, so that the last expression is ∞ and the result still holds.

Numerical experiments.

It is difficult to compare ǫ * sym and ǫsym on large problems directly because the symmetric backward error ǫ * sym is too expensive to compute via linear programming. However, recall that

ǫ * ≤ ǫ * sym ≤ ǫsym ,
where ǫ * is the unsymmetric backward error in [START_REF] Arioli | Solving sparse linear systems with sparse backward error[END_REF]. Both ǫ * and ǫsym can be computed cheaply. Therefore, we can investigate the ratio ǫsym /ǫ * , which is an upper bound on both ǫsym /ǫ * sym and ǫ * sym /ǫ * , to gain insight into the relationship between ǫsym and ǫ * sym , as well as the relationship between ǫ * sym and ǫ * . We ran tests on the 589 matrices from the University of Florida Sparse Matrix Collection [START_REF] Davis | The University of Florida sparse matrix collection[END_REF] that satisfied the following criteria:

• A is real and symmetric.

• A is not binary (i.e., some A ij is neither 0 nor 1).

• A is not structurally singular (i.e., some permuted diagonal is nonzero).

• 100 ≤ n ≤ 100000 and 4 ≤ nnz(A)/n ≤ 200.

• At most 10 10 flops are needed to factor a symmetric positive definite matrix with the off-diagonal nonzero structure of A. After choosing x as described below, we created b as follows (in MATLAB notation): We also computed approximations to ǫsym by solving N z = z using Gauss-Seidel (GS) and GMRES with the Gauss-Seidel preconditioner (GM) (see section 4). Columns C GS and C GM give the first steps k at which ǫsym -ǫ(k) sym ≤ 1 10 ǫsym .

Columns D GS and D GM give the number of steps until convergence is detected by the practical stopping criterion α (k) ≤ 1/3 (see ( 17) and ( 18)). The reported means are rounded to the nearest integer.

The results are remarkable. The computed ratios ǫsym /ǫ * are close to 1, indicating that in these tests ǫ * sym ≈ ǫsym and ǫ * ≈ ǫ * sym . Convergence of Gauss-Seidel and Gauss-Seidel-preconditioned GMRES to an order-of-magnitude estimate of the symmetric backward error is achieved within a few iterations, and our practical stopping criteria are also triggered very quickly, especially for Gauss-Seidel-preconditioned GMRES.

In the above experiments, the choice x i = sin(i) ensures that |b| > 0. In an attempt to create problems in which ǫsym /ǫ * ≫ 1, we also ran tests with sparse vectors b. In MATLAB notation, w = sum(abs(A),2) ./ sum(A˜=0,2); b = w .* randn(n,1) .* (randperm(n)'<=q*n); for different values of the density q. Thus b contains qn randomly positioned and generated entries. For each such b we picked x to be the computed solution of Ax = b, computed using an LDL T factorization of A with scaling and pivoting. (In order to efficiently compute the factorization and reliably compute such x, of the previous 589 matrices we restricted ourselves to testing only those 343 numerically nonsingular matrices of dimension at most 20000.) Finally, we applied sparsity-and symmetrypreserving perturbations to A and b as in the previous set of tests, so that once again up to rounding errors (A + E)x = b + f . The results for this choice of x for p = 6 and density q = 0.1 are summarized in Table 2. The values listed in column (24)/ǫ * sym include only the 221 problems out of 343 in which the bound (24) is finite. In these tests, the ratio ǫsym /ǫ * is as large as 18.91, and more iterations are needed to achieve and detect convergence. (The number of iterations was capped at 98. A value of 99 indicates that the iterations failed to converge to the desired tolerance within 98 iterations.) While not as good as those in Table 1, these results nevertheless support the idea that ǫ * sym ≪ ǫsym and ǫ * ≪ ǫ * sym are very uncommon occurrences that arise from some special structure in A, b, and/or x.

Although we have not been able to fully quantify this observation, we offer the following partial explanation. Suppose that N is strictly diagonally dominant. Then

ǫ * sym = z ∞ = N -1 z ∞ ≤ N -1 ∞ z ∞ ≤ N -1 ∞ ǫ * .
For ǫ * sym to be much larger than ǫ * , N -1 ∞ must be large. But from ( 14)

N -1 ∞ ≤ 1 min i |N ii | -j =i |N ij | ≤ max i (|A||x| + |b|) i |A ii ||x i | + |b i | = 1 + max i j =i |A ij x j | |A ii ||x i | + |b i | .
Thus a necessary condition for ǫ * ≪ ǫ * sym is

j =i |A ij x j | ≫ |A ii ||x i | + |b i |
for some i. Although such examples can be created artificially, they seem to be extremely rare in practice.

6. Concluding remarks. We have given an upper bound on the symmetric componentwise relative backward error that can be computed efficiently and in our numerical experiments is usually of the same order of magnitude as the true symmetric componentwise backward error. Therefore, we believe that our bound is suitable for use in practice. It also provides new insight into the relationship between the symmetric backward error and the unsymmetric one, showing that both are usually of the same order of magnitude. Whether the techniques presented in this note can be extended to solve backward error problems involving structural properties other than symmetry remains to be seen.

( 5 ) 6 )

 56 ∆A * = Z|A|S, ∆b * = -Z|b|. See [6, section 7.2] for a short proof. Analogously to (1), for symmetric A the symmetric componentwise relative backward error for Ax = b is defined as ǫ * sym = min ǫ,∆A,∆b ǫ : (A + ∆A)x = b + ∆b, |∆A| ≤ ǫ|A|, |∆b| ≤ ǫ|b|, ∆A = (∆A) T . (Since the perturbations ∆A = -A and ∆b = -b satisfy the constraints with ǫ = 1, we have 0 ≤ ǫ * sym ≤ 1 with ǫ * sym = 0 if and only if Ax = b. Because (

  ǫ,∆A,∆b ǫ : (A + ∆A)x = b + ∆b, |∆A| ≤ ǫE, |∆b| ≤ ǫf, ∆A = (∆A) T .

  -p) * randn(n,1) .* b ; b = b-f;

Table 1

 1 Test results for p = 4 and x i = s in(i).Thus, up to rounding error, x satisfies (A + E)x = b + f , where E and f are random perturbations of A and b that preserve their sparsity structure and E is symmetric. Since the componentwise relative errors |E ij |/|A ij | and |f j |/|b j | are of the order of magnitude 10 -p , we expect ǫ * sym to be of the same order. The results for p = 4 and x i = sin(i) are summarized in Table1, but the results were similar for other values of p. Because b is generated randomly, all values reported are the mean, median, or maximum over the problems of the medians or maxima over 101 runs for each of the problems. The value of ǫsym was computed by solving N z = z directly. The values listed in column (24)/ǫ * sym include only the 415 problems out of 589 in which the bound (24) is finite.

	Matrix	ǫ * /10 -p ǫsym/ǫ * (24)/ǫ * sym C GS D GS C GM D GM
	Mean (median)	2.97	1.47	1 × 10 15	1	3	1	2
	Median (median)	3.08	1.43	3 × 10 4	1	3	1	2
	Maximum (median)	4.80	2.06	3 × 10 17	2	17	2	5
	Median (maximum)	4.21	2.02	4 × 10 4	2	3	2	3
	Maximum (maximum)	6.09	2.93	4 × 10 17	9	17	7	7

Table 2

 2 Test results for p = 6 and b spa rse.Matrixǫ * /10 -p ǫsym/ǫ * (24)/ǫ * sym C GS D GS C GM D GM

	Mean (median)	3.42	1.85	7 × 10 14	7	32	3	9
	Median (median)	3.42	1.75	2 × 10 4	1	4	2	3
	Maximum (median)	4.82	8.18	2 × 10 17	99	99	99	99
	Median (maximum)	4.75	2.42	2 × 10 4	2	4	2	3
	Maximum (maximum)	7.05	18.91	3 × 10 17	99	99	99	99

Acknowledgments. The authors would like to thank Xiao-Wen Chang, Alicja Smoktunowicz, as well as two anonymous referees for their very helpful comments.

The work of the authors was conducted with the support of the "Fondation Sciences et Technologies pour l'Aéronautique et l'Espace (STAE)," within the "Réseau Thématique de Recherche Avancée (RTRA)," Toulouse, France.


E Bergou

S Gratton
email: serge.gratton@enseeiht.fr

L N Vicente

Levenberg-Marquardt Methods Based o n Probabilistic Gradient Models and Inexact Subproblem Solution, with A pplication to Data Assimilation

Keywords: Levenberg-Marquardt method, nonlinear leas t squares, regularization, random models, inexactness, variational data assimilation, Kalman filter/smoother, ensemble Kalman filter/smoother

The Levenberg-Marquardt algorithm is one of t he most popular algorithms for the solution of nonlinear least squares problems. Motivated by the problem structure in data assimilation, we consider in this paper the extension of the classical Levenberg-Marquardt algorithm to the scenarios where the linearized least squares subproblems are solved inexactly and/or the gradient model is noisy and accurate only within a certain probability. Under appropriate assumptions, we show that the modified algorithm converges globally to a first order stationary point with probability one. Our proposed approach is first tested on simple problems where the exact gradient is perturbed with a Gaussian noise or only called with a certain probability. It is then applied to an instance in variational data assimilation where stochastic models of the gradient are computed by the so-called ensemble methods.

controls the size of the step, making Gauss-Newton globally convergent (i.e., convergent to stationarity independently of the starting point). We found that the regularization term added to Gauss-Newton maintains the structure of the linearized least squares subproblems arising in data assimilation, enabling us to use techniques like ensemble methods while simultaneously providing a globally convergent approach.

However, the use of ensemble methods in data assimilation poses difficulties since it makes random approximations to the gradient. We thus propose and analyze a variant of the Levenberg-Marquardt method to deal with probabilistic gradient models. It is assumed that an approximation to the gradient is provided but only accurate with a certain probability. The knowledge of the probability of the error between the exact gradient and the model one, and in particular of its density function, can be used in our favor in the update of the regularization parameter.

Having in mind large-scale applications (as those arising from data assimilation), we then consider that the least squares subproblems formulated in the Levenberg-Marquardt method are only solved in some approximated way. The amount of inexactness in such approximated solutions (tolerated for global convergence) is rigorously quantified as a function of the regularization parameter, in a way that it can be used in practical implementations.

We organize this paper as follows. In section 2, a short introduction to the Levenberg-Marquardt method is provided. The new Levenberg-Marquardt method based on probabilistic gradient models is described in section 3. Section 4 addresses the inexact solution of the linearized least squares subproblems arising within Levenberg-Marquardt methods. We cover essentially two possibilities: conjugate gradients (CGs) and any generic inexact solution of the corresponding normal equations. The whole approach is shown to be globally convergent to first order critical points in section 5, in the sense that a subsequence of the true objective function gradients goes to zero with probability one. The proposed approach is numerically illustrated in section 6 with a simple problem, artificially modified to create (i) a scenario where the model gradient is a Gaussian perturbation of the exact gradient, and (ii) a scenario case where to compute the model gradient both exact/approximated gradient routines are available but the exact one (seen as expensive) is called only with a certain probability.

An application to data assimilation is presented in section 7 where the purpose is to solve the 4DVAR problem using the methodology described in this paper. For the less familiar reader, we start by describing the 4DVAR incremental approach (Gauss-Newton) and the ways to solve the resulting linearized least squares subproblems, in particular Kalman smoother and EnKS, the latter one leading to stochastic model gradients. We then show how our approach, namely, the Levenberg-Marquardt method based on probabilistic gradient models and an inexact subproblem solution, provides an appropriate framework for the application of the 4DVAR incremental approach using the EnKS method for the subproblems and finite differences for derivative approximation. Illustrative numerical results using the Lorenz-63 model as a forecast model are provided.

A discussion of conclusions and future improvements is given in section 8. Throughout this paper • will denote the vector or matrix ℓ 2 -norm. The notation [X; Y ] will represent the concatenation of X and Y as in MATLAB syntax.

2. The Levenberg-Marquardt method. Let us consider the following general nonlinear least squares problem [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] min

x∈R n f (x) = 1 2 F (x) 2 ,
where F : R n → R m is a (deterministic) vector-valued function, assumed continuously differentiable, and m > n. Our main probabilistic approach to deal with nonlinear least squares problems is derived having in mind a class of inverse problems arising from data assimilation, for which the function f to be minimized in ( 1) is of the form

(2) 1 2 x 0 -x b 2 B -1 + T i=1 x i -M i (x i-1 ) 2 Q -1 i + T i=0 y i -H i (x i ) 2 R -1 i ,
where (x 0 , . . . , x T ) corresponds to x and where the operators M i and H i and the scaling matrices will be defined in section 7.

The Gauss-Newton method is an iterative procedure where at each point x j a step is computed as a solution of the linearized least squares subproblem min

s∈R n 1 2 F j + J j s 2 ,
where F j = F (x j ) and J j = J(x j ) denotes the Jacobian of F at x j . The subproblem has a unique solution if J j has full column rank, and in that case the step is a descent direction for f . In the case of our target application problem (2), such a linearized least squares subproblem becomes min δx 0 ,...,δx T ∈R n

1 2 x 0 + δx 0 -x b 2 B -1 + T i=1 x i + δx i -M i (x i-1 ) -M ′ i (x i-1 )δx i-1 2 Q -1 i + T i=0 y i -H i (x i ) -H ′ i (x i )δx i 2 R -1 i , (3)
where (δx 0 , . . . , δx T ) corresponds to s (and the other details are given in section 7).

The Levenberg-Marquardt method [START_REF] Levenberg | A method for the solution of certain problems in least squares[END_REF][START_REF] Marquardt | An algorithm for least-squares estimation of nonlinear parameters[END_REF] (see also [START_REF] Osborne | Nonlinear least squares-the Levenberg algorithm revisited[END_REF]) was developed to deal with the rank deficiency of J j and to provide a globalization strategy for Gauss-Newton. At each iteration it considers a step of the form -(J ⊤ j J j + γ j I) -1 J ⊤ j F j , corresponding to the unique solution of min

s∈R n m j (x j + s) = 1 2 F j + J j s 2 + 1 2 γ 2 j s 2 ,
where γ j is an appropriately chosen regularization parameter. See [18, Notes and References of Chapter 10] for a brief summary of theoretical and practical aspects regarding the Levenberg-Marquardt method. The Levenberg-Marquardt method can be seen as a precursor of the trust-region method [START_REF] Conn | Trust-Region Methods[END_REF] in the sense that it seeks to determine when the Gauss-Newton step is applicable (in which case the regularization parameter is set to zero) or when it should be replaced by a slower but safer gradient or steepest descent step (corresponding to a sufficiently large regularization parameter). The comparison with trust-region methods can also be drawn by looking at the square of the regularization parameter as the Lagrange multiplier of a trust-region subproblem of the form min s∈R n (1/2) F j + J j s 2 s.t. s ≤ δ j , and in fact it was soon suggested in [START_REF] Moré | The Levenberg-Marquardt algorithm: Implementations and theory[END_REF] to update the regularization parameter γ j in the same form as the trust-region radius δ j . For this purpose, one considers the ratio between the actual reduction f (x j ) -f (x j + s j ) attained in the objective function and the reduction m j (x j ) -m j (x j + s j ) predicted by the model, given by

ρ j = f (x j ) -f (x j + s j ) m j (x j ) -m j (x j + s j ) .
Then, if ρ j is sufficiently greater than zero, the step is accepted and γ j is possibly decreased (corresponding to "δ j is possibly increased"). Otherwise the step is rejected and γ j is increased (corresponding to "δ j is decreased").

3. The Levenberg-Marquardt method based on probabilistic gradient models. We are interested in the case where we do not have exact values for the Jacobian J j and the gradient J ⊤ j F j (of the model m j (x j + s) at s = 0), but rather approximations which we will denote by J m j and g m j . We are further interested in the case where these model approximations are built in some random fashion. We will then consider random models of the form M j , where g M j and J M j are random variables, and use the notation m j = M j (ω j ), g m j = g M j (ω j ), and J m j = J M j (ω j ) for their realizations. Note that the randomness of the models implies the randomness of the current point x j = X j (ω j ) and the current regularization parameter γ j = Γ j (ω j ), generated by the corresponding optimization algorithm.

Thus, the model (where F m j represents an approximation to F j )

m j (x j + s) -m j (x j ) = 1 2 F m j + J m j s 2 + 1 2 γ 2 j s 2 - 1 2 F m j 2 = g ⊤ m j s + 1 2 s ⊤ J ⊤ m j J m j + γ 2 j I s is a realization of M j (X j + s) -M j (X j ) = g ⊤ M j s + 1 2 s ⊤ J ⊤ M j J M j + Γ 2 j I s.
Note that we subtracted the order zero term to the model to avoid unnecessary terminology. Our subproblem then becomes

(4) min s∈R n m j (x j + s) -m j (x j ) = g ⊤ m j s + 1 2 s ⊤ J ⊤ m j J m j + γ 2 j I s.
In our data assimilation applied problem (2), the randomness arises from the use of EnKS to approximately solve the linearized least squares subproblem [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF]. In fact, as we will see in section 7, quadratic models of the form 1/2( u 2 (B N ) -1 + Hu -D 2 R -1 ) will be realizations of random quadratic models 1/2( u 2

B -1 + Hu -D 2 R -1 )
, where u corresponds to s, and where it would be easy to see what are the realizations g m and J m and the corresponding random variables g M and J M .

We will now impose that the gradient models g M j are accurate with a certain probability regardless of the history M 1 , . . . , M j-1 . The accuracy is defined in terms of a multiple of the inverse of the square of the regularization parameter (as happens in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF] for trust-region methods based on probabilistic models where it is defined in terms of a multiple of the trust-region radius). As we will see later in the convergence analysis (since the regularization parameter is bounded from below), one can demand less here and consider just the inverse of a positive power of the regularization parameter. Assumption 3.1. Given constants α ∈ (0, 2], κ eg > 0, and p ∈ (0, 1], the sequence of random gradient models {g M j } is (p)-probabilistically κ eg -first-order accurate, for corresponding sequences {X j }, {Γ j }, if the events

S j = g M j -J(X j ) ⊤ F (X j ) ≤ κ eg Γ α j
satisfy the following submartingale-like condition (5)

p * j = P (S j |F M j-1 ) ≥ p,
where F M j = σ(M 0 , . . . , M j-1 ) is the σ-algebra generated by M 0 , . . . , M j-1 . Correspondingly, a gradient model realization g m j is said to be κ eg -first-order accurate if

g m j -J(x j ) ⊤ F (x j ) ≤ κ eg γ α j .
The version of Levenberg-Marquart that we will analyze and implement takes a successful step if the ratio ρ j between actual and predicted reductions is sufficiently positive (condition ρ j ≥ η 1 below). In such cases, and now deviating from classical Levenberg-Marquart methods and following [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF], the regularization parameter γ j is increased if the size of the gradient model is of the order of the inverse of γ j squared (i.e., if g m j < η 2 /γ 2 j for some positive constant η 2 > 0). Another relevant distinction is that we necessarily decrease γ j in successful iterations when g m j ≥ η 2 /γ 2 j . The algorithm is described below and generates a sequence of realizations for the above-mentioned random variables. Algorithm 3.1 (Levenberg-Marquardt method based on probabilistic gradient models).

Initialization

Choose the constants η 1 ∈ (0, 1), η 2 , γ min > 0, λ > 1, and 0 < p min ≤ p max < 1. Select x 0 and γ 0 ≥ γ min . For j = 0, 1, 2, . . . 1. Solve (or approximately solve) (4), and let s j denote such a solution. 2. Compute ρ j = f (x j )-f (x j +s j ) m j (x j )-m j (x j +s j ) . 3. Make a guess p j of the probability p * j given in [START_REF] Conn | Trust-Region Methods[END_REF] such that p min ≤ p j ≤ p max .

If ρ j ≥ η 1 , then set x j+1 = x j + s j and

γ j+1 =      λγ j if g m j < η 2 /γ 2 j , max γ j λ 1-p j p j , γ min if g m j ≥ η 2 /γ 2 j .
Otherwise, set x j+1 = x j and γ j+1 = λγ j .

If exact gradients are used (in other words, if g M j = J(X j ) ⊤ F (X j )), then one always has

p * j = P 0 ≤ κ eg Γ α j F M j-1 = 1,
and the update of γ in successful iterations reduces to γ j+1 = max{γ j , γ min } (when g m j ≥ η 2 /γ 2 j ), as in the more classical deterministic-type Levenberg-Marquart methods. In general one should guess p j based on the knowledge of the random error incurred in the application context. It is however pertinent to stress that the algorithm runs for any guess of p j ∈ (0, 1] such that p j ∈ [p min , p max ].

4. Inexact solution of the linearized least squares subproblems. Step 1 of Algorithm 3.1 requires the approximate solution of subproblem [START_REF] Bocquet | An iterative ensemble Kalman smoother[END_REF]. As in trust-region methods, there are different techniques to approximate the solution of this subproblem yielding a globally convergent step, and we will discuss three of them in this section. For the purposes of global convergence it is sufficient to compute a step s j that provides a reduction in the model as good as the one produced by the so-called Cauchy step (defined as the minimizer of the model along the negative gradient or steepest descent direction -g m j ).

4.1.

A Cauchy step. The Cauchy step is defined by minimizing m j (x j -tg m j ) when t > 0 and is given by ( 6)

s c j = - g m j 2 g ⊤ m j (J ⊤ m j J m j + γ 2 j I)g m j g m j .
The corresponding Cauchy decrease (on the model) is

m j (x j ) -m j (x j + s c j ) = 1 2 g m j 4 g ⊤ m j (J ⊤ m j J m j + γ 2 j I)g m j . Since g ⊤ m j (J ⊤ m j J m j + γ 2 j I)g m j ≤ g m j 2 ( J m j 2 + γ 2 j ), we conclude that m j (x j ) -m j (x j + s c j ) ≥ 1 2 g m j 2 J m j 2 + γ 2 j .
The Cauchy step ( 6) is cheap to calculate as it does not require solving any system of linear equations. Moreover, the Levenberg-Marquart method will be globally convergent if it uses a step that attains a reduction in the model as good as a multiple of the Cauchy decrease. Thus we will impose the following assumption on the step calculation. Assumption 4.1. For every step j and for all realizations m j of M j ,

m j (x j ) -m j (x j + s j ) ≥ θ f cd 2 g m j 2 J m j 2 + γ 2 j
for some constant θ f cd > 0.

Such an assumption asks from the step a very mild reduction on the model (a fraction of what a step along the negative gradient would achieve) and it can thus be seen as a sort of minimum first order requirement.

4.2.

A truncated-CG step. Despite providing a sufficient reduction in the model and being cheap to compute, the Cauchy step is a particular form of steepest descent, which can perform poorly regardless of the step length. One can see that the Cauchy step depends on J ⊤ m j J m j only in the step length. Faster convergence can be expected if the matrix J ⊤ m j J m j also influences the step direction.

Since the Cauchy step is the first step of the CG method when applied to the minimization of the quadratic m j (x j + s) -m j (x j ), it is natural to propose to run CG further and stop only when the residual becomes relatively small. Since CG generates iterates by minimizing the quadratic over nested Krylov subspaces, and the first subspace is the one generated by g m j (see, e.g., [START_REF] Nocedal | Numerical Optimization[END_REF]Theorem 5.2]), the decrease attained at the first CG iteration (i.e., by the Cauchy step) is kept by the remaining.

A step from inexact solution of normal equations.

Another possibility to approximately solve subproblem ( 4) is to apply some iterative solver (not necessarily CG) to the solution of the normal equations

J ⊤ m j J m j + γ 2 j I s j = -g m j .
An inexact solution s in j is then computed such that

(7) J ⊤ m j J m j + γ 2 j I s in j = -g m j + r j
for a relatively small residual r j satisfying r j ≤ ǫ j g m j . For such sufficiently small residuals we can guarantee a Cauchy decrease.

Assumption 4.2. For some constants β in ∈ (0, 1) and θ in > 0, suppose that r j ≤ ǫ j g m j and

ǫ j ≤ min θ in γ α j , β in γ 2 j J m j 2 + γ 2 j .
Note that we only need the second bound on ǫ j (see the above inequality) to prove the desired Cauchy decrease. The first bound on ǫ j will be used later in the convergence analysis. The following result is proved in the appendix. 5. Global convergence to first order critical points. We start by stating that two terms, that later will appear in the difference between the actual and predicted decreases, have the right order accuracy in terms of γ j . The proof is given in the appendix. Lemma 5.1. For the three steps proposed (Cauchy, truncated CG, and inexact normal equations), one has that

s j ≤ 2 g m j γ 2 j and |s ⊤ j (γ 2 j s j + g m j )| ≤ 4 J m j 2 g m j 2 + 2θ in g m j 2 min{1, γ 2-α min }γ 2+α j .
(Assumption 4.2 is assumed for the inexact normal equations step s j = s in j .) We proceed by describing the conditions required for global convergence.

Assumption 5.1. The function f is continuously differentiable in an open set containing L(x 0 ) = {x ∈ R n : f (x) ≤ f (x 0 )} with Lipschitz continuous gradient on L(x 0 ) and corre- sponding constant ν > 0.
The Jacobian model is uniformly bounded, i.e., there exists κ Jm > 0 such that J m j ≤ κ Jm for all j.

The next result is a classical one and essentially says that the actual and predicted reductions match each other well for a value of the regularization parameter γ j sufficiently large relatively to the size of the gradient model (which would correspond to a sufficiently small trust-region radius in trust-region methods). Lemma 5.2. Let Assumption 5.1 hold. Let also Assumption 4.2 hold for the inexact normal equations step s j = s in j . If x j is not a critical point of f and the gradient model g m j is κ egfirst-order accurate, and if

γ j ≥ κ j 1 -η 1 1 α with κ j = 1 + κ 2 Jm γ 2 min 2ν + 2κeg gm j + 2θ in + 8κ 2 Jm min{1, γ 2-α min }θ f cd , then ρ j ≥ η 1 .
Proof. Again we omit the indices j in the proof. Applying a Taylor expansion,

1 - ρ 2 = m(x) -f (x) + f (x + s) -m(x + s) + m(x) -m(x + s) 2[m(x) -m(x + s)] = s ⊤ J(x) ⊤ F (x) + R -s ⊤ g m -s ⊤ (J ⊤ m J m + γ 2 I)s -s ⊤ g m 2[m(x) -m(x + s)] = R + (J(x) ⊤ F (x) -g m ) ⊤ s -s ⊤ (J ⊤ m J m )s -s ⊤ (γ 2 s + g m ) 2[m(x) -m(x + s)]
,

where R ≤ ν s 2 /2. Now, using Lemma 5.1, Assumptions 4.1 and 5.1, and γ ≥ γ min ,

1 - ρ 2 ≤ ν 2 s 2 + κeg γ α s + J m 2 s 2 -s ⊤ (γ 2 s + g) θ f cd gm 2 Jm 2 +γ 2 ≤ 2ν gm 2 γ 4 + 2κeg gm γ 2+α + 4κ 2 J m gm 2 γ 4 + 4κ 2 J m gm 2 +2 gm 2 θ in min{1,γ 2-α min }γ 2+α θ f cd gm 2 γ 2 ( Jm 2 /γ 2 min +1) ≤ 1 + κ J m γ 2 min 2ν + 2κeg gm + 2θ in + 8κ 2 Jm min{1, γ 2-α min }θ f cd γ α ≤ κ γ α ≤ 1 -η 1 .
We have thus proved that ρ ≥ 2η 1 > η 1 .

One now establishes that the regularization parameter goes to infinity, which corresponds to the trust-region radius going to zero in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF].

Lemma 5.3. Let the second part of Assumption 5.1 hold (the uniform bound on J m j ). For every realization of the Algorithm 3.1, lim j→∞ γ j = ∞.

Proof. If the result is not true, then there exists a bound B > 0 such that the number of times that γ j < B happens is infinite. Because of the way γ j is updated one must have an infinity of iterations such that γ j+1 ≤ γ j , and for these iterations one has ρ j ≥ η 1 and g m j ≥ η 2 /B 2 . Thus,

f (x j ) -f (x j + s j ) ≥ η 1 [m j (x j ) -m j (x j + s j )] ≥ η 1 θ f cd 2 1 J m 2 + γ 2 g m j 2 ≥ η 1 θ f cd 2(κ 2 Jm + B 2 )
η 2 B 2 2 .
Since f is bounded from below by zero, the number of such iterations cannot be infinite, and hence we arrive at a contradiction. Now, if we assume that the gradient models are (p j )-probabilistically κ eg -first-order accurate, we can show our main global convergence result. First we will state an auxiliary result from the literature that will be useful for the analysis (see [START_REF] Durrett | Probability: Theory and Examples[END_REF]Theorem 5.3.1] and [8, Exercise 5.3.1]).

Lemma 5.4. Let G j be a submartingale, in other words, a set of random variables which are integrable (E(|G j |) < ∞) and satisfy E(G j |F j-1 ) ≥ G j-1 for every j, where F j-1 = σ(G 0 , . . . , G j-1 ) is the σ-algebra generated by G 0 , . . . , G j-1 and E(G j |F j-1 ) denotes the conditional expectation of G j given the past history of events F j-1 .

Assume further that there exists M > 0 such that |G j -G j-1 | ≤ M < ∞ for every j. Consider the random events C = {lim j→∞ G j exists and is finite} and

D = {lim sup j→∞ G j = ∞}. Then P (C ∪ D) = 1.
Theorem 5.1. Let Assumption 5.1 hold. Let also Assumption 4.2 hold for the inexact normal equations step s j = s in j .

Suppose that the gradient model sequence {g M j } is (p j )-probabilistically κ eg -first-order accurate for some positive constant κ eg (Assumption 3.1). Let {X j } be a sequence of random iterates generated by Algorithm 3.1. Then almost surely, lim inf j→∞ ∇f (X j ) = 0.

Proof. The proof follows the same lines as [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF]Theorem 4.2]. Let

W j = j i=0 1 p i 1 S i -1 ,
where S i is as in Assumption 3.1. Recalling p * j = P (S j |F M j-1 ) ≥ p j , we start by showing that {W j } is a submartingale:

E(W j |F M j-1 ) = W j-1 + 1 p j P (S j |F M j-1 ) -1 ≥ W j-1 . Moreover, min{1, 1/p j -1} ≤ |W j -W j-1 | ≤ max{(1 -p j )/p j , 1} ≤ max{1/p j , 1} = 1/p j . Since 0 < p min ≤ p j ≤ p max < 1, one has 0 < min{1, 1/p max -1} ≤ |W j -W j-1 | ≤ 1/p min .
Thus, from 0 < min{1, 1/p max -1} ≤ |W j -W j-1 |, the event {lim j→∞ W j exists and is finite} has probability zero, and using Lemma 5.4 and

|W j -W j-1 | ≤ 1/p min , one concludes that P (lim sup j→∞ W j = ∞) = 1.
Suppose there exist ǫ > 0 and j 1 such that, with positive probability, ∇f (X j ) ≥ ǫ for all j ≥ j 1 . Let now {x j } and {γ j } be any realization of {X j } and {Γ j }, respectively, built by Algorithm 3.1. By Lemma 5.3, there exists j 2 such that ∀j ≥ j 2 (8)

γ j > b ǫ = max 2κ eg ǫ 1 α , 2η 2 ǫ 1 2 , λ p-1 p γ min , κ ǫ 1 -η 1 1 α
, where

κ ǫ = 1 + κ 2 Jm γ 2 min 2ν + 4κeg ǫ + 2θ in + 8κ 2 Jm min{1, γ 2-α min }θ f cd .
For any j ≥ j 0 = max{j 1 , j 2 } two cases are possible. If 1 S j = 1, then, from [START_REF] Durrett | Probability: Theory and Examples[END_REF],

g m j -J(x j ) ⊤ F (x j ) ≤ κ eg γ α j < ǫ 2 ,
yielding g m j ≥ ǫ/2. From (8) we also have that g m j ≥ ǫ/2 ≥ η 2 /γ 2 j . On the other hand, Lemma 5.2, (8), and g m j ≥ ǫ/2 together imply that ρ j ≥ η 1 . Hence, from this and step 3 of Algorithm 3.1, the iteration is successful. Also, from g m j ≥ η 2 /γ 2 j and (8) (note that (1 -x)/x is decreasing in (0, 1]), γ is updated in step 3 as

γ j+1 = γ j λ 1-p j p j
.

Let now B j be a random variable with realization b j = log λ (b ǫ /γ j ). In the case 1

S j = 1, b j+1 = b j + 1 -p j p j .
If 1 S j = 0, then b j+1 ≥ b j -1, because either γ j+1 ≤ γ j and therefore b j+1 ≥ b j , or γ j+1 = λγ j and therefore b j+1 ≥ b j -1. Hence B j -B j 0 ≥ W j -W j 0 , and from P (lim sup j→∞ W j = ∞) = 1 one obtains P (lim sup j→∞ B j = ∞) = 1 which leads to a contradiction with the fact that B j < 0 happens for all j ≥ j 0 with positive probability.

6. A numerical illustration. The main concern in the application of Algorithm 3.1 is to ensure that the gradient model is (p j )-probabilistically accurate (i.e., p * j ≥ p j ; see Assumption 3.1) or at least to find a lower bound p min > 0 such that p * j ≥ p min . However, one can, in some situations, overcome these difficulties such as in the cases where the model gradient (i) is a Gaussian perturbation of the exact one, or (ii) results from using either the exact one (seen as expensive) or an approximation. In the former case we will consider a run of the algorithm under a stopping criterion of the form γ j > γ max .

6.1. Gaussian noise. At each iteration of the algorithm, we consider an artificial random gradient model, by adding to the exact gradient an independent Gaussian noise, more precisely we have g M j = J(X j ) ⊤ ∇F (X j ) +ε j , where (ε j ) i ∼ N (0, σ 2 j,i ) for i = 1, . . . , n. Let Σ j be a diagonal matrix with diagonal elements σ j,i , i = 1, . . . , n. It is known that

Σ j ε j 2 = n i=1 (ε j ) i σ j,i 2 ∼ χ 2 (n),
where χ 2 (n) is the chi-squared distribution with n degrees of freedom. To be able to give an explicit form of the probability of the model being κ eg -first-order accurate, for a chosen κ eg > 0, we assume also that the components of the noise are identically distributed, that is, σ j,i = σ j ∀i ∈ {1, . . . , n}. Because of the way in which γ j is updated in Algorithm 3.1, it is bounded by λ j γ 0 and, thus, Γ j ≤ min{λ j γ 0 , γ max }, where γ max is the constant used in the stopping criterion. One therefore has

p * j = P g M j -J(X j ) ⊤ F (X j ) ≤ κ eg Γ α j F M j-1 ≥ P Σ j ε j 2 ≤ κ eg σ j min{λ j γ 0 , γ max } α 2 F M j-1 .
Using the Gaussian nature of the noise ε j and the fact that it is independent from the filtration F M j-1 , we obtain

p * j ≥ CDF -1 χ 2 (n) κ eg σ j min{λ j γ 0 , γ max } α 2 def = pj , (9)
where CDF χ 2 (n) is the cumulative density function of a chi-squared distribution with n degrees of freedom. The numerical illustration was done with the following nonlinear least squares problem defined using the well-known Rosenbrock function

f (x, y) = 1 2 x -1 2 + 100 y -x 2 2 = 1 2 F (x, y) 2 .
The minimizer of this problem is (x * , y * ) ⊤ = (1, 1) ⊤ . Algorithm 3.1 was initialized with x 0 = (1.2, 0) ⊤ and γ 0 = 1. The algorithmic parameters were set to η 1 = η 2 = 10 -3 , γ min = 10 -6 , and λ = 2. The stopping criterion used is γ j > γ max , where γ max = 10 6 . We used α = 1/2, σ j = σ = 10 ∀j, and κ eg = 100 for the random gradient model.

Figure 1 depicts the average, over 60 runs of Algorithm 3.1, of the objective function values, the absolute errors of the iterates, and the percentages of successful iterations, using, across all iterations, the three choices p j = 1, p j = pj , and p j = p min . In the last case, p min is an underestimation of p * j given by

p min = CDF -1 χ 2 (n) κ eg σγ α max 2 = 5 • 10 -3 .

Table 1

For three different runs of Algorithm 3.1, the table shows the values of the objective function and relative error of the solution found for the three choices pj = 1, pj = pj, and pj = pmin = 5 • 10 -3 . The final objective function values and the relative final errors are shown in Table 1 for the first three runs of the algorithm. One can see that the use of p j = pj leads to a better performance than p j = p min (because pj ≥ p min is a better bound for p * j than p min is). In the case where p j = 1, Algorithm 3.1 exhibits a performance worse than for the two other choices of p j . The algorithm stagnated after some iterations, and could not approximate the minimizer with a decent accuracy. In this case, γ j is increasing along the iterations, and thus it becomes very large after some iterations while the step s j ∼ 1/γ 2 j becomes very small. Other numerical experiments (not reported here) have shown that, when the error on the gradient is small (σ ≪ 1), the two versions p j = pj and p j = 1 give almost the same results, and this is consistent with the theory because when σ → 0, from [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF],

Run

pj → CDF -1 χ 2 (n) (∞) = 1.
Note that, on the other extreme, when the error on the gradient is big (σ ≫ 1), version p j = pj approaches version p j = p min since pj ≃ p min .

6.2. Expensive gradient case. Let us assume that, in practice, for a given problem, one has two routines for gradient calculation. The first routine computes the exact gradient and is expensive. The second routine is less expensive but computes only an approximation of the gradient. The model gradient results from a call to either routine. In this section, we propose a technique to choose the probability of calling the exact gradient which makes our approach applicable. Algorithm 6.1 (Algorithm to determine when to call the exact gradient g M j ).

Initialization

Choose the constant p min ∈ (0, 1) (p min is the lower bound of all the probabilities p * j ). For a chosen probability pj such that pj ≥ p min 1. Sample a random variable U ∼ U ([0, 1/p j ]), independently from F M j-1 , and U ([0, 1/p j ]) is the uniform distribution on the interval [0, 1/p j ].

1.1 If U ≤ 1, compute g M j using the routine which gives the exact gradient. 1.2 Otherwise, compute g M j using the routine which gives an approximation of the exact gradient. Proof. By using inclusion of events, we have that

p * j = P g M j -J(X j ) ⊤ F (X j ) ≤ κ eg Γ α j F M j-1 ≥ P g M j -J(X j ) ⊤ F (X j ) = 0 F M j-1
and from Algorithm 6.1 we conclude that

P g M j -J(X j ) ⊤ F (X j ) = 0 F M j-1 ≥ P (U ≤ 1) = 1 1/p j ,
and thus p * j ≥ pj . The other inequality, pj ≥ p min , is imposed in the algorithm. For the experiments we use the same test function and the same parameters as in section 6.1. In step 1.2 of Algorithm 6.1, we set the model gradient g M j to the exact gradient of the function plus a Gaussian noise sampled from N (0, 10I). Across all iterations, we use Algorithm 6.1 to compute g M j with the three following choices of pj :

• pj = 1/10, i.e., at iteration j the model gradient coincides with the exact gradient with probability at least pj = 1/10. Moreover, we have p * j ≥ pj , where pj is the same as in ( 9), and thus one can choose p j = max{1/10, pj }.

• pj = 1/50, with the same analysis as before and one can choose p j = max{1/50, pj }.

• pj ≃ 0 (p j = 10 -10 in the experiment below), i.e., at iteration j the probability that the model gradient coincides with the exact gradient is very small. Thus one can choose p j = pj . Figure 2 depicts the average of the function values and the absolute error of the iterates over 60 runs of Algorithm 3.1 when using the three choices of the probability p j . As expected, the better the quality of the model is the more efficient Algorithm 3.1 is (fewer iterations are needed to "converge" in the sense of sufficiently reducing the objective function value and absolute error). We can clearly see that Algorithm 3.1 using the models for which p j = max{1/10, pj } provides a better approximation to the minimizer of the objective function than using the models for which p j = max{1/50, pj }, and this latter one is better than the case when p j = pj .

7. Application to data assimilation. Data assimilation is the process by which observations of a real system are incorporated into a computer model (the forecast) to produce an estimate (the analysis) of the state of the system. 4DVAR is the data assimilation method mostly used in numerical weather prediction centers worldwide. 4DVAR attempts to reconcile a numerical model and the observations, by solving a very large weighted nonlinear least squares problem. The unknown is a vector of system states over discrete points in time. The objective function to be minimized is the sum of the squares of the differences between the initial state and a known background state at the initial time and the differences between the actual observations and ones predicted by the model. 7.1. 4DVAR problem. We want to determine x 0 , . . . , x T , where x i is an estimator of the state X i at time i, from the background state

X 0 = x b + W b , W b ∼ N (0, B). The observations are denoted by y i = H i (X i ) + V i , V i ∼ N (0, R i ), i = 0, . . .

, T , and the numerical model by

X i = M i (X i-1 ) + W i , W i ∼ N (0, Q i ), i = 1, . . . , T
, where M i is the model operator at time i and H i is the observation operator at time i (both not necessarily linear). The random vectors W b , V i , W i are the noises on the background, on the observation at time i, and on the model at time i, respectively, and are supposed to be Gaussian distributed with mean zero and covariance matrices B, R i , and Q i , respectively. Assuming that the errors (the background, the observation, and the model errors) are independent from each other and uncorrelated in time [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF], the posterior probability function of this system (in other words, the pdf of X 0 , . . . , X T knowing y 0 , . . . , y T ) is proportional to [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF] exp

-1 2 x 0 -x b 2 B -1 + T i=1 x i -M i (x i-1 ) 2 Q -1 i + T i=0 y i -H i (x i ) 2 R -1 i
and therefore the maximizer of the posterior probability function estimator is defined to be the minimizer of the weak constraint 4DVAR problem [START_REF] Trémolet | Model-error estimation in 4D-Var[END_REF] defined as the minimizer of the function defined in [START_REF] Bell | The iterated Kalman smoother as a Gauss-Newton method[END_REF], which is the negative logarithm of (10).

Incremental 4DVAR.

To find the solution of the nonlinear least squares problem (2), one proceeds iteratively by linearization. At each iteration, one solves the auxiliary linear least squares subproblem defined in (3) for the increments δx 0 , . . . , δx T . Such an iterative process is nothing else than the Gauss-Newton method [START_REF] Bell | The iterated Kalman smoother as a Gauss-Newton method[END_REF] for nonlinear least squares, known in the data assimilation community as the incremental approach [START_REF] Courtier | A strategy for operational implementation of 4D-VAR, using an incremental approach[END_REF].

Denote

z i = δx i , z b = x b -x 0 , z = [z 0 ; . . . ; z T ], m i = M i (x i-1 ) -x i , d i = y i -H i (x i ), M i = M ′ i (x i-1 )
, and 3) becomes [START_REF] Khare | An investigation into the application of an ensemble Kalman smoother to high-dimensional geophysical systems[END_REF] min

H i = H ′ i (x i ). Then (
z∈R n(T +1) 1 2 z 0 -z b 2 B -1 + T i=1 z i -M i z i-1 -m i 2 Q -1 i + T i=0 d i -H i z i 2 R -1 i .
It is known that the solution of the linear least squares problem [START_REF] Khare | An investigation into the application of an ensemble Kalman smoother to high-dimensional geophysical systems[END_REF] is exactly the same as the Kalman smoother estimator for the following linear system (see [START_REF] Tshimanga | Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation[END_REF])

Z 0 = z b + W b , W b ∼ N (0, B), (12)
Z i = M i Z i-1 + m i + W i , W i ∼ N (0, Q i ), i = 1, . . . , T, (13)
d i = H i Z i + V i , V i ∼ N (0, R i ), i = 0, . . . , T. ( 14
)
For simplicity, we now rewrite the linear system ( 12)-( 14) as

Z = Z b + W, W ∼ N (0, B W ), ( 15
)
D = HZ + V, V ∼ N (0, R), (16)
where

Z = [Z 0 ; . . . ; Z T ] is the joint state of the states Z 0 , . . . , Z T , D = [d 0 ; d 1 ; . . . ; d T ], Z b = [z b ; M 1 z b + m 1 ; M 2 (M 1 z b + m 1 ) + m 2 ; . . . ; M T (• • • M 1 z b + m 1 • • • ) + m T ], H = diag(H 0 , . . . , H T ) is the joint observation operator, W = [W b ; M 1 W b + W 1 ; M 2 (M 1 W b + W 1 ) + W 2 ; . . . ; M T (• • • M 1 W b + W 1 • • • ) + W T ], B W = cov(W ), V = [V 0 ; V 1 ; . . . ; V T ], and R = cov(V ).
To simplify it even more, we make the change of variables U = Z -Z b , and then ( 15)-( 16) becomes

U ∼ N (0, B W ), D -HZ b = HU + V, V ∼ N (0, R),
and the linear least squares problem [START_REF] Khare | An investigation into the application of an ensemble Kalman smoother to high-dimensional geophysical systems[END_REF] becomes (with z replaced by u + Z b ) min

u∈R n(T +1) 1 2 u 2 B -1 W + D -HZ b -Hu 2 R -1 . ( 17
)
To solve problem [START_REF] Moré | The Levenberg-Marquardt algorithm: Implementations and theory[END_REF], we propose to use the EnKS as a linear least squares solver instead of the Kalman smoother. The ensemble approach is naturally parallelizable over the ensemble members. Moreover, the proposed approach uses finite differences from the ensemble, and no tangent or adjoint operators are needed (i.e., the method is free of derivatives).

Kalman and EnKS.

The Kalman smoother gives the expectation and the covariance of the state U (equivalently Z) knowing the data D, in other words it calculates U a = E(U |D) and P a = cov(U |D), and is described by

U a = K(D -HZ b ), P a = (I -KH)B W , K = B W H ⊤ (HB W H ⊤ + R) -1 .
For Z one has Z a = E(Z|D) = Z b + K(D -HZ b ). In the data assimilation community, the vector U a (equivalently Z a ) is called the analysis and the matrix K is called the Kalman gain.

The EnKS [START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF][START_REF] Evensen | Data Assimilation: The Ensemble Kalman Filter[END_REF] consists of applying Monte Carlo to generate an ensemble following N (0, B W ) and then use its corresponding empirical covariance matrix instead of B W to approximate U a . Let us denote by k the ensemble members index, running over k = 1, . . . , N , where N is the ensemble size. We sample an ensemble Ũ k from N (0, B W ) by first sampling w k b according to N (0, B), w k 1 according to N (0, Q 1 ), . . . , w k T according to N (0, Q T ), and then by setting Ũ k as follows:

Ũ k 0 = w k b , Ũ k 1 = M 1 w k b + w k 1 , . . . , Ũ k T = M T (• • • M 1 w k b + w k 1 • • • ) + w k T . Let Ũ k = [ Ũ k 0 ; Ũ k 1 ; . . . ; Ũ k T ] and Ū = 1 N N k=1 Ũ k and B N = 1 N -1 N k=1 ( Ũ k -Ū )( Ũ k -Ū ) ⊤
be the empirical mean and covariance of the ensemble Ũ k , respectively. One has

B N = CC ⊤ , where C = 1 √ N -1 [ Ũ 1 -Ū, Ũ 2 -Ū, . . . , Ũ N -Ū ].
We then build the centered ensemble U k = Ũ k -Ū . Note that the empirical mean of the ensemble U k is equal to zero and that its empirical covariance matrix is B N . Now one generates the ensemble U k,a as follows:

U k,a = U k + K N (D -HZ b -V k ), (18)
where V k is sampled from N (0, R), and

K N = B N H ⊤ (HB N H ⊤ + R) -1 .
In practice, the empirical covariance matrix B N is never computed or stored since to compute the matrix products B N H ⊤ and HB N H ⊤ only matrix-vector products are needed:

B N H ⊤ = 1 N -1 N k=1 U k U k ⊤ H ⊤ = 1 N -1 N k=1 U k h ⊤ k , HB N H ⊤ = H 1 N -1 N k=1 U k U k ⊤ H ⊤ = 1 N -1 N k=1 h k h ⊤ k , K N = 1 N -1 N k=1 U k h ⊤ k 1 N -1 N k=1 h k h ⊤ k + R -1
,

where h k = HU k = [H 0 U k 0 ; . . . ; H T U k T ]
. We denote by Ū a and V the empirical mean of the ensembles U k,a and V k , respectively. One has from ( 18)

Ū a = K N (D -HZ b -V ). ( 19
)
It is known that when N → ∞, Ū a → U a in L p (see [START_REF] Le Gland | Large sample asymptotics for the ensemble Kalman filter[END_REF][START_REF] Mandel | On the convergence of the ensemble Kalman filter[END_REF]) and, thus, asymptotically, Ū a is the solution of the linearized subproblem (17) (and Ū a + Z b is the solution of the linearized subproblem (11)).

7.4. The linearized least squares subproblems arising in EnKS. From [START_REF] Osborne | Nonlinear least squares-the Levenberg algorithm revisited[END_REF] we conclude that Ū a is the Kalman smoother estimator for the following system,

Ũ ∼ N (0, B N ), D = H Ũ + Ṽ , Ṽ ∼ N (0, R), where D = D -HZ b -V . (20)
Hence, for a large N (such that B N is invertible), Ū a is the solution of the following linear least squares problem [START_REF] Tshimanga | Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation[END_REF] min

u∈R n(T +1) 1 2 u 2 (B N ) -1 + Hu -D 2 R -1 .
From the above derivation, we conclude that when we use the EnKS (until now with exact derivatives) to approximate the solution of the linearized subproblem [START_REF] Khare | An investigation into the application of an ensemble Kalman smoother to high-dimensional geophysical systems[END_REF], what is obtained is the solution of the linear least squares problem [START_REF] Tshimanga | Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation[END_REF]. The least squares model in ( 21) can be seen, in turn, as a realization of the following stochastic model,

(22) 1 2 u 2 B -1 + Hu -D 2 R -1 ,
where B -1 and D are random variables, with realizations (B N ) -1 and D, respectively. Both the incremental method and the method which approximates the solution of the linearized subproblem (11) using EnKS may diverge. Convergence to a stationary point of ( 2) can be recovered by controlling the size of the step, and one possibility to do so is to consider the application of the Levenberg-Marquardt method as in Algorithm 3.1. As in [START_REF] Mandel | Hybrid Levenberg Marquardt and weak constraint ensemble Kalman smoother method[END_REF], at each step, a regularization term is then added to the model in [START_REF] Tshimanga | Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation[END_REF],

(23) m(x + u) = 1 2 u 2 (B N ) -1 + Hu -D 2 R -1 + γ 2 u 2 ,
which corresponds to adding a regularization term to the model [START_REF] Zupanski | Maximum likelihood ensemble filter: Theoretical aspects[END_REF] (24)

M (x + u) = 1 2 u 2 B -1 + Hu -D 2 R -1 + Γ 2 u 2 .
We now provide the details about the solution of (23). For this purpose let (25)

P N = (I -K N H)B N .
Note that by using the Sherman-Morrison-Woodbury formula one has (26)

P N = (B N ) -1 + H T R -1 H -1 ,
in other words, P N is the inverse of the Hessian of model in [START_REF] Tshimanga | Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation[END_REF].

Proposition 7.1. The minimizer of the model (23

) is u * = Ū a -P N (P N + (1/γ 2 )I n ) -1 Ū a .
Proof. Since Ū a is the solution of problem ( 21), a Taylor expansion around Ū a of the model in [START_REF] Tshimanga | Limited-memory preconditioners, with application to incremental four-dimensional variational data assimilation[END_REF] gives

1 2 u 2 (B N ) -1 + Hu -D 2 R -1 = 1 2 Ū a 2 (B N ) -1 + H Ū a -D 2 R -1 + u -Ū a 2 (P N ) -1 .
Hence, the minimizer of the model ( 23) is the same as the minimizer of

1 2 Ū a 2 (B N ) -1 + H Ū a -D 2 R -1 + u -Ū a 2 (P N ) -1 + γ 2 u 2
and thus given by

u * = (P N ) -1 + γ 2 I -1 (P N ) -1 Ū a . (27)
By using the Sherman-Morrison-Woodbury formula, one has

(P N ) -1 + γ 2 I -1 = P N -P N P N + (1/γ 2 )I n -1 P N ,
which together with (27) concludes the proof.

7.5. Derivative-free LM-EnKS. The linearized model (LM) and observation operators appear only when acting on a given vector, and therefore they could be efficiently approximated by finite differences. The linearized observation operator H i = H ′ i (x i ) appears in the action on the ensemble members and can be approximated by

H i δx i = H ′ i (x i )δx i ≃ H i (x i + τ δx i ) -H i (x i ) τ ,
where τ > 0 is a finite differences parameter. Originally, in EnKS, to avoid the derivatives of H i , the quantity

H i δx i = H i (x i +δx i -x i ) is approximated by H i (x i +δx i )-H i (x i )
, which is equivalent to using finite differences with the parameter τ = 1. The LM M 1 = M ′ 1 (x 0 ) appears in the action on a given vector (in this case z b ), and so do the remaining ones M 3 , . . . , M T . Such actions can be approximated by finite differences in the following way:

M 1 z b = M ′ 1 (x 0 )z b ≃ M 1 (x 0 + τ z b ) -M 1 (x 0 ) τ , M 2 (M 1 z b + m 1 ) = M ′ 2 (x 1 )(M 1 z b + m 1 ) ≃ M 2 (x 1 + τ (M 1 z b + m 1 )) -M 2 (x 1 ) τ ≃ M 2 (x 1 + M 1 (x 0 + τ z b ) -M 1 (x 0 ) + τ m 1 ) -M 2 (x 1 ) τ .
Since our approach is derivative free, we replace all the derivatives of the model and of the observation operators by approximation by finite differences. The quantities using derivatives

then become ĥk = H 0 x 0 + τ U k 0 -H 0 (x 0 ) τ ; . . . ; H T x T + τ U k T -H T (x T ) τ ≃ h k , KN = 1 N -1 N k=1 U k ĥ⊤ k 1 N -1 N k=1 ĥk ĥ⊤ k + R -1 ≃ K N , (28) Ẑb = z b ; M 1 (x 0 + τ z b ) -M 1 (x 0 ) τ + m 1 ; . . . ≃ Z b , (29) ĤZ b = H 0 (x 0 + τ z b ) -H 0 (x 0 ) τ ; . . . ≃ HZ b , Û a = KN (D -ĤZ b -V ) ≃ Ū a , P N = B N -KN 1 N -1 N k=1 ĥk U k ⊤ ≃ P N , (30) û * = Û a -P N P N + (1/γ 2 )I n -1 Û a ≃ u * . (31)
Since û * is an approximation to u * using finite differences for derivatives, there exists a constant M > 0, which depends on the second derivatives of the model and observation operators, such that e ≤ M τ , where e = u * -û * . Moreover, the minimizer u * of the weighted least squares model ( 23) is the solution of the normal equations

B N -1 + H T R -1 H + γ 2 I u * = H T R -1 D,
where H T R -1 D = ∇m(x) = g m , and thus

B N -1 + H T R -1 H + γ 2 I û * = g m -B N -1 + H T R -1 H + γ 2 I e,
and so û * can be seen as an inexact solution of the normal equations, with a residual equal to

r = -(B N ) -1 + H T R -1 H + γ 2 I e.
We have seen that the solution of the normal equations can be inexact as long as Assumption 4.2 is met. The residual r is then required to satisfy r ≤ ǫ g m for some ǫ > 0, to fulfill the global convergence requirements of our Levenberg-Marquardt approach, and for this purpose we need the following assumption. The Jacobian of the observation operator H is uniformly bounded, i.e., there exists κ H > 0 such that H ′ i (x i ) ≤ κ H for all i ∈ {0, . . . , T } and for all iterations j. We note that the iteration index j has been omitted from the notation of this section until now. In fact, the point x has been denoting the iterate x j . Proposition 7.2. Under Assumption 7.1, if the finite differences parameter τ is such that

(32) τ ≤ ǫ g m M (B N ) -1 + κ 2 H R -1 + γ 2 ,
then r ≤ ǫ g m .

Proof. One has

r ≤ (B N ) -1 + H T R -1 H + γ 2 I e ≤ (B N ) -1 + κ 2 H R -1 + γ 2 M τ ≤ ǫ g m .
Now, from (24) the gradient of the stochastic model is g M j = -H ⊤ R -1 D and from ( 17) the exact gradient of the function to be minimized in problem [START_REF] Bell | The iterated Kalman smoother as a Gauss-Newton method[END_REF] is -H ⊤ R -1 (D -HZ b ). Thus,

p * j = P H ⊤ R -1 (D -HZ b -D) ≤ κ eg Γ α j F M j-1 .
But we know that D -HZ b -D = V = (1/N ) N i=1 V i , where V i are independently and identically distributed and follow N (0, R), and thus D -

HZ b -D ∼ N (0, R/N ) and R -1 (D - HZ b -D) ∼ R -1/2 √ N N (0, I). Thus p * j ≥ P κ H R -1/2 √ N N (0, I) ≤ κ eg Γ α j F M j-1 = P N (0, I) ≤ κ √ N Γ α j F M j-1 ,
where κ = κeg κ H R -1/2 . Since Γ j ≤ min{λ j γ 0 , γ max },

(33) p * j ≥ CDF -1 χ 2 (m)   κ √ N min{λ j γ 0 , γ max } α 2   def = pj ,
where m = T i=0 m i , m i is the size of y i , and γ max is the tolerance used in the stopping criterion. Note that lim N →∞ pj = 1, thus lim N →∞ p * j = 1, and hence when N → ∞ the gradient approximation using ensemble converges almost surely to the exact gradient.

We are now ready to propose a version of Algorithm 3.1 for the solution of the 4DVAR problem (2) when using EnKS as the linear solver. Algorithm 7.1 (Levenberg-Marquardt method based on probabilistic gradient models for data assimilation).

Initialization

Choose the constants η 1 ∈ (0, 1), η 2 , γ min , γ max > 0, and λ > 1. Select x 0 and γ 0 ∈ [γ min , γ max ]. Choose all the parameters related to solving the 4DVAR problem (2) using EnKS as the linear solver.

For j = 0, 1, 2, . . . and while γ j ≤ γ max 1. Let x = x j . Choose τ satisfying (32). Compute the increment û * using (31) and set z * = û * + Ẑb , where Ẑb is computed as in (29). Let s j = z * . 2. Compute ρ j = f (x j )-f (x j +s j ) m j (x j )-m j (x j +s j ) , where f is the nonlinear least squares model in [START_REF] Bell | The iterated Kalman smoother as a Gauss-Newton method[END_REF] and m j is the model (23). 3. If ρ j ≥ η 1 , then set x j+1 = x j + s j and

γ j+1 =      λγ j if g m j < η 2 /γ 2 j , max γ j λ 1-p j p j , γ min if g m j ≥ η 2 /γ 2 j ,
where p j = pj is computed as in (33).

Otherwise, set x j+1 = x j and γ j+1 = λγ j .

7.6. Derivative-free LM-EnKS in practice. In sections 7.4-7.5, we have assumed that the ensemble size N was large enough for the empirical covariance matrix B N to be invertible (holding Proposition 7.1). However, for the EnKS to be relevant the ensemble size has to be smaller than the dimension of the state space. In this subsection, we explain how we circumvent this problem in practice. The theoretical extension of our method to small ensemble sizes as well as its performance for large and realistic problems is the subject of future research.

For small values of the ensemble size (N < n), the matrix B N is no longer invertible. In particular, the Sherman-Morrisson-Woodbury formula is no longer applicable, as it was before to establish (25)-( 26) in terms of (B N ) -1 . However, following the spirit of (26), we could think of using pseudoinverses instead, and approximating the matrix P N defined in (25) by

P N = (B N ) † + H T R -1 H † .
In practice what we do is simply to replace inverses by pseudoinverses in all calculations, namely, in (28) and in (31).

Another concern when using a small ensemble size is how to ensure Assumption 3.1 (gradient model being (p j )-probabilistically accurate). When N is sufficiently large, we have shown that the formula (33) provides a value of p j that satisfies the assumption. This formula can, however, still be used in practice for small values of N . 7.7. Computational experiments with Lorenz-63 model. To evaluate the performance of Algorithm 7.1 for data assimilation, we will test it using the classical twin experiment technique used in the data assimilation community. This technique consists on fixing an initial true state (denoted by truth 0 ) and then to integrate it over time using the model to obtain the true state at each time i (denoted by truth i ). We then build the data y i by applying the observation operator H i to the truth at time i and by adding a Gaussian perturbation N (0, R i ). Similarly, the background x b is sampled from the Gaussian distribution with mean truth 0 and covariance matrix B. Then we try to recover the truth using the observations and the background.

For the 4DVAR problem (2), we consider the Lorenz-63 model, a simple dynamical system with chaotic behavior. The Lorenz equations are given by the nonlinear system dx dt = -σ(x -y), dy dt = ρx -y -xz, and dz dt = xy -βz, where x = x(t), y = y(t), z = z(t), and σ, ρ, β are parameters. The state at time t is X t = (x(t), y(t), z(t)) ⊤ ∈ R 3 . This nonlinear system is discretized using a fourth order Runge-Kutta method. The parameters σ, ρ, β are chosen as 10, 28, and 8/3, respectively. The initial truth is set to (1, 1, 1) ⊤ and the truth at time i to truth i = M(truth i-1 ) + W i , where W i is sampled from N (0, Q i ) and M is the model obtained by discretization of the Lorenz-63 model. The model error covariance is given by Q i = σ 2 q I, where σ q = 10 -4 . The background mean x b is sampled from N (truth 0 , B). The background covariance is B = σ 2 b I, where σ b = 1. The time step is chosen as dt = 0.11 > 0.01. (Note here that the model at time t + 1, as a function of the model at time t, becomes more nonlinear as dt increases, and this justifies having chosen dt larger than in [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF].) The time windows length is T = 40. The observation operator is H i = 10I. At each time i, the observations are constructed as follows:

y i = H i (truth i ) + V i , where V i is sampled from N (0, R), R = σ 2
r I, and σ r = 1. Following the spirit of Assumption 4.2, the finite difference parameter is set as

τ j = min    10 -3 , ǫ j g m j M (B N ) † + κ 2 H R -1 + γ 2 j (B N ) †    ,
where the value of 1 is chosen for the unknown constants M and κ H (see Assumption 7.1). In this experimental framework, the model gradient is given by g m j = -H ⊤ R -1 D = 10 D, where D is computed according to [START_REF] Trémolet | Model-error estimation in 4D-Var[END_REF]. Then, following the spirit of Assumption 4.2, ǫ j is chosen as

ǫ j = min θ in γ α j , β in γ 2 j κ 2 Jm + γ 2 j ,
where β in = 1/2, θ in = 1, and α = 0.5. The unknown constant κ Jm (see Assumption 5.1) is set to 1.

The basic algorithmic parameters are set to η 1 = η 2 = 10 -6 , γ min = 10 -5 , γ max = 10 6 , and λ = 8. The initial regularization parameter is γ 0 = 1. Finally, we set κ = 1 in the calculation of pj given in (33).

Figure 3 depicts the plots of the objective function values for one run of Algorithm 7.1, using the choices p j = pj and p j = 1 and four ensemble sizes N = 4, 40, 80, 400. A single run shows well the behavior of the algorithm on this problem, thus there is no need to take averages over several runs. For all ensemble sizes used, the version using p j = 1 stagnated after some iterations, and could not approximate the minimizer with a decent accuracy. One can see that the version with p j = pj performs much better than the one with p j = 1, regardless of the size of the ensemble. As expected, the larger this size, the better is the accuracy of the final solution found (which can be further confirmed in Table 2). These results illustrate the importance of using probability p j = pj to update the regularization parameter γ.

Table 2

The table shows the final values of the objective function found for the two versions pj = 1 and pj = pj and the four ensemble sizes. 8. Conclusions. In this paper we have adapted the Levenberg-Marquardt method for nonlinear least squares problems to handle the cases where the gradient of the objective function is subject to noise or only computed accurately within a certain probability. The gradient model was then considered random in the sense of being a realization of a random variable, and assumed first order accurate under some probability p * j (see [START_REF] Conn | Trust-Region Methods[END_REF]). Given the knowledge of a lower bound p j for this probability (see Assumption 3.1), we have shown how to update the regularization parameter of the method in such a way that the whole approach is almost surely globally convergent. The analysis followed similar steps as in the theory in [START_REF] Bandeira | Convergence of trust-region methods based on probabilistic models[END_REF]. The main difficulty in the application of Algorithm 3.1 is to ensure that the models are indeed (p j )-probabilistically accurate, but we presented a number of practical situations where this is achievable.

The last section of the paper was devoted to the well-known 4DVAR problem in data assimilation. We have shown that a lower bound for the probability of first order accuracy can also be provided here (to be used in our Levenberg-Marquardt framework) when using the EnKS method for the formulation and solution of the corresponding linearized least squares subproblems.

We have also covered the situation where the linearized least squares problems arising in the Levenberg-Marquardt method are solved inexactly, which then encompasses a range of practical situations, from inexactness in linear algebra to inexactness in derivatives. This is particularly useful in the 4DVAR application to accommodate finite differences of the nonlinear operators involved.

A number of issues need further and deeper investigation, in particular, the study of the performance of our approach when applied to large and realistic data assimilation problems.

After we had submitted our paper, Bocquet and Sakov [START_REF] Bocquet | An iterative ensemble Kalman smoother[END_REF] extended their previous approach [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF] to 4DVAR and used finite difference approximations for the tangent operators, similarly to our paper. Bocquet and Sakov [START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF][START_REF] Bocquet | An iterative ensemble Kalman smoother[END_REF] nest the minimization loop for the 4DVAR objective function inside the EnKS and minimize over the span of the ensemble, rather than nesting EnKS as a linear solver inside the 4DVAR minimization loop over the full state space, as we do. Moreover, they use a classical version of the Levenberg-Marquardt method to perform their minimization without any control or assumption on the derivative approximations arising from the use of ensembles. Their method was designed for strong-constraint 4DVAR, i.e., for the case Q i = 0 ∀i.

Appendix.

Proof of Lemma 4.1. In the proof we will omit the indices j. One has m(x) -m(x + s in ) = -g ⊤ m s in -min{1, γ 2-α min }γ 2+α

. If s = s cg is obtained by truncated CG, then there exists an orthogonal matrix V with a first column given by -g m / g m and such that

s cg = V V ⊤ (J ⊤ m J m + γ 2 I)V -1 V T g m = V V ⊤ J ⊤ m J m V + γ 2 I -1 g m e 1 ,
where e 1 is the first vector of the canonical basis of R n . From the positive semidefiniteness of V ⊤ J ⊤ m J m V , we immediately obtain s cg ≤ g m /γ 2 . To prove the second inequality we apply the Sherman-Morrisson-Woodbury formula, to obtain

s cg = V 1 γ 2 I - 1 γ 4 (J m V ) ⊤ I + (J m V )(J m V ) ⊤ γ 2 -1 (J m V ) g m e 1 .
Since V e 1 = -g m / g m ,

γ 2 s cg + g m = - 1 γ 2 V (J m V ) ⊤ I + (J m V )(J m V ) ⊤ γ 2 -1 (J m V ) g m e 1 .
Now, from the fact that (J m V )(J m V ) ⊤ /γ 2 is positive semidefinite, the norm of the inverse of I + (J m V )(J m V ) ⊤ /γ 2 is no greater than one, and thus (since V is orthogonal)

γ 2 s cg + g m ≤ J m 2 g m

Lemma 4 . 1 .

 41 Under Assumption 4.2, an inexact step s in j of the form[START_REF] Courtier | A strategy for operational implementation of 4D-VAR, using an incremental approach[END_REF] achieves Cauchy decrease if it satisfies Assumption 4.1 with θ f cd = 2(1 -β in ).

  (a) Average of function values. (b) Average of absolute error of iterates. (c) Average percentage of successful iterations.

Figure 1 .

 1 Figure 1. Average results of Algorithm 3.1 for 60 runs when using probabilities pj = 1 (dotted line), pj = pj (solid line), and pj = pmin (dashed line). The x-axis represents number of iterations.

Lemma 6 . 1 .

 61 If we use Algorithm 6.1 to compute the model gradient at the jth iteration of Algorithm 3.1, then we have p * j ≥ pj ≥ p min . (a) Average of function values. (b) Average of absolute error of iterates.

Figure 2 .

 2 Figure 2. Average results of Algorithm 3.1 for 60 runs when using probabilities pj = pj (solid line), pj = max{1/10, pj} (dotted line), and pj = max{1/50, pj} (dashed line). The x-axis represents number of iterations.

Assumption 7 . 1 .

 71 The approximation û * of u * satisfies e ≤ M τ , where e = u * -û * , for some constant M > 0.

  (a) Objective function values (size N = 4). (b) Objective function values (size N = 40). (c) Objective function values (size N = 80). (d) Objective function values (size N = 400).

Figure 3 .

 3 Figure 3. Results of one run of Algorithm 7.1, using probabilities pj = 1 (dotted line) and pj = pj (solid line), for different ensemble sizes. The x-axis represents number of iterations.

γ 2 .m 2 g m 2 γ 4 ≤ 4 J m 2 g m 2 + 2θ in g m 2 min{1, γ 2

 2244222 Finally (recalling γ ≥ γ min ),|(s cg ) ⊤ (γ 2 (s cg ) + g m )| ≤ s cg γ 2 s cg + g m ≤ J -α min }γ 2+α

The research of the third author was supported by FCT under grants PTDC/MAT/116736/2010 and PEst-C/MAT/UI0324/2011. †