Serge Gratton
email: serge.gratton@enseeiht.fr

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

IMPROVED BOUNDS FOR SMALL-SAMPLE ESTIMATION *

Keywords: statistical estimation, small-sample estimation, norm estimation, condition number estimation

We derive improved error bounds for small-sample statistical estimation of the matrix Frobenius norm. The bounds rigorously establish that small-sample estimators provide reliable orderof-magnitude estimates of norms and condition numbers, for matrices of arbitrary rank, even when very few random samples are used.

1. Introduction. The Frobenius norm of a matrix B ∈ R m×n is in principle straightforward to compute. In many applications, however, one is interested in the Frobenius norm of a matrix whose entries are either too expensive to compute explicitly or not available at all. For example, B might be the inverse of a known matrix A of large dimensions, and explicitly computing a factorization of A that would reveal A -1 might be too expensive. Another example is condition number estimation, in which the matrix B is the (unknown) matrix representation of the Fréchet derivative of a given matrix function; see, e.g., [START_REF] Laub | Applications of statistical condition estimation to the solution of linear systems[END_REF][START_REF] Laub | Fast condition estimation for a class of structured eigenvalue problems[END_REF][START_REF] Laub | Statistical condition estimation for the roots of polynomials[END_REF][START_REF] Baboulin | Statistical estimates for the conditioning of linear least squares problems[END_REF][START_REF] Gratton | Stochastic conditioning of matrix functions[END_REF][START_REF] Baboulin | Computing least squares condition numbers on hybrid multicore/GPU systems[END_REF][START_REF] Diao | Structured condition numbers of structured Tikhonov regularization problem and their estimations[END_REF][START_REF] Diao | Small sample statistical condition estimation for the total least squares problem[END_REF], to cite only a few examples. For a more thorough discussion of other applications, we refer to [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF] and the references therein.

In this paper we investigate some properties of two randomized estimators of the Frobenius norm of a matrix B ∈ R m×n . The estimators only require B through the action of a matrix-vector product: u ← Bv. We are interested in "small-sample" estimation as defined by Gudmundsson, Kenney, and Laub [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF][START_REF] Kenney | Small-sample statistical condition estimates for general matrix functions[END_REF]. Specifically, we wish to prove that, with high probability, a realization of the estimator is close in magnitude to B F even when very few random samples are used. In other words,

(1.1) Prob B F τ ≤ estimator k (B) ≤ τ B F ≥ 1 -δ(k, τ ),
where δ(k, τ ) → 0 as the number of samples k and/or the tolerance τ > 1 increases, and δ(k, τ ) is close to 0 even for small k and for τ not too large. We present error bounds that do not deteriorate as the rank of B increases, thereby settling a conjecture from [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF]. Our results also lead to an improvement of probabilistic relative error bounds from [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF] and [START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF].

2. The Gaussian and GKL estimators. First we recall the Gaussian estimator, as well as the small-sample statistical estimator of Gudmundsson, Kenney, and Laub [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF][START_REF] Kenney | Small-sample statistical condition estimates for general matrix functions[END_REF], which we abbreviate as the GKL estimator. Definition 2.1. Let Z be an n × k matrix with k ≤ n whose entries are mutually independent standard normal random variables, and let Q be the n×k "thin" Q factor in the QR decomposition of Z. The Gaussian estimator of B F is defined as

(2.1) ψ k (B) = 1 √ k BZ F ,
while the GKL estimator is defined as

(2.2) η k (B) = √ n √ k BQ F .
The GKL estimator has been widely used for norm and condition number estimation. It is shown in [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF] that if B has rank r,

(2.3) Prob B F τ ≤ η k (B) ≤ τ B F ≥ 1 -r 1 -I( kτ 2 n ) + I( k nτ 2 ) ,
where

I(α) = 1 if α ≥ 1 and
I(α) = α 0 t (k-2)/2 (1 -t) (n-k-2)/2 dt 1 0 t (k-2)/2 (1 -t) (n-k-2)/2 dt is a regularized Beta function if 0 < α < 1.
When the rank of B is close to 1, the right-hand side of (2.3) increases fast toward 1 as k and/or τ increase. On the other hand, when r is large, the bound becomes very pessimistic. Furthermore, because r is not generally known a priori, for prediction purposes the best one can do is use (2.3) with r = min{m, n}, leading to a very pessimistic bound. It is conjectured in [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF] that the bound (2.3) holds without the factor r. It is this conjecture which motivated the present study.

Next we turn to the Gaussian estimator in (2.1). It is straightforward to verify that ψ 2 k (B) is a linear combination of independent chi-squared variables with k degrees of freedom, denoted χ 2 i (k) below:

(2.4)

ψ 2 k (B) = 1 k r i=1 σ 2 i χ 2 i (k),
where σ 1 ≥ • • • ≥ σ r > 0 are the nonzero singular values of B. From this we have

E{ψ 2 k (B)} = B 2 F , Var{ψ 2 k (B)} = 2 k r i=1 σ 4 i ≤ 2 k B 4 F .
Recall that for small-sample estimation, we are restricted to k ≪ n and we are interested in order-of-magnitude estimation. The following proposition establishes that in this setting the two estimators η k (B) and ψ k (B) are essentially equivalent. Proposition 2.2. For any matrix B ∈ R m×n and any τ > 0, the estimators ψ k (B) and η k (B) defined in (2.1) and (2.2) satisfy

Prob 1 - √ k + τ √ n ≤ ψ k (B) η k (B) ≤ 1 + √ k + τ √ n ≥ 1 -2 exp - τ 2 2 . Proof. Because Z = QR with Q T Q = I k , ψ k (B) = √ n √ k BQR F √ n ≤ √ n √ k BQ F R 2 √ n = η k (B) R 2 √ n = η k (B) Z 2 √ n .
On the other hand,

η k (B) = √ n √ k BQRR -1 F ≤ √ n √ k BZ F R -1 2 = ψ k (B) √ n σ min (Z) ,
where σ min (•) denotes the smallest singular value. It follows that

σ min (Z) √ n ≤ ψ k (B) η k (B) ≤ Z 2 √ n .
To complete the proof we use a result of Gordon [START_REF] Gordon | Some inequalities for Gaussian processes and applications[END_REF][START_REF] Gordon | Gaussian and almost spherical sections of convex bodies[END_REF] popularized by Davidson and Szarek [6]:

for any τ > 0, √ n - √ k -τ ≤ σ min (Z) ≤ Z 2 ≤ √ n + √ k + τ holds with probability at least 1 -2 exp(-τ 2 /2).
Proposition 2.2 is not very surprising: the fact that Gaussian and random orthonormal matrices behave the same way is well known and dates at least as far back as [START_REF] Marcus | Random Fourier Series with Applications to Harmonic Analysis[END_REF].

In the following section we present new bounds of the form (1.1) for ψ k (B). By Proposition 2.2, these bounds also essentially apply for η k (B) when k ≪ n. Our bounds are stated in terms of the stable rank of B,

(2.5) ρ(B) = B 2 F B 2 2 .
The concept of stable rank dates at least as far back as [START_REF] Bourgain | Invertibility of 'large' submatrices with applications to the geometry of Banach spaces and harmonic analysis[END_REF][START_REF] Bourgain | On a problem of Kadison and Singer[END_REF]; see, for example, [START_REF] Rudelson | Sampling from large matrices: An approach through geometric functional analysis[END_REF][START_REF] Tropp | An introduction to matrix concentration inequalities, Found[END_REF] for many examples of its use in linear algebra. In contrast to (2.3), our bounds are useful even when the matrix B has a large rank: they do not deteriorate with increasing r but improve with increasing ρ.

3. Improved convergence bounds.

3.1. Chernoff bounds. The Gaussian estimator ψ k (B) is one of several randomized estimators surveyed in [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF] and [START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF]. In our notation, [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF]Theorem 5.2] shows that for ǫ < 0.1,

(3.1) Prob ψ 2 k (B) -B 2 F ≤ ǫ B 2 F ≥ 1 -2 exp - k 20 ǫ 2 .
Thus, for any ǫ ∈ (0, 0.1) and δ ∈ (0, 1), at least k = 20ǫ -2 ln(2/δ) samples are required to ensure that

Prob ψ 2 k (B) -B 2 F ≤ ǫ B 2 F ≥ 1 -δ.
In [START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF]Theorem 3], this result is improved to

(3.2) k = 8ρ -1 ǫ -2 ln(2/δ).
The value of k in (3.2) can be very large when ǫ is small because of the ǫ -2 factor. Recall, however, that in the setting of small-sample estimation one does not require B F with a great deal of accuracy and one is limited to using a small value of k.

The following theorem applies in this setting.

Theorem 3.1. Let B ∈ R m×n be a matrix of rank r and stable rank ρ whose Frobenius norm is to be estimated, and let ψ k (B) be defined by (2.1). For any τ > 1 and k ≤ n,

Prob{ψ k (B) ≥ τ B F } ≤ exp - kρ 2 τ -1 2 , (3.3) Prob ψ k (B) ≤ B F τ ≤ exp - kρ 4 (τ 2 -1) 2 τ 4 , (3.4) and Prob B F τ < ψ k (B) < τ B F ≥ 1 -exp - kρ 2 τ -1 2 -exp - kρ 4 (τ 2 -1) 2 τ 4 . (3.5)
Proof. We use the same idea as in [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF] and [START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF]. Namely, we use the fact that

Prob B F τ < ψ k (B) < τ B F = 1 -Prob{ψ k (B) ≥ τ B F } -Prob{ψ k (B) ≤ B F /τ } (3.6)
and bound each probability in the right-hand side separately by using a Chernoff approach. For any t > 0, by Markov's inequality we have

Prob{ψ k (B) ≥ τ B F } = Prob exp(tkψ 2 k (B)) ≥ exp(tkτ 2 B 2 F ) ≤ E exp(tkψ 2 k (B)) exp(tkτ 2 B 2 F ) . (3.7)
In the above, the numerator is the moment generating function (MGF) of the random variable kψ 2 k (B). Recall from (2.4) that this is a linear combination of independent chi-squared variables with k degrees of freedom. It is known (see, e.g., [18, section 3.2]) that its MGF is defined for |t| < 1 2σ 2 1 and satisfies

E exp(tkψ 2 k (B)) = r i=1 (1 -2σ 2 i t) -k/2 = exp - k 2 r i=1 ln(1 -2σ 2 i t) = exp   k 2 r i=1 ∞ j=1 (2σ 2 i t) j j   ,
where

σ 1 ≥ • • • ≥ σ r > 0 are the nonzero singular values of B. Pick t = µ 2σ 2 1
for some µ ∈ (0, 1). Then from (3.7),

Prob{ψ k (B) ≥ τ B F } ≤ exp   k 2 r i=1 ∞ j=1 µ j σ 2j i jσ 2j 1   exp -µkρτ 2 2 = exp   - kρ 2 µ(τ 2 -1) + k 2 r i=1 ∞ j=2 µ j σ 2j i jσ 2j 1   .
In the above,

r i=1 ∞ j=2 µ j σ 2j i jσ 2j 1 ≤ r i=1 ∞ j=2 µ j σ 2j i σ 2j 1 = r i=1 µ 2 σ 4 i /σ 4 1 1 -µσ 2 i /σ 2 1 ≤ µ 2 ρ 1 -µ , which gives Prob{ψ k (B) ≥ τ B F } ≤ exp - kρ 2 µ(τ 2 -1) - µ 2 1 -µ .
It is straightforward to verify that over all µ ∈ (0, 1), the smallest bound is achieved when µ = 1 -τ -1 , leading to (3.3).

Similarly, for any t ∈ 0, 1/(2σ 2 1 ) we have

Prob{ψ k (B) ≤ B F /τ } = Prob exp(-tkψ 2 k (B)) ≥ exp(-tk B 2 F /τ 2 ) ≤ E exp(-tkψ 2 k (B)) exp(tk B 2 F /τ 2 ) = exp   k 2 r i=1 ∞ j=1 (-1) j (2σ 2 i t) j j + tk B 2 F τ 2   . Picking t = µ 2σ 2 1
for some µ ∈ (0, 1), we obtain

Prob{ψ k (B) ≤ B F /τ } ≤ exp   - kρ 2 µ(1 -τ -2 ) + k 2 r i=1 ∞ j=2 (-1) j µ j σ 2j i jσ 2j 1   .
In the above,

r i=1 ∞ j=2 (-1) j µ j σ 2j i jσ 2j 1 ≤ r i=1 µ 2 σ 4 i 2σ 4 1 ≤ µ 2 ρ 2 , leading to Prob{ψ k (B) ≤ B F /τ } ≤ exp - kρ 2 µ(1 -τ -2 ) - µ 2 2 .
It is straightforward to verify that µ = 1 -τ -2 gives the optimal bound (3.4). Substituting (3.3) and (3.4) into (3.6) completes the proof.

We can use Theorem 3.1 to determine the number of samples k required to estimate B F to within a given factor τ with probability at least 1 -δ. Corollary 3.2. In the notation of Theorem 3.1, for any τ ≥ 2 and δ ∈ (0, 1), at least

k = 4ρ -1 τ 4 (τ 2 -1) 2 ln 2 δ samples are required to ensure that Prob B F τ < ψ k (B) < τ B F ≥ 1 -δ.
Copyright © by SIAM. Unauthorized reproduction of this article is prohibited.

Proof. If τ ≥ 2, then (τ -1) 2 ≥ (τ 2 -1) 2 2τ 4 . Thus from (3.5) Prob B F τ < ψ k (B) < τ B F ≥ 1 -2 exp - kρ 4 (τ 2 -1) 2 τ 4 .
The result follows by setting the above equal to 1 -δ and solving for k.

Theorem 3.1 can also be used to improve the constant factor in known relative error bounds.

Corollary 3.3. In the notation of Theorem 3.1, for any ǫ ∈ (0, 1 2 ),

Prob ψ k (B) -B F ≤ ǫ B F ≥ 1 -2 exp - kρ 2 ǫ 2 .
Therefore, for any ǫ ∈ (0, 1 2 ) and δ ∈ (0, 1), at least

(3.8) k = 2ρ -1 ǫ -2 ln 2 δ
samples are required to ensure that

Prob ψ k (B) -B F ≤ ǫ B F ≥ 1 -δ.
Proof. Note that

ψ k (B) -B F ≤ ǫ B F ⇐⇒ (1 -ǫ) B F ≤ ψ k (B) ≤ (1 + ǫ) B F . Taking τ = 1 + ǫ in (3.3) and τ = 1/(1 -ǫ) in (3.4) we obtain Prob ψ k (B) -B F ≤ ǫ B F ≥ 1 -exp - kρ 2 ǫ 2 -exp - kρ 2 ǫ 2 (2 -ǫ) 2 2 ≥ 1 -2 exp - kρ 2 ǫ 2 .
The result follows by setting the above equal to 1 -δ and solving for k.

As with (3.2), the bound (3.8) on the number of samples inevitably depends on ǫ -2 . However, the factor 8 in (3.2) is reduced to a factor 2. In other words, four times fewer samples than predicted by (3.2) are in fact required.

A combined bound.

We now return to the result given in Theorem 3.1. For fixed k and τ , the bound (3.5) increases very fast toward 1 as the stable rank ρ increases. For example, suppose that B ∈ R m×n with n = r = 1000 and B F = 10 B 2 . Even with small values of k = τ = 2, the bound states that

Prob B F 2 < ψ 2 (B) < 2 B F ≥ 1 -10 -12 .
This is much closer to 1 than the lower bound (2.3), even if, as conjectured in [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF], (2.3) holds without the factor r. (In this example, the lower bound (2.3) is negative, while ignoring the factor r in (2.3) gives 0.7610 to four digits.) On the other hand, when ρ is close to 1, the two exponential terms (3.3) and (3.4) exhibit different behavior. The bound (3.3) on the probability of overestimation remains quite small, even for small values of k and τ . For example, in the extreme case

B F = B 2 , with k = τ = 3, we have exp - kρ 2 τ -1 2 = exp(-6) ≤ 0.0025.
Clearly this term further decreases with increasing k and τ . The bound (3.4) on the probability of underestimation may be quite large for small k when ρ is small. Once again in the extreme case

B F = B 2 , with k = τ = 3, we have exp - k 2 B 2 F B 2 2 (τ 2 -1) 2 2τ 4 = exp - 16 27 ≈ 0.5529.
Furthermore, although this term is a decreasing function of k and of τ , it does not tend to 0 with increasing τ . Consequently, if ρ is close to 1, the bound (3.5) may be quite pessimistic when k is not very large and τ is very large, i.e., if we are interested in a very rough estimate of B F using very few samples. The following theorem extends Theorem 3.1 by providing a bound that is useful when ρ is close to 1 and/or τ is very large.

Theorem 3.4. In the notation of Theorem 3.1,

Prob B F τ < ψ k (B) < τ B F ≥ 1 -exp - kρ 2 (τ -1) 2 -ω(B, k , τ ), (3.9)
where

ω(B, k, τ ) = min exp - kρ 4 (τ 2 -1) 2 τ 4 , P kρ τ 2

and

(3.10)

P (α) = α/2 0 t (k-2)/2 exp(-t) dt ∞ 0 t (k-2)/2 exp(-t) dt is a regularized gamma function.
Proof. We only modify the bound (3.4) on underestimation from Theorem 3.1. Recall from (2.4) that

kψ 2 k (B) = r i=1 σ 2 i χ 2 i (k),
where the random variables χ 2 i (k) are mutually independent. Therefore,

Prob ψ k (B) ≤ B F τ = Prob r i=1 σ 2 i χ 2 i (k) ≤ k B 2 F τ 2 ≤ Prob r i=1 σ 2 i χ 2 i (k) ≤ k B 2 F τ 2 = r i=1 Prob χ 2 i (k) ≤ k B 2 F τ 2 σ 2 i = r i=1 P k B 2 F τ 2 σ 2 i ≤ P kρ τ 2 . (3.11)
The rest of the proof is as in Theorem 3.1. Since we have two upper bounds (3.4) and (3.11) for the probability of underestimation, the mininimum of the two is also a valid bound.

In contrast to (3.4), the bound (3.11) decreases toward 0 as τ increases, and is an increasing function of ρ. Thus, (3.4) is more suitable when ρ is large, while (3.11) is preferable when ρ is close to 1 or τ is very large.

4. Numerical experiments. We give some numerical experiments to illustrate the behavior of our bounds as a function of k, τ , and ρ.

First we test whether our bound in Theorem 3.4 is descriptive by comparing it to empirical probabilities. For each value of k and τ tested, we compute 10 5 realizations of ψ k (B) and count the number of times that Results are given in Table 4.1, rounded to four digits. We write 1, as opposed to 1.0000, when the empirical probability is exactly 1, i.e., when all 10 5 realizations of ψ k (B) satisfy (4.1). We observe that if ρ is close to 1, the lower bound (3.9) can be pessimistic when k < 3 and/or τ < 3. However, when k ≥ 3 and τ ≥ 3, the bound gives a good indication of the actual probability. In case (a) the empirical probability is sometimes slightly lower than the lower bound. This is merely the effect of sampling error.

Conclusion.

The purpose of the present paper is to derive a bound of the form (1.1) for small-sample estimation. Our main result is the combined bound in Theorem 3.4. We derived this bound for the Gaussian estimator ψ k (B), but by Proposition 2.2 it also essentially holds for the GKL estimator η k (B) when k ≪ n, which is typically how the estimators are used in practice.

Unlike the original small-sample bounds from [START_REF] Gudmundsson | Small-sample statistical estimates for matrix norms[END_REF], but similarly to those in [START_REF] Avron | Randomized algorithms for estimating the trace of an implicit symmetric positive semi-definite matrix[END_REF][START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF], our bounds hold independently of the rank of B. Instead, they depend on ρ, the stable rank of B. Uniform bounds that hold for all B are readily obtained by replacing ρ by its lower bound 1 (or its upper bound min{m, n} in (3.11)).

Even when ρ is very close to 1, our bounds behave similarly to the bounds for the case rank(B) = 1 in (2.3), on which the small-sample estimation theory is based. Similarly to [START_REF] Roosta-Khorasani | Improved bounds on sample size for implicit matrix trace estimators[END_REF]Theorem 3], our bounds increase very quickly toward 1 as the stable rank of B increases. If B F is significantly larger than B 2 , say, B F ≥ 5 B 2 , then the probability of obtaining a good order-of-magnitude estimate is overwhelmingly close to 1. Thus, small-sample estimation works generally for matrices B not only of rank 1 but of any rank.

3 / 2 i- 1 (

 321 (4.1) B F τ ≤ ψ k (B) ≤ τ B F holds, divided by 10 5 . Without loss of generality we test matrices B that are square, diagonal, and of size n = 100. The diagonal elements σ i of B are chosen to obtain different values of the stable rank ρ, as follows: (a) ρ = 1 : σ 1 = 1 and σ i = 0 for i > 1. (b) ρ ≈ 1.2 : σ 1 = 10π/ √ 126 and σ i = 1/(i -1) for i > 1. (c) ρ ≈ 4 : σ i = max 10 -10 , √ d) ρ = 25 : σ 1 = 1 and σ i = 24/(n -1) for i > 1.

Table 4 . 1

 41 Empirical probability of (4.1) (top) and bound in (3.9) (bottom).

			(a) Case ρ = 1
	❍ ❍ ❍ ❍ τ k	1	2	3	4	5
	2	0.5717 0.0105	0.7608 0.4109	0.8551 0.6383	0.9071 0.7745	0.9393 0.8579
	3	0.7357 0.6035	0.8948 0.8765	0.9546 0.9512	0.9784 0.9783	0.9894 0.9899
	5	0.8410 0.8411	0.9600 0.9608	0.9893 0.9893	0.9970 0.9970	0.9988 0.9991
	10	0.9209 0.9203	0.9900 0.9900	0.9986 0.9986	0.9998 0.9998	0.9999 1.0000
			(b) Case ρ ≈ 1.2
	❍ ❍ ❍ ❍ τ k	1	2	3	4	5
	2	0.72360 0.03609	0.8849 0.4406	0.9455 0.6607	0.9729 0.7875	0.9869 0.8631
	3	0.9214 0.6249	0.9877 0.8664	0.9976 0.9389	0.9995 0.9697	0.9999 0.9845
	5	0.9955 0.8259	1 0.9528	1 0.9859	1 0.9956	1 0.9986
	10	1 0.9125	1 0.988	1 0.9982	1 0.9997	1 1.0000
			(c) Case ρ ≈ 4
	❍ ❍ ❍ ❍ τ k	1	2	3	4	5
	2	0.9878 0.2949	0.9995 0.6570	1 0.8125	1 0.8943	1 0.9399
	3	0.9998 0.5459	1 0.7941	1 0.9066	1 0.9576	1 0.9808
	5	1 0.6892	1 0.8521	1 0.9370	1 0.9749	1 0.9900
	10	1 0.8415	1 0.9608	1 0.9893	1 0.9970	1 0.9991
			(d) Case ρ = 25
	❍ ❍ ❍ ❍ τ k	1	2	3	4	5
	2	1 0.9703	1 0.9991	1 1.0000	1 1.0000	1 1.0000
	3	1 0.9928	1 0.9999	1 1.0000	1 1.0000	1 1.0000
	5	1 0.9968	1 1.0000	1 1.0000	1 1.0000	1 1.0000
	10	1 0.9978	1 1.0000	1 1.0000	1 1.0000	1 1.0000


Anthony Fillion
email: anthony.fillion@enpc.fr

Marc Bocquet

Serge Gratton

Quasi-static ensemble variational data assimilation: a theoretical and numerical study with the iterative ensemble Kalman smoother

The analysis in nonlinear variational data assimilation is the solution of a non-quadratic minimization. Thus, the analysis efficiency relies on its ability to locate a global minimum of the cost function. If this minimization uses a Gauss-Newton (GN) method, it is critical for the starting point to be in the attraction basin of a global minimum. Otherwise the method may converge to a local extremum, which degrades the analysis. With chaotic models, the number of local extrema often increases with the temporal extent of the data assimilation window, making the former condition harder to satisfy. This is unfortunate because the assimilation performance also increases with this temporal extent. However, a quasi-static (QS) minimization may overcome these local extrema. It accomplishes this by gradually injecting the observations in the cost function. This method was introduced by Pires et al. (1996) in a 4D-Var context.

We generalize this approach to four-dimensional strongconstraint nonlinear ensemble variational (EnVar) methods, which are based on both a nonlinear variational analysis and the propagation of dynamical error statistics via an ensemble. This forces one to consider the cost function minimizations in the broader context of cycled data assimilation algorithms. We adapt this QS approach to the iterative ensemble Kalman smoother (IEnKS), an exemplar of nonlinear deterministic four-dimensional EnVar methods. Using low-order models, we quantify the positive impact of the QS approach on the IEnKS, especially for long data assimilation windows. We also examine the computational cost of QS implementations and suggest cheaper algorithms.

1 Introduction

Context

Data assimilation (DA) aims at gathering knowledge about the state of a system from acquired observations. In the Bayesian framework, this knowledge is represented by the posterior probability density function (pdf) of the system state given the observations. A specificity of sequential DA is that observations are not directly available; they become available as time goes by. Thus, the posterior pdf should be regularly updated.

In order to do so, one usually proceeds in two steps: the analysis and the propagation (or forecast). During the analysis step, a background pdf is used as a prior together with the observation likelihood to construct the (often approximate) posterior pdf, following Bayes' theorem. During the propagation step, this posterior pdf is propagated in time with the model to yield the prior pdf of the next assimilation cycle.

In general these posterior and prior pdfs are not easily computable. In the Kalman filter, assumptions are notably made on the linearity of operators, to keep these pdfs Gaussian; this way, they are characterized by their mean and covariance matrix. Linear algebra is then sufficient to enforce both Bayes' theorem and the propagation step into operations on means and covariances.

However, with nonlinear models, the Kalman filter assumptions do not hold as the posterior and prior pdfs are not Gaussian anymore. A possibility in this case is to enforce Gaussianity with approximations. This requires the selection of mean and covariances intended for the Gaussian surrogate pdfs. With the Kullback-Leibler divergence, the best Gaus-

Ensemble variational methods

By contrast, four-dimensional ensemble variational (EnVar) schemes allow one to perform both a nonlinear variational analysis and a propagation of dynamical errors via the ensemble (see chap. 7 of [START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF]. The improvement brought by QS minimizations to these schemes has been suggested and numerically evaluated in Bocquet andSakov (2013, 2014) and [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF]. This motivates a more complete analytical and numerical investigation.

The iterative ensemble Kalman smoother (IEnKS) (Bocquet and Sakov, 2014;[START_REF] Asch | Data assimilation: methods, algorithms, and applications[END_REF] is the archetype of such a four-dimensional nonlinear EnVar scheme, where the ensemble parts of the algorithm are deterministic. Using low-order models (usually toy-models), it was shown to significantly outperform 4D-Var, the EnKF or the ensemble Kalman smoother in terms of accuracy.

The IEnKS improves the DA cycling by keeping track of the pdfs' mean and covariance matrix. To do this, a Laplace approximation is used to replace the posterior mean and covariance matrix with the minimizer of the cost function and an approximation of the inverse Hessian at the minimizer, respectively. These moments are then used to update the ensemble statistics. The updated ensemble is next propagated to estimate the prior mean and covariance matrix. Hence, it is also critical for the IEnKS to locate the global minimum of the cost function.

Here, we are interested in the application of the QS minimization to the IEnKS. One of the variant of the IEnKS called the multiple data assimilation (MDA) IEnKS was shown (Bocquet and Sakov, 2014) to be quasi-static by design and can be seen as the EnVar generalization of the sequential QSVA by [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF]. It was first tested in Bocquet and Sakov (2013). However, the MDA IEnKS is a specific variant of the IEnKS whereas we see here the IEnKS as an exemplar of deterministic nonlinear 4D EnVar methods. [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF] have applied QSVA numerically to a collection of hybrid and EnVar techniques on the Lorenz 1963 model [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF], where they vary the magnitude of nonlinearity. Nonetheless the focus of their study was not cycling and the transfer of information from one cycle to the next, which is critical to EnVar methods. They also showed that the ensemble transform Kalman smoother outperforms all of the benchmarked methods. We have furthermore shown that this smoother was systematically outperformed by the IEnKS by design, which was also demonstrated on numerics (Bocquet and Sakov, 2013). This strengthens our claim that the IEnKS can be used here as an exemplar of deterministic nonlinear 4D EnVar methods.

Outline

The rest of the paper is organized as follows. In Sect. 2, the performance dependency of 4D-Var and IEnKS algorithms sian approximation of a pdf is achieved by equating the mean and covariances (see, e.g., [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF]. However, the integrations necessary to evaluate these moments are also prohibitive.

In the 4D-Var algorithm (see, e.g., Lorenc, 2014, and references therein), the Laplace approximation gives us a way to work around the problem by replacing the posterior mean with the presumed unique minimizer of the cost function over all input values. A model propagation is then sufficient to estimate the prior pdf mean. This global approach calls for efficient global optimization routines. However, in practice, solving a global optimization problem is challenging when the number of unknowns is large, and local methods focused on finding a minimizer over an open subset like Gauss-Newton are often preferred (see, e.g., [START_REF] Björck | Numerical methods for least squares problems[END_REF].

Unfortunately, the Gauss-Newton method's ability to locate the global minimum depends on the minimization starting point and the cost function properties. Furthermore, missing this global minimum is likely to cause a quick divergence (from the truth) of the sequential DA method. Thus, it is critical for the assimilation algorithm to keep the minimization starting point in a global minimum basin of attraction.

Quasi-static variational data assimilation

Keeping the minimization starting point in a global minimum basin of attraction is constraining because, with a chaotic model, the number of local minima may increase exponentially with the data assimilation window (DAW) time extent L [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF][START_REF] Swanson | Four-dimensional variational assimilation and predictability in a quasi-geostrophic model[END_REF]. Unfortunately, assuming a perfect, chaotic -and hence unstable -model, the assimilation performs best for the longest time extents. Several strategies have been investigated to go beyond this restriction. [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] propose the quasi-static (QS) minimization in a 4D-Var context: as the observations are progressively added to the cost function, the starting point (or first guess) of the 4D-Var minimization is also gradually updated. This led to the method known as QSVA for quasistatic variational assimilation. [START_REF] Ye | Systematic variational method for statistical nonlinear state and parameter estimation[END_REF] propose the gradual increase of the model error covariances in the weakconstraint 4D-Var cost function in a minimization over an entire trajectory; this way the model nonlinearity is gradually introduced into the cost function (see also [START_REF] Judd | Gradient free descent: shadowing, and state estimation using limited derivative information[END_REF]. They also propose to parallelize this minimization over multiple starting points to increase the chance to locate the global minimum.

On the one hand, 4D-Var benefits from the QS approach to approximate the posterior and prior means. On the other hand, with traditional 4D-Var, the prior covariance matrix is taken as static. This is appropriate when only one cycle of assimilation is considered. But this limits the dynamical transfer of error statistics from one cycle to the next, for instance when [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] propose gradually moving their fixed-lag data assimilation window in order to build a sequential QSVA. on the DAW parameters is investigated. The emphasis is on the transfer of information from one cycle to next, which distinguishes the IEnKS from 4D-Var. In order to do this, a brief presentation of 4D-Var and the IEnKS algorithms is given. We then define a measure of performance for assimilation algorithms. This definition is used to give analytic expressions for the accuracy of both algorithms with a linear, diagonal, autonomous model. This quantifies the impact of cycling on the algorithms. After these preliminaries, the nonlinear, chaotic case is studied. In Sect. 3, we provide and describe the algorithms of the quasi-static IEnKS (IEnKS QS ). Section 4 is dedicated to numerical experiments with two Lorenz low-order models and to the improvement of the numerical efficiency. Conclusions are given in Sect. 5.

We emphasize that the algorithmic developments of this study are not meant to improve either high-dimensional or imperfect model data assimilation techniques. Even if [START_REF] Miller | Advanced Data Assimilation in Strongly Nonlinear Dynamical Systems[END_REF] show some similarities between the perfect and imperfect settings of a model of intermediate complexity, model error would generally forbid the use of very long DAWs as sometimes considered in this study. Instead, the objective of this paper is to better understand the interplay between chaotic dynamics, ensemble variational data assimilation schemes and their cycling, irrespective of whether they could be useful in high-dimensional systems.

The data assimilation window and assimilation performance

After reviewing 4D-Var (in a constant background matrix version) and the IEnKS algorithms, we investigate the dependency of assimilation performance on the DAW key parameters. This illustrates the cycling improvement brought in by the IEnKS compared to 4D-Var. We will see that, with chaotic models, the longer the DAW is, the better the accuracy of these algorithms. Which highlights the QSVA relevance for cycled data assimilation. The evolution and observation equations of the system are assumed to be in the following form:

y l = H (x l ) + ε l , (1a)
x l+1 = M (x l ) , (1b)
where the unknown state x l at time t l is propagated to t l+1 with the model resolvent M : R m → R m . The model is assumed to be perfect (no errors in Eq. ( 1b)) and autonomous (M does not depend on time). The observation operator H : R m → R d relates the state x l to the observation vector y l . The observation errors (ε l ) l ≥ 0 are assumed to be Gaussian with mean 0 ∈ R d and covariance matrix R ∈ R d×d ; they are uncorrelated in time.

t 0 x 0 t K y K t L y L L S observations
Figure 1. Schematic of a DAW. The state variable at t 0 is x 0 , the observation vector y K at time t K is the first of the DAW to be assimilated and the observation y L at present time t L is the last to be assimilated. These observations, possibly observation vectors, are represented by black dots.

4D-Var and IEnKS algorithms

Both 4D-Var and the IEnKS use a variational minimization in their analysis step. The objective of this minimization is to locate the global maximum of the posterior pdf p x 0 |y L:K of the system past state x 0 given the observations y L:K = y K , . . ., y L at times t L:K = [t K , . . ., t L ] ∈ R S .

The system state and observations are seen as random vectors with values in R m and R d , respectively. The posterior pdf quantifies how our knowledge on the state x 0 changes with realizations of y L:K . Thus, its maximum is the most probable state after assimilating the observations. The DAW is displayed in Fig. 1. The parameters K and L are the time index of the DAW first and last assimilated observation batch, respectively. The number of observation vectors used within the DAW is L -K + 1.

To specify this posterior pdf, we have to make further assumptions on x 0 .

The initial state x 0 is assumed to be Gaussian with mean

x b 0 ∈ R m and covariance matrix B ∈ R m×m : p (x 0 ) = N x 0 |x b 0 , B . (2)
Following these assumptions, an analytic expression can be obtained for the posterior smoothing pdf p x 0 |y L:K at the first cycle, or for the cost function associated with this pdf. The latter is defined as

G x 0 |y L:K = -ln p x 0 |y L:K . (3)
The notation G is used rather than the traditional J to refer to an exact cost function, i.e., a cost function defined from an exact posterior pdf. Bayes' theorem consequently yields the following:

G x 0 |y L:K = -ln p (x 0 ) -ln p y L:K |x 0 + ln p y L:K , (4)
and with the Gaussian assumption applied to the background and observation errors we find

G x 0 |y L:K = 1 2 x b 0 -x 0 2 B -1 + c 0 + 1 2 L l=K y l -H • M l (x 0 ) 2 R -1 , (5)
where x 2 A = x T Ax is the norm of x associated with a symmetric positive definite matrix A; c 0 is a normalization constant ensuring e -G(x 0 |y L:K ) dx 0 = 1; M l stands for l compositions of M.

The propagation corresponds to a time shift of S time steps. Thus, at the kth assimilation cycle, the posterior pdf is p x kS |y kS+L:K . Using Bayes' theorem the kth cycle cost function is

G x kS |y kS+L:K = -ln p x kS |y (k-1)S+L:K + 1 2 L l=K y kS+l -H • M l (x kS ) 2 R -1 + c kS , (6)
whereln p x kS |y (k-1)S+L:K is the background term. If the model and observation operators are nonlinear, an analytical expression for this prior is not accessible and one needs to approximate it. The 4D-Var and the IEnKS algorithms are solutions based on distinct approximation strategies. The 4D-Var cost function at the kth cycle, based on the static error covariance matrix B, is defined by

J x kS ; y kS+L:kS+K , x b kS = 1 2 x b kS -x kS 2 B -1 + 1 2 L l=K y kS+l -H • M l (x kS ) 2 R -1 . (7)
The analysis of 4D-Var consists in by minimizing Eq. ( 7), yielding x a kS at t kS . Because this cost function depends on realizations of the random observations, x a kS is also a random variable. This analysis is then propagated at time t (k+1)S with the resolvent of the model to produce the next cycle background state:

x b (k+1)S = M S x a kS . (8)
In general, G x 0 |y L:K and J x 0 ; y L:K , x b 0 only coincide at the first cycle of an assimilation because of the assumption outlined by Eq. (2). Subsequently, the background term of the 4D-Var cost function 1

2 x b kS -x kS 2 B -1
is a Gaussian approximation of the exact background term ln p x kS |y (k 1)S L K . By definition, the background er-B

where observations are assimilated only once, and the multiple data assimilation (MDA) version where they are assimilated several times. We focus on the SDA version in the theoretical development. The MDA version, which can be seen as the first published quasi-static EnVar scheme is used in the numerical experiments for comparison. Note that, for the SDA IEnKS, the number of observations is

L -K + 1 = S.
At the kth cycle of the IEnKS, the background ensemble at t kS is obtained by a propagation from the previous cycle. The ensemble members are the columns of the matrix E b kS , which is seen as a random matrix with values in R m×n . It is used to estimate the prior mean and covariance matrix as follows:

E x kS |y (k-1)S+L:K ≃ x b kS , (9)
C x kS |y (k-1)S+L:K ≃ X b kS X bT kS , (10)
where E and C are the expectation and covariance operators, respectively; x b kS and X b kS are the empirical mean and normalized anomaly of E b kS , respectively:

x b kS = E b kS 1 n n , (11)
X b kS = E b kS I n - 1 n 1 T n n √ n -1 , (12)
with 1 n = [1, . . ., 1] T ∈ R n a vector of ones and I n is the identity of R n . Note that Eqs. ( 9) and ( 10) are approximations because of sampling errors. If the state vector is of the form

x kS = x b kS + X b kS w kS , (13)
with w kS ∈ R n the control variable in the ensemble space, the IEnKS cost function is defined in the ensemble space by

J w kS ; y kS+L:kS+K , E b kS = 1 2 w kS 2 + 1 2 L l=K y kS+l -H • M l x b kS + X b kS w kS 2 R -1 , (14)
where • = • I with I the identity matrix. The analysis of the IEnKS consists in minimizing this cost function. It yields an updated ensemble E a kS at time t kS verifying

x a kS = x b kS + X b kS w a kS , (15)
X a kS = X b kS ∇ 2 J w a kS ; y kS+L:kS+K , E b kS -1/2 U, (16)
with x a kS and X a kS representing the empirical mean and normalized anomalies of the analyzed ensemble, respectively, w a kS the cost function minimizer, ∇ 2 J w a kS ; y kS+L:kS+K , E b kS its Hessian at this minimum usually approximated in the Gauss-Newton method by

∇ 2 J w a kS ; y kS+L:kS+K , E b kS ≃ I n + L l=K F T l R -1 F l , (17)
ror covariance

- matrix + :
of the traditional 4D-Var cost function is the same for each cycle -this is not the case for the IEnKS.

The IEnKS (Bocquet and Sakov, 2014) is an ensemble method with a variational analysis. Two versions of the algorithm exist: the singular data assimilation (SDA) version Analysis Propagation Analysis

x b 0 , E b 0 x a 0 , E a 0 y K y L x b S , E b S x a S , E a S y K+S y L+S M S : 0 Figure 2.
Chaining of the 4D-Var and IEnKS first two cycles with S = 3, L = 4. The first 4D-Var analysis uses the background x b 0 at t 0 and the observations y L:K to give the analysis x a 0 at t 0 . It is propagated S steps forward in time to produce the new background at time t S where another analysis can be performed. The IEnKS does the same but with an ensemble. The dashed rectangle symbolizes the current DAW, black dots represent the observations assimilated in the current cycle, gray dots represent already assimilated observations and white dots represent observations not assimilated.

with

F l = dF l dw kS w a kS ,
where

F l : w kS -→ H • M l x b kS + X b
kS w kS , in order to avoid computing the model second derivatives. The exponent -1/2 refers to the unique symmetric definite positive inverse square root of a symmetric definite positive matrix and U ∈ R n×n an orthogonal matrix such that U1 n = 1 n . This analyzed ensemble is also propagated to time t (k+1)S to produce the next cycle background ensemble

E b (k+1)S = M S E a kS . ( 18
)
The cycle is completed by using the cost function Eq. ( 14) with time indexes incremented by S in the next analysis. 4D-Var and the IEnKS are displayed in Fig. 2.

Performance of assimilation

In order to evaluate the efficiency of 4D-Var and the IEnKS with DAW parameters S and L, two measures of accuracy are investigated here: the usual empirical RMSE, and a theoretical counterpart.

At the kth cycle, the algorithm generates at time t kS an analysis x a kS from the observations. This analysis is propagated with the model l steps forward in time to yield the analysis x a kS+l meant to approximate the system's true state x kS+l . A traditional measure of assimilation performance is the root mean square error (RMSE). It is defined by

RMSE = 1 √ m x kS+l -x a kS+l . (19)
The RMSE takes different names depending on when it is computed. If l = L, it is called the filtering RMSE; if 0 ≤ l ≤ L -1, it is the smoothing RMSE with lag Ll. In the following, the smoothing RMSE will correspond to the one with (maximum) lag L.

The RMSE rigorously depends on the random variable realizations, and thus it is also a random variable. In our numerical experiments, as is usually the case, the RMSE is averaged over the cycles to mitigate this variability:

aRMSE N = 1 N N -1 k = 0 1 √ m x kS+l -x a kS+l . (20)
Let us assume that there is a random couple x ∞S+l , x a ∞S+l whose distribution is invariant and ergodic with respect to the shift transformation:

T : x kS+l , x a kS+l -→ x (k+1)S+l , x a (k+1)S+l . (21)
Then by Birkhoff's ergodic theorem (see [START_REF] Walters | An Introduction to Ergodic Theory[END_REF] the sequence (aRMSE N ) N converges when N → ∞ and its limit aRMSE verifies

aRMSE = 1 √ m E x ∞S+l -x a ∞S+l , (22)
where the expectation E is taken over p x ∞S+l , x a ∞S+l . In this case, the aRMSE measures the long-term impact of the cycling on the assimilation accuracy. This limit is difficult to exploit algebraically. That is why, in the theoretical developments, we prefer the expected MSE (eMSE), denoted by P as follows:

P kS+l = E x kS+l -x a kS+l 2 ,
where the expectation is taken over p x kS+l , x a kS+l . In the following subsection, we will focus on the long term impact of the cycling on P . Simplifying assumptions will be made to express P ∞S+l = lim k→∞ P kS+l as a function of S and L.

Performance in the linear, diagonal, autonomous case

In order to obtain analytical expressions of the eMSE for 4D-Var and the IEnKS, we make drastic simplifying assumptions.

First, the model is assumed to be the resolvent of a linear, diagonal, autonomous ordinary differential equation. Thus, it can be expressed as

M l (x + δx) = M l (x) + M l δx, (23)
where M = diag(α i ) i=1..m is diagonal and does not depend on x. We further assume H = hI m , B = bI m , R = rI m , where I m is the identity matrix of R m and h, r, b > 0. With these assumptions, Appendix A provides an expression for the 4D-Var asymptotic eMSE in the univariate case. The generalization to the diagonal multivariate case is obtained by summing www.nonlin-processes-geophys.net/25/315/2018/ Nonlin. Processes Geophys., 25,[315][316][317][318][319][320][321][322][323][324][325][326][327][328][329][330][331][332][333][334]2018 up the eMSEs of each direction:

P 4D-Var ∞S+l = m i=1      ∞ if i ≥ 1 b 2 L K,i α 2(S-l) i i 1 -i otherwise , (24a) L K,i = h 2 r α 2(L+1) i -α 2K i α 2 i -1 , ( 24b
) i = α 2S i 1 + b L K,i 2 . ( 24c
)
The case i ≥ 1 means that too much credit is given to the background variance, which is approximated for all cycles by the constant b in our 4D-Var scheme. Therefore, the information carried by the observations is not sufficient to mitigate the exponential growth of errors in the propagation.

Concerning the IEnKS, the anomalies are assumed to be full rank to avoid any complication due to singular covariance matrices. Moreover, the linearity of the model is employed to express the background statistics:

x b (k+1)S = M S x a kS , (25a)
X b (k+1)S = M S X a kS , (25b)
x b 0 = x b 0 , (25c)
X b 0 = B 1/2 . ( 25d
)
This way, sampling errors are avoided as they are not the focus of this study. The background ensemble E b kS becomes a notational shortcut for the pair x b kS , X b kS . This simplified IEnKS is actually a Kalman smoother [START_REF] Cosme | Smoothing Problems in a Bayesian Framework and Their Linear Gaussian Solutions[END_REF][START_REF] Bocquet | Four-dimensional ensemble variational data assimilation and the unstable subspace[END_REF]. With these assumptions, Appendix B gives an expression for the IEnKS asymptotic eMSE in the univariate case. The optimality of the IEnKS eMSE is also proven. The generalization to the diagonal multivariate case is also obtained by summing up the eMSEs of each direction:

P IEnKS ∞S+l = m i=1    0 if |α i | ≤ 1 r h 2 α 2(L-l) i α 2 i -1 α 2 i otherwise . ( 26
)
In the following, we study the eMSE dependency on the DAW parameters, S and L, for both algorithms. We focus on a bivariate case with α 1 = 1.2, α 2 = 0.8 in order to have one stable and one unstable direction; h = b = r = 1. Using Eqs. ( 24) and ( 26), the asymptotic smoothing and filtering eMSEs are displayed as a function of L, S (Fig. 3). The eMSE components on the stable and unstable directions are also shown. Those graphs are interpreted in the following.

Specifically, the eMSE expression for the IEnKS is of the form

P IEnKS ∞S+l = c 1 α 2(l-L) 1 , (27)
where c 1 does not depend on S, L, l. The contribution to P IEnKS ∞S+l on the stable direction is zero. It depends only on the lag Ll; it does not depend on S. Moreover, the filtering eMSE (l = L) is constant: the propagation compensates for the analysis. The smoothing eMSE (l = 0) decreases exponentially with L.

Concerning 4D-Var, we assume K to be fixed and S → ∞ to give an asymptotic eMSE expression:

P 4D-Var ∞S+l = P IEnKS ∞S+l (1 + o (1)) + c 2 α 2l 2 (1 + o (1)) , (28)
where c 2 is constant with S, l and o (1) → 0 when S → ∞.

The unstable component is close to the IEnKS overall eMSE.

The biggest difference with it concerns the eMSE on the stable component. The inexact background variance modeling adds a detrimental term to the eMSE. To qualify the long term impact of the cycling on the errors, the filtering eMSE is more instructive than the smoothing eMSE. Indeed, the smoothing eMSE is improved with L as it adds future observations (with respect to the analysis time) in the DAW. This improvement dominates the potentially detrimental impact of the Gaussian background approximations. In the filtering eMSE, this improvement is balanced by the propagation of the analysis at the end of the DAW. In Fig. 3, the filtering eMSE stable component is mitigated such that it has little effect. However, the parameter S improves the filtering eMSE on the unstable component. The bigger S becomes, the closer P 4D-Var ∞S+l is to P IEnKS ∞S+l , which is optimal (cf. Appendix B). A qualitative explanation is that to assimilate the same numbers of observations, a 4D-Var using high values of S needs fewer cycles; therefore, it less often relies on the background approximation, making the analysis more trustworthy.

Figure 4 displays the asymptotic eMSEs of both algorithms as a function of the lag for S = L = 5; these curves are similar to those of [START_REF] Trevisan | Four-dimensional variational assimilation in the unstable subspace and the optimal subspace dimension[END_REF]. Concerning 4D-Var, the eMSE can be written as a function of the lag Ll as follows:

P 4D-Var ∞S+l ≃ c 1 α -2lag 1 + c 2 α 2(L-lag) 2 . ( 29
)
Thus, the unstable component of the eMSE is an exponentially decreasing function of the lag and the stable component

This expression shows that the eMSE components on the stable directions are null. Indeed one expects the IEnKS to be at least more efficient than a free-run, whose errors in the stable directions tend to zero. 1 This is not the case for 4D-Var since the static background covariance matrix introduces spurious variance in the stable directions as seen in Eq. ( 24).

In [START_REF] Trevisan | Four-dimensional variational assimilation in the unstable subspace and the optimal subspace dimension[END_REF], 4D-Var error variances in the stable directions are forced to zero to improve the accuracy of the assimilation.

1 The fact that the errors lie in the unstable subspace is more general. [START_REF] Bocquet | Four-dimensional ensemble variational data assimilation and the unstable subspace[END_REF]) balances the gain in eMSE due to the assimilation of future observations. Thus, the filtering eMSE is not improved by the assimilation of observations distant in time; it is affected by the background pdf approximation. To assimilate the same number of observations, algorithms using high values of S need fewer cycles. Therefore, they less often rely on the background approximation. That is why the 4D-Var filtering performance is improved with S. For the IEnKS as in Sect. 2.3, the prior pdf approximation is exact so that there is no dependence on S. This is no longer true in the nonlinear case. Hence, with a similar reasoning, one actually expects an improvement on the IEnKS filtering performance with S.

As a matter of fact, it has been shown that with a nonlinear chaotic model, the filtering accuracy increases with L in most cases (see Bocquet and Sakov, 2014, and Sect. 4 of the present paper).

is an exponentially increasing function of the lag. The sum therefore decreases when the unstable component is dominant and increases when the stable component is dominant.

In this section, we studied the accuracy of both the cycled IEnKS and 4D-Var eMSEs as a function of the DAW parameters under linear, autonomous, diagonal assumptions. We found that the DAW parameter L improves the smoothing eMSE and S improves the filtering eMSE. These properties will be discussed and numerically investigated in a nonlinear context in the next section.

Performance in the nonlinear, chaotic case

The results of Sect. 2.3 foster the use of the largest possible L to improve the 4D-Var and IEnKS smoothing eMSEs. For filtering, the error propagation at the end of the DAW

Multiple local minima

However, with a chaotic model, [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] show that the 4D-Var cost function number of local extrema increases with L, making minimization problematic. We show in this section that the IEnKS cost function suffers from the same problem.

This behavior will be illustrated with the Lorenz 95 (L95) model [START_REF] Lorenz | Optimal Sites for Supplementary Weather Observations: Simulation with a Small Model[END_REF]. It represents a midlatitude zonal circle of the atmosphere and is described by a set of m nonlinear differential equations:

dx j dt = x j +1 -x j -2 x j -1 -x j + F, (30)
where x j is the j -th modulo m component of x, m = 40 and F = 8. This equation is integrated using a fourth-order Runge-Kutta scheme with a time step of δt = 0.05. The dynamics of L95 are chaotic; the L95 model's largest Lyapunov exponent is λ ≃ 1.7.

Figure 5 shows a typical IEnKS cost function profile in one verge towards any local extremum, take many iterations or even diverge. However, if the cost function is quadratic, the global minimum is reached in one iteration. The cost function non-quadraticity is induced by the model nonlinearity.

In the following, we will give a heuristic argument yielding a bound on the S parameter beyond which the GN method probably misses the global minimum in the configuration

H = B = R = I m .
To some extent, our argument can be seen as an improvement on the notion of useful DAW length introduced by [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] beyond which the performance gain is negligible. In contrast with the useful DAW length, we account in the following for both the cycling and the nonlinearity.

Effective data assimilation window length

First, the GN convergence properties are drastically simplified. We assume the method converges to the global minimum if, and only if, the minimization starting point is in a neighborhood of the global minimizer where the IEnKS cost function is almost quadratic.

Unfortunately, this minimizer is unknown because the cost function depends on realizations of many random variables. In order to eliminate this variability, [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] introduced a so-called error-free cost function. We prefer an averaged cost function J ∞S in this particular study defined by

J ∞S (w) = lim N→∞ 1 N N -1 k = 0 J w; y kS+L:kS+K , E b kS . (31)
Relying on an ergodicity assumption, Appendix C proves that this averaged cost function verifies

J ∞S (w) = 1 2 w 2 + dS 2 + 1 2 L l=K E δx b ∞S+l 2 , ( 32a
)
δx b ∞S+l =M l x b ∞S + X b ∞S w -M l (x ∞S ) , (32b)
where the ergodic random variables x ∞S , X b ∞S and x ∞S have been defined in Sect. 2.2 and Appendix C. As seen in Eq. ( 32), a sufficient condition for the starting point w = 0 to be in a neighborhood of the global minimizer where the cost function is assumed almost quadratic requires that x b ∞S be in a neighborhood of x ∞S where all the M l K≤l≤L are almost linear.

In the univariate case, if the model behavior is almost linear and unstable, we can use Eq. ( 26) to estimate the terms in the sum in Eq. ( 32) at the starting point w = 0:

E δx b ∞S+l 2 ≃ α 2(S+l-L) α 2 -1 α 2 , ( 33
)
where α is the linear part of the model and the extra S accounts for the propagation. But in a necessarily bounded direction of the analyzed ensemble space for multiple values of the DAW parameters. The system is observed at every time step and

H = B = R = I m .
The curves have more and more local extrema as L increases. The curves with the highest amplitudes of the ripples are found for small values of S. Indeed, an averaging effect may settle in as the number of observations increases.

This hilly shape causes minimization problems. A possible minimization procedure for the IEnKS is the Gauss-Newton (GN) algorithm (e.g., [START_REF] Björck | Numerical methods for least squares problems[END_REF]. GN is not a global procedure meaning that, depending on the starting point and the cost function properties, the algorithm can con- physical system, the right-hand side of Eq. ( 33) cannot grow indefinitely with l + S. Such model saturation imposes

E δx b ∞S+l 2 ≤ B, (34)
with B representing a bound. Hence, Eqs. ( 33) and ( 34) yield the following inequality on S:

S ≤ S max , (35a)
S max = ln (B) -ln 1 -α -2 2 ln (α) . ( 35b
)
We choose l = L because it corresponds to the most constraining case. When Eq. ( 35a) is violated, x b ∞S departs from x ∞S such that the nonlinearities of M L are significant.

To apply this inequality to the L95 model we choose the following:

α = lim N→∞ 1 N N-1 k = 0 σ dM dx (x k ) , (36)
σ being the mean of the singular values greater than 1. This corresponds to an average of the error amplification by dM dx (x ∞S ) in the local unstable subspace. We also choose the average squared norm between two long trajectories for the bound B:

B = lim N→∞ 1 mN N -1 k = 0 x kS -M kS x b 0 2 . ( 37
)
This quantity is greater than E δx b ∞S+l 2 because the IEnKS asymptotic performance is at least better than a free run. From the values of α and B we find S max = 14.

Figure 6 shows the filtering and smoothing aRMSEs of an IEnKS L = S with L95 as a function of S; the performance strongly deteriorates for S > 16, which is remarkably consistent with our estimation.

Figure 6 also shows another difference with the linear case: the IEnKS filtering aRMSE depends on S. The former discussion on the local extrema explains this dependency for large values of S.

However, for small values of S, the decreasing aRMSE has not been explained. This is a consequence of the Gaussian background approximation. At each cycle, the IEnKS uses the background ensemble E b kS to estimate the first two moments of the background pdf and makes the following ap- Gauss-Newton IEnKS L = S as a function of S with the L95 model over 5 × 10 5 cycles (logarithmic scale). We used the finite-size version (Bocquet, 2011) to account for sampling errors and avoid the need for inflation, H = R = B = I m , n = 20. proximation:

p x kS |y (k-1)S+L:K ≃ N x kS |x b kS , X b kS X b kS T . ( 38
)
This procedure is directly applicable to the IEnKS cost function minimization. The left panel in Fig. 7 is a schematic of a QS minimization and Algorithm 1 gives the pseudo-code of an IEnKS with a QS minimization. The new parameters (L q ) q<N Q control the number of observations added at each minimization, where N Q is the total number of batches of observation vectors.

The first three lines initialize the minimization starting point, the ensemble mean and anomaly matrix. The for loop in lines 4-23 repeats the QS minimization. The while loop in lines 6-22 is the Gauss-Newton minimization. Lines 7 and 8 center the ensemble on the current minimizer. Lines 9 and 10 initialize the cost function gradient and the approximate Hessian. The for loop in lines 12-18 compute the observation terms of the cost function, the gradient and the approximate Hessian. Lines 13 and 14 use a finite difference formula to compute the tangent linear and adjoint of w -→ H • M l x b 0 + X b 0 w . Lines 15 and 16 use this adjoint to update the gradient and approximate the Hessian. Lines 19 and 20 solve the linear system of the Gauss-Newton algorithm to update the current minimizer. When GN convergence is reached, this minimizer will be used as a starting point for the next QS minimization. Line 24 updates the ensemble. Line 25 propagates the updated ensemble to the next assimilation cycle. Figure 5 illustrates the QS scheme on a single analysis, as described in the caption.

To cycle the scheme, the DAW is then shifted with a small S. This ensures minimal cost function deformation, since few vectors of observation enter and few leave the DAW. This new cost function is then minimized using the forecast of the preceding minimizer as a starting point.

To keep this cost function deformation statistically consistent, Bocquet and Sakov (2014) advocate the use of S = 1 and assimilate one observation vector at the end of the DAW, which avoids multiple assimilation of the same observations. It is also the easiest way to ensure inequality (Eq. 35a). However, this method has been shown to be suboptimal because of the frequent Gaussian background approximations. Moreover, it ultimately fails to be QS since there is only one distant observation vector per analysis.

An alternative is to keep S = 1 but use all observations in the DAW and, consequently, relax the condition K = L-S + 1. This way, observations are assimilated several times. This is done in the MDA IEnKS (Bocquet and Sakov, 2014). To keep the statistics consistent, at least in the linear/Gaussian case, the observations error covariances should be adequately altered. Hence, the MDA IEnKS is truly a QS scheme and makes a good reference scheme for our numerical experiments. However, the multiple assimilation of the observations introduces spurious correlations in the nonlinear/non-Gaussian case, which entail sub-optimality. A scheme similar to the IEnKS QS has also been successfully used in [START_REF] Carrassi | Estimating model evidence using data assimilation[END_REF] to compute model evidence. Indeed, the efficient computation of model evidence as an integral over the state space depended on the proper identification of a global Because the model is nonlinear, this pdf is unlikely to be Gaussian. Therefore, this approximation results in a loss of information and the more it is used, the farther from G the IEnKS cost function is. This is exactly what happens when S is small: to assimilate the same number of observations, the IEnKS uses more cycles so that it relies more on the Gaussian background approximation.

Quasi-static algorithms

We have seen in the previous section that the effective DAW length is constrained by the cost function non-quadraticity. In this section we review and propose algorithms able to overcome these minimization issues and reach longer DAWs.

Quasi-static algorithms were introduced by [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] in a 4D-Var context. The idea behind QSVA is to control the way observations are added to the cost function in order to keep the minimization starting point in the basin of attraction containing the global minimizer. The method is carried out by repeating the minimization with an increasing number of observations: the first minimization is performed using the cost function with a single observation vector, then the number of observation vectors is increased and another minimization can be performed with the former minimizer as the new starting point. The process is then repeated until all observations are accounted for.

t 0 t L0 y K y L

Gauss Newton

Starting point:

w a(0) 0
Minimizer:

w a(1) 0 t 0 t L1 y K y L

Gauss Newton

Starting point:

w a(1) 0
Minimizer:

w a(2) 0 t 0 t L2 y K y L

Gauss Newton

Starting point:

w a(2)
Minimizer:

w a(3) N Q steps t 0 y K t L0 y L j 0 GN loops
Starting point:

w a(0) 0
Minimizer:

w a(1) 0 t 0 t L1 y K y L j 1 GN loops
Starting point:

w a(1) 0
Minimizer:

w a(2) 0 t 0 t L2 y K y L j 2 GN loops
Starting point:

w a(2)
Minimizer:

w a(3) N Q steps Figure 7
. Schematics of an IEnKS QS minimization (a) and an IEnKS QC minimization (b). The rectangle contains the observations to be assimilated and the black dots represent the observations used in the current minimization. The "Gauss-Newton" surrounded by arrows represents the iterations of the Gauss-Newton procedure. The number of quasi-static steps is N Q = 3. The flow of observations is controlled by the parameters (L 0 , L 1 , L 2 ) = (3, 5, 6). For the QC IEnKS, the number of GN iterations is controlled by the (j 0 , j 1 , j 2 ) parameters.

maximum of the integrand. However, its implementation was based on the update of the ensemble whenever an observation batch is added to the cost function, which is not as numerically efficient as the scheme presented here.

The success of the QS minimization lies in the fact that, when an observation is successfully assimilated, the eMSE is reduced. Thus, the analysis probability mass concentrates around the true state. The analysis is then more likely to be in a neighborhood of the true state where the model is linear. The cost function non-quadraticity can then be increased by adding a new term in it. This is confirmed by the following argument: let P (q) be the IEnKS asymptotic eMSE at the qth step of a QS minimization. With the notations and assumptions of Sect. 2.3, i.e., in a linear context, we have the following recurrence relation:

P (q+1) -1 = P (q) -1 + L q+1 L q +1 , (39a)
P (-1) = α 2S α 2 -1 α 2 L N Q -1 +1 . ( 39b
)
propagation of errors does not exceed the bound:

α 2L q+1 P (q) ≤ B. (40)
This implies

L q+1 ≤ L q + S max , (41a)
L 0 ≤ L N Q -1 + S max -S, (41b)
which yields S ≤ N Q S max . Therefore, the QS minimizations allow for a N Q times longer DAW.

Unfortunately, these QS minimizations are very expensive. Indeed, they add a third outer loop repeating N Q GN minimizations. The GN iterations used to compute the intermediate starting points give unnecessary precision; all that is required for these starting points is to be in a neighborhood of x ∞S where the model is almost linear. Thus, one can restrain the number of intermediate GN loops and save the full convergence to the last minimization. This is done in the quasi-convergent IEnKS (IEnKS QC ) with the parameters (j q ) q<N Q . These parameters correspond to the numbers of GN loops in the intermediate QS minimizations. They are typically equal to 1 except for the last one. Algorithm 2 gives the pseudo-code of the IEnKS QC and the right panel in Fig. 7 is a schematic for it. Thus, P (q+1) < P (q) and we can increase the cost function non-quadraticity by adding new terms in it as long as the Algorithm 1 One cycle of the IEnKS QS Require: E b 0 the background ensemble at t 0 ; λ the inflation; (L q ) q<NQ a list of DAW time indexes; ǫ the finite differences step; δ, j max GN end of loop parameters; 1 = (1, . . . , 1)

T ∈ R n and I is the identity of R n .

1: w a 0 = 0 2: xb 0 = E b 0 1/n 3: X b 0 = λ E b 0 -xb 0 1 T 4: for q = 0...N Q -1 do 5: j = 0 6: repeat 7: xa 0 = xb 0 + X b 0 w a 0 8: E 0 = xa 0 1 T + ǫX b 0 9: ∇J = (n -1)w a 0 10: ∇2 J = (n -1)I
11:

E K = M K (E 0 ) 12:
for l = K, ..., L q do 13:

ȳl = H (E l ) 1/n 14: Y l = H (E l ) -ȳl 1 T /ǫ 15: ∇J = ∇J -Y T l R -1 (y l -ȳl ) 16: ∇2 J = ∇2 J -Y T l R -1 Y l 17: E l+1 = M (E l ) 18:
end for 19:

solve ∇2 Jδw = ∇J 20:

w a 0 = w a 0 -δw 21: j = j + 1 22: until ||δw|| ≤ δ or j ≥ j max 23: end for 24: E a 0 = x a 0 1 T + √ n -1X b 0 ∇2 J -1/2 25: E b S = M S (E a 0 )
4 Numerical experiments with low-order models

In the following, we perform numerical experiments with the Lorenz 1963 (L63) and Lorenz 1995 (L95) models. L95 has already been presented in Sect. 2.4. L63 [START_REF] Lorenz | Deterministic Nonperiodic Flow[END_REF]) is a simplified model for atmospheric convection. It is defined by the following ordinary differential equations:

dx dt = σ (y -x) , (42a)
dy dt = ρx -y -xz, (42b)
dz dt = xy -βz. (42c)
The IEnKS parameters are ǫ = 10 -4 , δ = 10 -3 , j max = 20, N Q = 1, L 0 = L. For the IEnKS QS , the QS parameters are L 0 , L 1 , . . ., L N Q -1 = K + 0, 1, . . ., N Q -1 × S-1 N Q -1 . For the IEnKS QC , the QS parameters are the same and, in addition, j 0 , . . ., j N Q -2 , j N Q -1 = (1, . . ., 1, 20). Sampling errors are systematically accounted for using the IEnKS finitesize version (Bocquet, 2011;[START_REF] Bocquet | Combining inflation-free and iterative ensemble Kalman filters for strongly nonlinear systems[END_REF]Bocquet et al., 2015) which avoids the need for inflation and its costly tuning. Finally the aRMSE is averaged over a number of cycles which is determined by the number of observations assimilated. We use 5 × 10 5 observation vectors for L95 and 5 × 10 6 observation vectors for L63.

Unlike [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF], our numerical experiments neither address increasing nonlinearity, nor do they address the use of climatological background error covariance matrices. Instead, we focus exclusively on the IEnKS performance dependence on the DAW key parameters L, S and the number N Q of QS minimizations. [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF] numerically evaluates the QSVA approach with hybrid and EnVar techniques, similarly to Bocquet and Sakov (2013), and confirms the findings of [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] albeit in a En-Var context. In terms of accuracy, [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF] show that the ensemble transform Kalman smoother (ETKS) outperforms all hybrid schemes in their numerical experiment. Since the IEnKS systematically outperforms the ETKS in all These equations are integrated using a fourth-order Runge-Kutta scheme with a time step of δt = 0.01 and (σ,ρ,β) = (10,28,8/3). The dynamics of the L63 model are chaotic, with a largest Lyapunov exponent given by λ ≃ 0.91.

Both models are assumed perfect. The truth run is generated from a random state space point. The initial ensemble is generated from the truth with B = I m where m = 40,3 and n = 20,3 for L95 and L63, respectively. Observation vectors are generated from the truth with H = R = I m every t = 0.05 for L95 and every t = 0.02 for L63. A burn-in period of 5 × 10 3 × t is enforced in both cases.

Algorithm 2 One cycle of the IEnKS QC Require: E b 0 the background ensemble at t 0 ; λ the inflation; (L q ) q<NQ a list of DAW time indexes; (j q ) q<NQ the number of intermediate GN loops; ǫ the finite differences step; δ GN end of loop parameter; 1 = (1, . . . , 1)

T ∈ R n and I is the identity of R n ; 1: w a 0 = 0 ∇2 J = (n -1)I

2: xb 0 = E b 0 1/n 3: X b 0 = λ E b 0 -xb 0 1 T 4: for q = 0...N Q -1 do
11:

E K = M K (E 0 )
12:

for l = K, ..., L q do 13:

ȳl = H (E l ) 1/n 14: Y l = H (E l ) -ȳl 1 T /ǫ 15: ∇J = ∇J -Y T l R -1 (y l -ȳl ) 16: ∇2 J = ∇2 J -Y T l R -1 Y l 17: E l+1 = M (E l ) 18:
end for 19:

solve ∇2 Jδw = ∇J 20:

w a 0 = w a 0 -δw 21:

j = j + 1 22:
until ||δw|| ≤ δ or j ≥ j q 23: end for 24:

E a 0 = x a 0 1 T + √ n -1X b 0 ∇2 J -1/2 25: E b S = M S (E a 0 )
lated and condensed in the background approximation. Thus, the performance gain with the DAW length comes from the precision of this Gaussian background approximation; a precision that the linearized theory is not able to provide, attested to by the filtering (l = L) performance independence of Eq. ( 26) with the DAW parameters.

Figure 10 compares the smoothing aRMSE (first column), the filtering aRMSE (second column) and the number of ensemble propagations (third column) of the IEnKS QS (N Q = S, S = L), the IEnKS (S = L) and the IEnKS-MDA (S = 1), for both L95 and L63 as a function of L. The number of ensemble propagations is the total number of ensemble propagations in units of t divided by the total number of assimilated observation vectors. For L < 20, all three algorithms show smoothing and filtering performance improvements with L. For L > 20, the IEnKS filtering and smoothing RMSEs increase because of the multiple local extrema. For L < 40, the IEnKS QS has smaller aRMSEs than the IEnKS-MDA. Because the IEnKS QS is SDA by design, it does not suffer from sub-optimality related to multiple assimilations and nonlinearity. Moreover, the IEnKS QS always requires less propagations of the ensemble, which improves the computational cost. However, for L > 40 with the L95 model, the quasi-static approach cannot sustain the nonlinearity any-conditions (Bocquet andSakov, 2013, 2014) as long as the DAW length is not excessively long (for a chaotic model), then one concludes that our RMSEs would be systematically equal to or smaller than those reported for any hybrid scheme in [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF].

Figures 8 and9 show the aRMSE of the IEnKS and IEnKS QS for both L95 and L63 as a function of S and L. The smoothing and filtering performance of the IEnKS increases for small values of L, S then decreases for high values of L,S. This is due to the appearances of local minima. As noted by [START_REF] Goodliff | Comparing hybrid data assimilation methods on the Lorenz 1963 model with increasing non-linearity[END_REF] with L63, the QS variant allows one to reach much longer DAWs, and improves the performance. However, some limits to the L95 model method are visible with when, for L = 50, best smoothing aRMSEs are reached for S < 50. However, the IEnKS QS filtering performance is invariant with L and improves with S as in the 4D-Var filtering performance shown in Fig. 3. As suggested by [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF], one may be tempted to estimate the useful DAW length from past observations beyond which the performance gain is negligible. However, their study estimated useful DAW length in a one-cycle 4D-Var context with a focus on the filtering RMSE. Hence, this estimation is not directly relevant for a cycled IEnKS QS . By contrast, in this study, a lot of observations have already been assimi- Quasi-static minimizations lead slowly but surely to the global minimum by repeated cost function minimizations. As the DAW length L is gradually increased, the starting point of the minimization remains in the global minimum basin of attraction. For most S, L couples, the quasi-static IEnKS turns out to be a more accurate substitute for the multiple data assimilation IEnKS (IEnKS-MDA).

Unfortunately, this method (IEnKS QS ) adds an outer loop which could significantly increase the numerical cost. Precision on intermediate minima being superfluous, one can limit the intermediate Gauss-Newton number of loops. The unavoidable space increments required to minimize the nonquadratic cost function are thus reported in time in the quasiconvergent IEnKS (IEnKS QC ).

We did not focus on the applicability of the methods to high-dimensional and imperfect models. In particular, we considered very long DAWs, which, even if of high mathematical interest or for low-order reliable models, is less relevant for significantly noisy models. However, we know from [START_REF] Swanson | Four-dimensional variational assimilation and predictability in a quasi-geostrophic model[END_REF], that the perfect model results are expected to extend to the imperfect model case provided that the growth rate of the model error is similar to that of the leading Lyapunov vectors of the model. This is likely to also apply to a (strong-constraint) IEnKS QS . Finally, an extension to this work would therefore consist in investigating the same ideas but using a weak constraint 4D-Var and IEnKS [START_REF] Trémolet | Accounting for an imperfect model in 4D-Var[END_REF]Sakov and Bocquet, 2018;Sakov et al., 2018).

Data availability. No data sets were used in this article.

those ensemble propagations have different behavior for the L95 and L63 models when the minimizations fail to find the global minimum. For the L95 model, the number of ensemble propagations is high meaning that the minimization takes a lot of iterations and fails to converge. For the L63 model, the number of ensemble propagations is low indicating that minimizations converge but to a non-global extrema.

Conclusions

In this paper, we have extended the study of [START_REF] Pires | On extending the limits of variational assimilation in nonlinear chaotic systems[END_REF] on quasi-static variational data assimilation, focused on 4D-Var technique, to cycled data assimilation schemes and specifically four-dimensional nonlinear ensemble variational techniques, an exemplar of this being the iterative ensemble Kalman smoother (IEnKS).

The long term impact of cycling was first theoretically investigated in a linear context for 4D-Var and the IEnKS, then numerically for the IEnKS in a nonlinear context. The way information is propagated between data assimilation cycles indeed explains the difference between 4D-Var and the IEnKS. Both reveal performance improvements with the DAW parameter S, which counts the number of observation vectors within the DAW, as well as the time shift between cycles. This is a consequence of the Gaussian background approximation; the larger S is, the less the assimilation relies on it.

However, it is observed that this improvement has a limit in the chaotic, perfect model case. The cost function global minimum basin of attraction appears to shrink with increasing L. This causes the Gauss-Newton procedure to miss the cost function global minimum, which deteriorates the assimilation performance. where c k+1 is a constant independent from x (k+1)S and w (k+1)S . Hence G x (k+1)S |y (k+1)S+L:K is Gaussian with moments given by Eqs. (B1) and (B2).

The conditional variance Eq. ( B2) is therefore related to the IEnKS performance by the total law of expectation:

P IEnKS kS = E x a kS -x kS 2 ,
= E E x a kS -x kS 2 |y kS+L:K ,

= E V x kS |y kS+L:K , = E X a kS 2 . (B7)
Then, from Eq. ( 16) we get the recurrence relation: Thus, X a kS is not random and P IEnKS kS = X a kS 2 . Equation (B8) shows that the sequence of inverse IEnKS eMSEs is arithmetico-geometric:

P IEnKS kS -1 α -2S P IEnKS (k-1)S -1 + L K , (B9a)
P IEnKS -S -1 = b -1 , (B9b)
where the notation Eq. (A7) has been used. Properties of arithmetico-geometric sequences allow one to obtain the IEnKS asymptotic eMSE:

P IEnKS ∞S = 0 if |α| ≤ 1, r h 2 α 2L α 2 -1 α 2 otherwise, ( B10
)
and the generalization to asymptotic eMSEs with lag Ll is straightforward:

P IEnKS ∞S+l = 0 if |α| ≤ 1, r h 2 α 2(L-l) α 2 -1 α 2 otherwise. ( B11
)
Let us now show that the IEnKS eMSE is optimal. Let x a kS (y kS+L:K ) be the 4D-Var analysis or any other function of y kS+L:K . A bias-variance decomposition (e.g., [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF] of this estimator yields E x kS ,y kS+L:K x kS -x a kS 2 = E y kS+L:K V x kS x kS |y kS+L:K + E y kS+L:K E x kS x kS |y kS+L:K -x a In the multivariate, diagonal case the algebra can be conducted on each direction independently. Thus, the eMSE in this case is the sum of the univariate eMSEs of each direction.

Appendix C: Expression of the averaged cost function

The IEnKS averaged cost function J ∞S is the N goes to ∞ limit of

1 N N-1 k = 0 J w; y kS+L:kS+K , E b kS = 1 2 w 2 (C1) + 1 2N L l=K N-1 k = 0 y kS+l -H • M l x b kS + X b kS w 2 R -1 .
Expanding the squared norm around H • M l (x kS ) using Eq. (1) gives

1 2 y kS+l -H • M l x b kS + X b kS w 2 R -1 = 1 2 ε kS+l 2 R -1 + 1 2 δy kS+l 2 R -1 + ε T kS+l R -1 δy kS+l , (C2)
where δy kS+l = H • M l (x kS ) -H • M l x b kS + X b kS w . We assume that random variables exist ε ∞S , x ∞S , E a ∞S whose distribution is invariant and ergodic with respect to the shift transformation:

T : ε kS , x kS , E a kS -→ ε (k+1)S , x (k+1)S , E a (k+1)S . (C3) Then because the (ε kS ) k are mutually independent, independent from the x kS , E a kS k and identically distributed, p ε ∞S , x ∞S , E a ∞S = p (ε 0 ) p x ∞S , E a ∞S . By Birkhoff's ergodic theorem (see [START_REF] Walters | An Introduction to Ergodic Theory[END_REF] we get 1 Special issue statement. This article is part of the special issue "Numerical modeling, predictability and data assimilation in weather, ocean and climate: A special issue honoring the legacy of Anna Trevisan (1946Trevisan ( -2016)". It is not associated with a conference.

N N-1 k = 0 ε kS+l 2 R -1 → N→∞ E ε 0 2 R -1 = d, (

Figure 3 .

 3 Figure3. The asymptotic smoothing (a) and filtering (b) eMSEs of 4D-Var and the IEnKS. eMSEs are displayed as functions of L, S with their components on the stable and unstable subspaces. The L parameter is on the abscissa axis for 4D-Var and on the ordinate axis for the IEnKS. The S parameter is on the abscissa axis for the IEnKS and on the ordinate axis for 4D-Var.

Figure 4 .

 4 Figure 4. The asymptotic eMSEs as a function of the lag (lag = 0 is the filtering performance and lag = L is the smoothing performance). The superposing curves have been slightly translated for better readability, S = L = 5, h = b = r = 1.

Figure 5 .

 5 Figure5. Cost functions of the IEnKS projected in one direction of the analyzed ensemble (hence centered and normalized) with various DAW parameters S, L. A quasi-static minimization can be visualized from these panels. It begins with the bottom-left cost function; the orange square is the starting point and the green dot is at the minimum. From the bottom-left cost function up to the top-right cost function, batches of 9 then 10 observation vectors are progressively added to the DAW, and the minimizer (green dot) is updated accordingly.

Figure 6 .

 6 Figure6. Smoothing aRMSE (a) and filtering aRMSE (b) of a Gauss-Newton IEnKS L = S as a function of S with the L95 model over 5 × 10 5 cycles (logarithmic scale). We used the finite-size version(Bocquet, 2011) to account for sampling errors and avoid the need for inflation, H = R = B = I m , n = 20.

Figure 8 .Figure 9 .

 89 Figure 8. IEnKS (lower triangles) and IEnKS QS (upper triangles, N Q = S) smoothing and filtering aRMSEs as a function of L and S with the L95 model. The L parameter is on the abscissa axis for the IEnKS and on the ordinate axis for the IEnKS QS . The S parameter is on the abscissa axis for the IEnKS QS and on the ordinate axis for the IEnKS. For readability, smoothing RMSE beyond 0.10 and filtering RMSE beyond 0.20 are in the same color and the scale is logarithmic.

Figure 10 .Figure 11 .

 1011 Figure10. IEnKS-MDA (S = 1), IEnKS QS and IEnKS (S = L, N Q = S) filtering and smoothing aRMSEs and number of ensemble propagations as a function of L for the L95 and L63 models (logarithmic scale). The ensemble propagations are in units of t. For instance, with L = 50, a single ensemble propagation through the DAW counts for 50.

  Eqs. (B1) and (B2) for a k ≥ 0. Because x (k+1)S = M S (x kS ) with M affine, one getsE x (k+1)S |y kS+L:K = M S x a kS = x b (k+1)S , (B4) V x (k+1)S |y kS+L:K = α 2S X a kS 2 = X b (k+1)S ,(B5)using the conditional expectation properties and Eqs. (25a) and (25b). This result together with Eq. (6) and assumption Eq. (13) yields G x (k+1)S |y (k+1)S+L:K (B6)= J w (k+1)S ; y (k+1)S+L:(k+1)S+K , E b (k+1)S + c k+1 ,

  δy ∞S+l = H • M l (x ∞S ) -H • M l x b ∞S + X b ∞S w and x b ∞S ,X b ∞S are respectively the mean and normalized anomaly of M S E a ∞S . Finally, www.nonlin-processes-geophys.net/25/315/2018/ Nonlin. Processes Geophys., 25, 315-334, 2018 A. Fillion et al.: Quasi-static ensemble variational data assimilation Competing interests. The authors declare that they have no conflict of interest.

  and for both models. Before this point, the algorithms fail to find the global minimum and the RMSE is close to the climatological variance. After this point, the algorithms succeed in finding the global minimum and the RMSE is low. For the IEnKS QS with the L95 model, this point can be estimated using results of Sect. 3 by S/S max ≃ 3.6, in remarkable agreement with Fig.11. This point comes later for the IEnKS QC but it demands less ensemble propagations making this algorithm numerically more efficient than the IEnKS QS . However more and the IEnKS QS aRMSE degrades. Hence, the IEnKS-MDA S = 1 still has the best performance. With the L63 model, the IEnKS QS is always better than the IEnKS-MDA suggesting that L could still be increased.Figure11compares the smoothing aRMSE (first column), the filtering aRMSE (second column) and the number of ensemble propagations (third column) of the IEnKS QS and IEnKS QC as a function of the N Q parameter. The IEnKS QS aRMSE quickly decreases after a point for both algorithms

	Smoothing aRMSE	0.047 0.054 0.063 0.079 0.098 0.175		L95		Filtering aRMSE	0.156 0.160 0.163 0.167 0.190		L95		Number of ensemble propagation	28 47 103 186 323 585 919 1902		L95
		0	10	20	30	40	50	0	10	20	30	40	50	0	10	20	30	40	50
				L						L						L
	Smoothing aRMSE	0.036 0.044 0.055 0.079 0.233		L63		Filtering aRMSE	0.069 0.082 0.108 0.244		L63		Number of ensemble propagation	37 145 264 435 872 1316 2254		L63
		0	10	20	30	40	50	0	10	20	30	40	50	0	10	20	30	40	50
				L						L						L
							IEnKS-MDA		IEnKS QS		IEnKS

Acknowledgements. The authors are grateful to Selime Gürol for her comments and suggestions on the manuscript. The authors are also grateful to two reviewers, Carlos Pires and an anonymous reviewer, for raising very interesting questions and for their suggestions, and to Alberto Carrassi, acting as editor, for his remarks and suggestions. CEREA is a member of Institut Pierre-Simon Laplace (IPSL).

Edited by: Alberto Carrassi Reviewed by: Carlos Pires and one anonymous referee

univariate case

The objective of this Appendix is to establish a recurrence relation between the 4D-Var eMSE of each cycle. From this relation we will get an expression for the 4D-Var asymptotic eMSE.

We assume m = 1 and M l (x) = α l x. At the kth cycle, the gradient and Hessian of the 4D-Var cost function Eq. ( 7) are

where y kS+L:kS+K is temporarily denoted y : . Because of the operators' linearity, J is quadratic with respect to x kS and convex. Hence, its minimizer x a kS exists and is characterized by the null gradient equation:

With an exact Taylor expansion around the state x kS we obtain the following:

Note that ∇ 2 J x kS ; y : , x b kS , given by Eq. (A2), is not random and does not depend on x kS , y : , x b kS . That is why it is simply noted ∇ 2 J . Using Eqs. ( 1) and (A4), we have

The random variable x b kS -x kS = α S x a (k-1)S -x (k-1)S is independent from the errors ε kS+K , . . ., ε kS+L . Thus, taking the expectation of the square of Eq. (A5) gives the following expression for the eMSE P 4D-Var kS of 4D-Var at time t kS :

Introducing the notation

we obtain

Thus, P 4D-Var kS k is an arithmetico-geometric sequence. Its limit P 4D-Var ∞S depends on the value of in the following way:

The generalization to any asymptotic eMSE with lag Ll is straightforward:

In the multivariate, diagonal case the algebra can be conducted on each direction independently. The eMSE in this case is the sum of the univariate eMSEs of each direction.

Appendix B: Performance of the IEnKS in the linear, univariate case

The objective of this appendix is to establish a recurrence relation between the IEnKS eMSE of each cycle. From this relation we will get an expression for the IEnKS asymptotic eMSE. First, it is proven by recurrence that for all k ≥ 0, G (x kS |y kS+L:K ) is Gaussian with moments

where V is the variance of a random variable and x a kS , X a kS are defined by Eqs. ( 15) and ( 16). Because of the assumptions Eqs. (25c) and (25d) with Eq. ( 13) one gets G (x 0 |y L:K ) = J w 0 ; y L:K , E b 0 + c 0 , (B3)

where c 0 is a constant independent from x 0 and w 0 . Hence, G (x 0 |y L:K ) is Gaussian and its moments are given by Eqs. (B1) and (B2) with k = 0. Now, assume G (x kS |y kS+L:K ) is Gaussian with moments given by