Denis Arzelier
email: arzelier@laas.fr

Florent Bréhard
email: orent.brehard@univ-lille.fr

Mioara Joldes
email: joldes@laas.fr

Marc Mezzarobba

Rounding Error Analysis of an Orbital Collision Probability Evaluation Algorithm

We present an error analysis of an algorithm due to Serra et al. (Journal of Guidance Control and Dynamics, 2016) for computing the orbital collision probability in the short term encounter model. The algorithm reduces the numerical computation of the collision probability to that of the sum of a series whose coecients are produced by a linear recurrence relation, and is specically designed to avoid cancellation issues in the evaluation of the sum. The authors derived a bound on the method error arising from the truncation of the series, and observed experimentally that the computation of the terms of the sum is numerically stable, but did not study the evaluation error. Here we give a rigorous bound on the accumulated rounding error when Serra et al.'s algorithm is implemented in oating-point arithmetic. For a unit roundo u and a truncation order N , the bound is of the form (N + A)u + o(u) where A is an explicit constant depending on the problem parameters and o(u) stands for explicitly bounded small terms compared to u. Our analysis is based on the observation that the generating series of the errors aecting each individual term is solution to a perturbed form of a dierential equation satised by the Laplace transform of a function related to the collision probability.

1

Introduction

Due to the drastic increase in the space debris number during the last decades, collision avoidance has become a usual and necessary procedure for many active satellites. The uncertainty aecting the measured data characterizing an encounter is a strong incentive to rely on a probability of collision as the decision variable to trigger a possible avoidance maneuver. When modeling conjunctions, two main paradigms the short-term and the long-term encounters are widely accepted and implemented in the eld of orbital collision risk assessment [START_REF] Chan | Spacecraft Collision Probability[END_REF]. The rst one is most frequently used in practice and assumes that the relative velocity between the two objects is suciently high, so that the encounter time is short. In this framework, the orbital collision probability is modeled as a 2-D integral on a disk, which can be eciently evaluated using an approximation by a power series. The corresponding algorithm of Serra et al. [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] was implemented in Floating Point (FP) arithmetic and has been used in practice by the French Space Agency (CNES) for ground space surveillance operations. More recently, an on-board implementation was successfully tested on an experimental satellite [START_REF] Thomassin | Asteria : Autonomous collision risks management[END_REF].

While the parameters of this algorithm are only estimations of physical quantities, it is however important to provide guaranties about the accuracy and reliability of its numerical implementation. This can be seen by analogy to the need for accurate implementations of special functions (like erf, Airy, Bessel, etc.) used in calculations of other physical phenomena. One would like to estimate and bound independently the numerical evaluation and truncation error for such a mathematical function, compared to other model errors.

With this in mind, the mathematics for this problem were well-studied (truncation error bounds, positivity of the coecients), but the round-o error analysis was so-far ignored. This was probably due to the diculty of the task, since it involves a loop, which implements the evaluation of a linearly recursive sequence.

It is known that the naive rounding error analysis of such recurrences can result in overestimation of the bounds [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] because rounding errors generated in the evaluation of the loop typically cancel out to a large extent, instead of purely adding up. Taking into account this phenomenon usually involves a careful study of the propagation of local errors in following steps of the algorithm, implying complicated manipulations of nested sums and yielding opaque expressions. To alleviate this issue, the main idea of the recent work [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] is to encode as generating series both the sequence of local errors committed at each step and that of global errors resulting from the accumulation of local errors. While far from classical in the context of rounding-error analysis, this technique proves to be very adequate for studying algorithms which originate from numerical methods implementing truncated series approximations, with coecients satisfying linear recurrences.

Dierent alternative methods are now briey recalled. Firstly, unrolling a linear recurrence can be seen as a special case of solving a triangular (banded Toeplitz) system of linear equations. Therefore, a rst result based on [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]Chap. 8], bounds the maximum relative error for evaluating n terms of an order m recurrence by the product of mu

1-mu

with the condition number of the associated n by n matrix. In this sense, in [START_REF] Barrio | [END_REF], a more rened analysis gives a rst-order bound (meaning that the terms of order O(u 2 ) are omitted). However, one has to resort to more complicated formulas, expressed in terms of quantities that may be dicult to estimate (inverting the associated triangular matrix, computing the so-called reverse homogeneous recurrence for instance).

A complementary class of approaches concerns the use of static error analyzers, which automatically provide sound (and often formally-proven) error bounds on FP rounding errors (see for instance [START_REF] Appel | Vcoat2: Floating-point error analysis in coq[END_REF][START_REF] Solovyev | Rigorous estimation of oating-point round-o errors with symbolic Taylor expansions[END_REF][START_REF] Titolo | An abstract interpretation framework for the round-o error analysis of oating-point programs[END_REF] and references therein for existing software). While these tools are aimed at generic numerical codes, they are not very ecient for handling a very large number of loop iterations due to the intrinsic high depth of the expression graph.

For instance, one of the currently fastest tools, SATIRE [START_REF] Das | Scalable yet rigorous oating-point error analysis[END_REF], reports a minimum execution time of 50s for unrolling 70 iterations of the Lorenz system.

By comparison, the algorithm analyzed in this article sometimes requires hundreds (or even thousands) of iterations. Furthermore, the parameters involved have rather high ranges and we would like an error bound which depends explicitly on these parameters, without additional runs of the program.

Finally, let us also mention that a basic automatic evaluation in interval arithmetic highly overestimates the bounds as the iterations directly reuse the previously-computed values.

All in all, we believe that the adaptation of [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] to this particular algorithm, oers a good remedy to these limitations and that the mathematical tools employed herein may be of interest to the rounding-analysis spectrum of methods.

The structure and contributions of this article are as follows. Firstly, we recall in Section 2 the description of a Laplace transform technique from [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF],

which proves that the terms of the recursive sequence implemented as the main loop of the algorithm are the coecients of a series solution to a simple rst-order Linear Dierential Equation (LDE). From that, the main contribution of this article is to make heavy use of this equation to interpret the individual rounding errors on each term as the coecients of another series whose analytic behavior is essentially similar, up to a factor proportional to the roundo unit. This is the key point for deriving realistic total relative error bounds. To do that, a preliminary step in Section 3 is a classical rounding error analysis for the loop-independent parameters and the body of the loop (local errors).

Then, we bound the global errors accumulating when executing the main loop in Section 4. There, we make the key observation that the generating series in the Laplace plane associated to global errors is solution to the same previously mentioned LDE, but with an inhomogeneous term generated by the local errors. Working with series in the Laplace plane allows for suciently simple closed-form formulas, even if it sometimes means performing some crafty term-by-term majorizations. This technique allows for the computation of explicit a priori rounding error bounds depending on the input parameters, without any additional runs or restrictive imposed ranges on the parameters. Given a unit roundo u and a truncation order N , the bound proposed in Theorem 1 is of the form (N + A)u + o(u), where A is an explicit constant depending on the parameters and o(u) stands for explicitly bounded small terms compared to u.

Finally, practical aspects are considered: in Section 5 we comment on how the analyzed implementation simulates an increased exponent range, as to avoid overows in practice; then, a numerical validation of the quality of the bound is provided on a range of examples in Section 6.

Computing the Orbital Collision Probability

In this section, the main steps of the mathematical derivations used to build the reviewed algorithm are briey reminded. This is particularly useful since the proposed rounding error analysis partly exploits the very same ingredients.

The short-term encounter model

The short-term encounter model (whose complete set of assumptions is recalled in detail in [START_REF] Chan | Spacecraft Collision Probability[END_REF] or [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF]) for the computation of the probability of collision between two spherical objects mainly consists in assuming that the relative trajectory is a straight line during the encounter and in projecting it onto the encounter plane dened to be perpendicular to the relative velocity vector. Let (x, y) denote the mean coordinates of the relative position of the secondary object with respect to the primary object in the encounter frame (see [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] for its denition). The relative position uncertainty is described by the following bivariate Gaussian density function

ρ(x, y) = 1 2πσ x σ y exp - 1 2 (x -x m ) 2 σ 2 x + (y -y m ) 2 σ 2 y ,
where (x m , y m ) is the mean position of the secondary object relative to the primary object in the covariance frame and σ x , σ y ∈ R + * are the standard deviations of the relative coordinates in the encounter plane. The probability of collision is then given by a two-dimensional integral parameterized by the radius R of the combined spherical object:

P(R) =
x 2 +y 2 ⩽R ρ(x, y)dxdy.

(

) 1
This integral is the cumulative density function of the random variable Ξ =

X 2 + Y 2 (i.e., P(R) = Pr {Ξ ⩽ R 2 }) where X ▷ N (x m , σ 2
x ), Y ▷ N (y m , σ 2 y ) are independent normal random variables. By rescaling Ξ, we get P(R) = g [START_REF] Appel | Vcoat2: Floating-point error analysis in coq[END_REF] where the function g : R + → R + is dened as g(ξ) = P(R √ ξ).

Several methods for computing the integral (1) have been proposed in the aerospace literature (see for instance [START_REF] García-Pelayo | Series for collision probability in short-encounter model[END_REF][START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] and references therein). Here we focus on the algorithm of [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF], where one obtains a convergent series expansion of g(ξ) by considering its Laplace transform. Various versions of this idea have been applied to both central and non-central quadratic forms; see for instance [START_REF] Mathai | Quadratic forms in random variable: Theory and Applications[END_REF]Chapter 4] and references therein.

Algorithm

Due to the cancellation phenomenon occurring when summing the terms of a series of dierent signs and similar magnitude, the direct evaluation of the power series expansion of g(ξ) is only practical for small values of R 2 . The idea of [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] to remedy this, inspired by [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF], is to introduce a preconditionner of the form Π(ξ) = exp(pξR 2 ) and consider the expansion of the function f = Πg, which then has nonnegative coecients. Assuming without loss of generality that 0 < σ y ⩽ σ x , a good choice is p = 1/(2σ 2 y ).

We now summarize how the authors of [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] obtain a linear-time algorithm for computing the rst N terms of the series f . We use the following notation:

p = 1 2σ 2 y , ϕ = 1 - σ 2 y σ 2 x , ω x = x 2 m 4σ 4 x , ω y = y 2 m 4σ 4
y , and observe that 0 ⩽ ϕ < 1, ω x ⩾ 0, ω y ⩾ 0.

The rst step is to compute the Laplace transform of the preconditioned function in closed form. For |λ| > p, one has

L f (λ) = L g (λ -pR 2 ) = πR 2 ρ(0, 0) exp ωyR 2 λ + ωxR 2 λ-pϕR 2 λ(λ -pϕR 2 )(λ -pR 2 ) .
The power series expansion of f (ξ) is the termwise inverse Laplace transform of the expansion at innity of L f (λ) (see [START_REF] Widder | An introduction to transform theory[END_REF]Chap. 9], [16, Chap. 2.14]), and the coecients of the latter are the same as those of the Taylor expansion of

L f (λ -1 ) = πR 2 ρ(0, 0)λ 2 exp ω y R 2 λ -ωx pϕ - ωx pϕ(pϕR 2 λ-1) 1 -pϕR 2 λ(1 -pR 2 λ) . Since the rst two coecients of L f (λ -1 ) are zero, let f (λ) = λ -2 L f (λ -1 ). Letting f (ξ) = ∞ n=0 c n ξ n+1 (note the n + 1), the series expansion of f reads f (λ) = ∞ n=0 c n (n + 1)!λ n .
(

) 2
The second step is to derive a linear recurrence relation satised by the coecients c n . For this, one uses the fact that f (λ) is solution to a LDE. Indeed, starting from the denition of f and taking logarithmic derivatives, one has

f ′ (λ) = φ(λ) f (λ), f (0) = πR 2 ρ(0, 0), (3)
φ(λ) = ω y R 2 + pϕR 2 2(1-pϕR 2 λ) + pR 2 1-pR 2 λ + ω x R 2 (1-pϕR 2 λ) 2 = P (λ) Q(λ) with Q(λ) = (1-pϕR 2 λ) 2 (1-pR 2 λ). (4)
The coecients of the polynomials P and Q alternate in sign: we write P (λ) = P 0 -

P 1 λ + P 2 λ 2 -P 3 λ 3 and Q(λ) = 1 -Q 1 λ + Q 2 λ 2 -Q 3 λ 3 where P i , Q i ⩾ 0.
Lemma 1. The sequence (c n ) satises the linear recurrence

nc n - Q 1 (n -1) + P 0 n + 1 c n-1 + Q 2 (n -2) + P 1 (n + 1)n c n-2 - Q 3 (n -3) + P 2 (n + 1)n(n -1) c n-3 + P 3 (n + 1)n(n -1)(n -2) c n-4 = 0, (5)
for n ⩾ 4, with initial terms c 0 , . . . , c 3 given in Algorithm 1.

Proof sketch. The LDE (3) is equivalent to

Q(λ)λ f ′ (λ) -P (λ)λ f (λ) = 0. (6)
One can check that, for any series f satisfying (2), one has

λ k f (λ) = +∞ n=k c n-k (n + 1) . . . (n -k + 2) (n + 1)!λ n , λ( f ) ′ (λ) = +∞ n=0 (nc n )(n + 1)!λ n . (7)
Using these identities repeatedly, one obtains

Q(λ)λ f ′ (λ) -P (λ)λ f (λ) = ∞ n=0 F n (c)(n + 1)!λ n ,
where F n (c) is exactly the left-hand side of ( 5), with the additional convention that c -1 = c -2 = c -3 = 0 and terms whose denominator vanishes are ignored. It follows using ( 6) that F n (c) = 0 for all n ⩾ 0. For n ⩾ 4, this gives the desired recurrence. For n = 1, 2, 3, one obtains the expressions for c 1 , c 2 , c 3 appearing in the algorithm, and similarly for c 0 since c 0 = f (0) = πR 2 ρ(0, 0).

Summarizing, we have

P(R) = g(1) = exp -pR 2 ∞ n=0 c n ,
where the coecients c n are given by Lemma 1. Algorithm 1 is a procedure for evaluating this expression. Our next goal is to analyse the eect of rounding errors on this procedure.

Local Rounding Error Bounds

In this section, we describe the employed FP setting and provide the local error analysis.

FP arithmetic setting

We assume that Algorithm 1 is implemented in radix- Input: Parameters: σ x , σ y , x m , y m ; combined object radius: R; number of terms: N . Output: P 0:N truncated series approximation of P.

1: p = 1 2σ 2 y ; ϕ = 1 -σy σx 2 ; ω x = x 2 m 4σ 4 x ; ω y = y 2 m 4σ 4 y ; 2: Q 1 = pR 2 (2ϕ + 1); Q 2 = p 2 R 4 ϕ(ϕ + 2); Q 3 = p 3 R 6 ϕ 2 ; 3: P 0 = p ϕ 2 + 1 + ω x + ω y R 2 ; 4: P 1 = pϕ(ϕ+5) 2 + ω x + ω y (2ϕ + 1) pR 4 ; 5: P 2 = 3 2 p + ω y (ϕ + 2) p 2 R 6 ϕ; 6: P 3 = p 3 ω y R 8 ϕ 2 ; 7: c 0 = 1 2σxσy exp -1 2 x 2 m σ 2 x + y 2 m σ 2 y R 2 ; 8: c 1 = P 0 2 c 0 ; c 2 = Q 1 +P 0 6 c 1 -P 1 12 c 0 ; 9: c 3 = 2Q 1 +P 0 12 c 2 -Q 2 +P 1 36 c 1 + P 2 72 c 0 ; 10: s = c 0 + c 1 + c 2 + c 3 11: for n = 4 to N -1 do 12: c n = Q 1 (n-1)+P 0 (n+1)n c n-1 -Q 2 (n-2)+P 1 (n+1)n 2 c n-2 + Q 3 (n-3)+P 2 (n+1)n 2 (n-1) c n-3 - P 3 (n+1)n 2 (n-1)(n-2) c n-4 ; 13: s = s + c n ;
14: end for 15: return P 0:N = exp -pR 2 s. tie-breaking rule). In particular, our arithmetic obeys the standard error models

[8, Chap. 2.2.] RN(x) = x(1 + r 1 ) = x/(1 + r 2 ), with |r 1 |, |r 2 | ⩽ u, (8)
with rounding unit u = 2 -t . In addition, multiplications by powers of two are exact. This setting correctly models real-life IEEE-754 arithmetic provided that no overows or underows occur. It turns out that, when implemented in binary64 arithmetic, Algorithm 1 can easily encounter overows for realistic values of the input. We comment in Section 5 on how the implementation by the authors of [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF] simulates an increased exponent range, making the assumption of an unbounded exponent range legitimate for the error analysis.

For deniteness, we also assume that composite expressions are evaluated from left to right: for instance, a + bcd is computed as RN(a + RN(RN(bc)d)), the power operation is implemented according to the formulas

x 2 = x•x, x 3 = x • x 2 , x 4 = (x 2 ) 2 , x 6 = (x 2 ) 3 , and x 8 = (x 4 ) 2 ,
the exponential function used at steps 7 and 15 is faithfully rounded, implying that the corresponding relative error is bounded by 2u.

These assumptions are not critical and our bounds easily adapt to slightly dierent implementations.

We denote by x the computed value of a quantity x. To express the relation between x and x, we use the θ k and γ k notation of [START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF]Chap. 3]. In short, each occurrence of the symbol θ k denotes a potentially dierent quantity of the form θ k = k i=1 (1 + r i ) ±1 -1 with |r i | ⩽ u for all i. Assuming ku < 1, one has |θ k | ⩽ γ k where γ k is dened as ku/(1 -ku) and satises in particular γ k = ku + O(u 2 ) as u → 0.

Bounds for loop-independent parameters

Let us rst bound the rounding errors occurring in lines 17 of Algorithm 1.

For instance, using Equation ( 8), one has

p = RN(1/RN(σ y • σ y ))/2 = p(1 + θ 2 ), (9)
which gives an absolute error bound of |p -p| ⩽ γ 2 p. Similar bounds for the other parameters are summarized below.

Param. p ω x ω y ϕ Q 1 Q 2 Q 3 Abs. Err. γ 2 p γ 5 ω x γ 5 ω y γ 4 γ 9 Q ♯ 1 γ 18 Q ♯ 2 γ 24 Q ♯ 3 Param. P 0 P 1 P 2 P 3 c 0 e -pR 2
Abs. Err. 1, where

γ 10 P ♯ 0 γ 18 P ♯ 1 γ 27 P ♯ 2 γ 32 P ♯ 3 e 0 c 0 τ e -pR 2
P ♯ i := P i {ϕ ← 1}, Q ♯ i := Q i {ϕ ← 1} denote the values for ϕ = 1 of the P i and Q i : Q ♯ 1 = 3pR 2 , Q ♯ 2 = 3p 2 R 4 , Q ♯ 3 = p 3 R 6 , P ♯ 0 = 3 2 p + ω x + ω y R 2 , P ♯ 1 = (3p + ω x + 3ω y ) pR 4 , P ♯ 2 = 3 2 p + 3ω y p 2 R 6 , P ♯ 3 = p 3 ω y R 8 , e 0 = exp 1 2 x 2 m σ 2 x + y 2 m σ 2 y γ 4 (1 + γ 6 ) -1, τ = exp pR 2 γ 2 (1 + γ 2 ) -1. (10)
Proof. Similarly to Equation ( 9) one obtains ωx = ω x (1 + θ 5 ) and ωy = ω y (1 + θ 5 ).

Concerning ϕ, rstly observe that since 0 < σ y ⩽ σ x , one has 0 < σ y /σ x ⩽ 1 and because RN preserves inequalities, 0 < ã := RN(RN(σ y /σ x )

• RN(σ y /σ x )) ⩽ 1, implying that 0 ⩽ b := 1 -ã < 1. Now, φ = RN(b), so that 0 ⩽ φ ⩽ 1 and |ϕ -φ| ⩽ |ϕ -b| + |b -φ| ⩽ |σ 2 y /σ 2 x -ã| + u ⩽ γ 3 + γ 1 ⩽ γ 4 .
Regarding Q 1 , one has Q1 = pR 2 (2 φ + 1)(1 + θ 6 ), and hence, using the previous bounds on ϕ and |ϕ -φ|,

|Q 1 -Q1 | = |pR 2 (2ϕ + 1) -pR 2 (2 φ + 1)(1 + θ 6 )| ⩽ 2pR 2 |ϕ -φ| + pR 2 |(2 φ + 1)|γ 6 ⩽ 2pR 2 γ 4 + 3pR 2 γ 6 ⩽ 3pR 2 γ 9 .
The last inequality is readily obtained by applying the rules given in Lemma

Local error analysis

Let us now turn to the computation of c n (steps 7 12 of Algorithm 1). We denote by cn the computed value of c n , and we call local absolute error on c n the absolute rounding error ε n generated at the corresponding step of the algorithm. In other words, for n ⩾ 4, we set

ε n := cn - Q 1 (n -1) + P 0 (n + 1)n cn-1 - Q 2 (n -2) + P 1 (n + 1)n 2 cn-2 + Q 3 (n -3) + P 2 (n + 1)n 2 (n -1) cn-3 - P 3 (n + 1)n 2 (n -1)(n -2) cn-4 (11)
(where all operations are mathematically exact). We then have the following bound on |ε n |.

Proposition 2. The local error introduced at iteration n at step 12 of Algorithm 1 satises

|ε n | ⩽ γ Q ♯ 1 (n -1) + P ♯ 0 (n + 1)n |c n-1 | + Q ♯ 2 (n -2) + P ♯ 1 (n + 1)n 2 |c n-2 | + Q ♯ 3 (n -3) + P ♯ 2 (n + 1)n 2 (n -1) |c n-3 | + P ♯ 3 (n + 1)n 2 (n -1)(n -2) |c n-4 | ,
where γ = γ 40 .

Proof. The coecient c n is computed as cn = (((t

1 + t 2 )(1 + θ 1 ) + t 3 )(1 + θ 1 ) + t 4 )(1 + θ 1 ) with t i = (-1) i+1 ( Qi (n -i)(1 + θ 1 ) + Pi-1 )(1 + θ 1 ) d i (n)(1 + θ i ) -1 cn-i (1 + θ 2 ), where d 1 (n) = (n+1)n, d 2 (n) = (n+1)n 2 , .
. . are the denominators appearing in [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF], and Q 4 = 0. We thus have

cn = 4 i=1 ( Qi (n -i)(1 + θ 9 ) + Pi-1 (1 + θ 8 ) d i (n) .
Substituting into [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF], we obtain

ε n = 4 i=1 (-1) i+1 cn-i d i (n) ( Qi -Q i + Qi θ 9 )(n -i) + ( Pi-1 -P i-1 + Pi-1 θ 8 ) . According to Table 1, we have | Qi -Q i | ⩽ Q ♯ i γ 24 (with the convention that Q ♯ 4 = 0) and | Qi θ 9 | ⩽ (1 + γ 24 )Q ♯ i γ 9 , so that | Qi -Q i + Qi θ 9 | ⩽ Q ♯ i γ 33 . Similarly, we have | Pi-1 -P i-1 + Pi-1 θ 8 | ⩽ P ♯
i-1 γ 40 , and the result follows.

Since the formulas used for computing c 1 , c 2 , c 3 correspond to truncated instances of the recurrence, (11) also applies for ε 1 , ε 2 , ε 3 , if terms with a zero denominator are ignored. With this convention the bound from Prop. 2 holds for all n ⩾ 1.

Global Rounding Error Bounds

Let us apply the generating series approach of [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] to the rounding error analysis of the main loop.

Global error modeling

The local errors ε n build up and lead to a global (absolute) error

δ n := c n -cn ,
that is the main quantity we need to control. For doing so it is convenient to encode the sequences (δ n ) and (ε n ) as coecients in the generating series δ(λ) = +∞ n=0 (n + 1)!δ n λ n and ε(λ) = +∞ n=0 (n + 1)!ε n λ n in the same manner as (2).

From Equations ( 5) and ( 11), we have

nδ n = Q 1 (n -1) + P 0 (n + 1) δ n-1 - Q 2 (n -2) + P 1 (n + 1)n δ n-2 + Q 3 (n -3) + P 2 (n + 1)n(n -1) δ n-3 - P 3 (n + 1)n(n -1)(n -2)
δ n-4 -nε n .

After multiplying this relation by λ n and summing over n, we obtain, using the identities [START_REF] Gawronski | Reduced cancellation in the evaluation of entire functions and applications to the error function[END_REF], a LDE satised by the series δ:

Q(λ)( δ) ′ (λ) -P (λ) δ(λ) = ε′ (λ). (12)
Comparing with Equation ( 6), we see that δ satises the same rst-order LDE as the generating series f of the (exact) coecients c n , except for the right hand side which now depends on the local errors ε n . Only bounds are available for these, so we need to work with dierential inequalities.

ε′ (λ) ≪ γ Q ♯ (λ)φ(λ) + P ♯ (λ) f (λ)+ Q ♯ (λ)| δ| ′ (λ) + P ♯ (λ)| δ|(λ) , (13)
with

Q ♯ (λ) = Q ♯ 1 λ + Q ♯ 2 λ 2 + Q ♯ 3 λ 3 and P ♯ (λ) = P ♯ 0 +P ♯ 1 λ+P ♯ 2 λ 2 +P ♯ 3 λ 3 +P ♯ 4 λ 4 .
Proof. It follows by using the inequality |c n | ⩽ c n + |δ n | in the bound on ε n obtained from Proposition 2, which is to be multiplied by nλ n , and summed over n.

Located in the Laplace plane, Eqs. ( 12) and ( 13) allow for deriving bounds on δ(λ). They are obtained as solutions of order-1 LDE. But there is still a need to bound the inverse Laplace transform δ(ξ) = ∞ n=0 δ n ξ n+1 . In particular, we need to bound the total sum of absolute rounding errors |δ|(1). This is done by an ad-hoc majorization of convolution terms. For clarity, we provide in Table 2 a synthesis of forthcoming notations in this twofold view.

A simplied bound

For the sake of exposition, we rst prove a simplied, not fully rigorous error bound obtained by neglecting the terms involving γ δ(λ) in [START_REF] Solovyev | Rigorous estimation of oating-point round-o errors with symbolic Taylor expansions[END_REF], which are of order O(u 2 ):

Q(λ) δ′ (λ) -P (λ) δ(λ) ≪ γ Q ♯ (λ)φ(λ) + P ♯ (λ) f (λ).
To further simplify this equation, we denote φ ♯ (λ) := φ(λ){ϕ ← 1} and use 0

≪ φ(λ) ≪ φ ♯ (λ) and 0 ≪ Q(λ) -1 ≪ Q(λ) -1 {ϕ ← 1} = (1 -pR 2 λ) -3
(this follows directly from Equation ( 4)) to obtain

| δ| ′ (λ) ≪ φ(λ)| δ|(λ) + γ ψ(λ) f (λ), with ψ(λ) := Q ♯ (λ)φ ♯ (λ) + P ♯ (λ) (1 -pR 2 λ) 3 ≫ 0. ( 14
)
Solving this dierential inequality gives rise to the following simplied bound.

Proposition 4. Under the simplied model above, the total rounding error accumulated while computing f (1) satises

+∞ n=0 |c n -c n | = |δ|(1) ⩽ (e 0 + γC)f (1), 'Real' (ξ)
'Laplace at ∞' (λ) with e 0 given in [START_REF] Mathai | Quadratic forms in random variable: Theory and Applications[END_REF] Proof. Since all the series on the right-hand side of ( 14) have nonnegative coecients, Lemma 6.5 in [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] implies that | δ|(λ) ≪ ∆(λ) where ∆(λ) is the solution with ∆(0) = e 0 f (0

f (ξ) = ∞ n=0 c n ξ n+1 f (λ) = ∞ n=0 c n (n + 1)!λ n f (λ) = λ -2 L f (λ -1 ) Q f ′ -P f = 0 δ(ξ) = ∞ n=0 δ n ξ n+1 δ(λ) = ∞ n=0 δ n (n + 1)!λ n δ(λ) = λ -2 L δ (λ -1 ) |δ|(1) = ∞ n=0 |c n -cn | Q δ′ -P δ = ε′ ε(ξ) = ∞ n=0 ε n ξ n+1 ε(λ) = ∞ n=0 ε n (n + 1)!λ n ε′ ≪ γ Q ♯ f ′ + P ♯ f + Q ♯ δ′ + P ♯ δ δ(ξ) ≪ ∆(ξ) δ ≪ ∆ ∆′ = φ ∆ + γ Q Q ♯ f ′ + ∆′ + P ♯ f + ∆ ∆ = e 0 f + ∞ k=0 γ k k! Ψ * . . . * Ψ k times * f ∆ = ê f ê = e 0 + (1 + e 0 ) ∞ k=0 γ k k! Ψk Ψ(ξ) ⩽ W (ξ)e pR 2 ξ ∆ = e 0 f + γΨ * f Order 1 approx of ê ê = e 0 + γ Ψ |δ|(1) ⩽ ∆(1)
) ⩾ |δ(0)| of the LDE ∆′ (λ) = φ(λ) ∆(λ) + γ ψ(λ) f (λ). (16)
Using f as a solution of the homogeneous part of ( 16), shows that the series Ψ(ξ) can be bounded as Ψ(ξ) ≪ W (ξ)e pR 2 ξ , where W (ξ) is an explicit polynomial of degree 3 in ξ with nonnegative coecients and 1 0 W (τ )dτ is equal to the constant C dened in [START_REF] Titolo | An abstract interpretation framework for the round-o error analysis of oating-point programs[END_REF]. It follows that

∆(λ) = e 0 + γ Ψ(λ) f (λ), Ψ(λ) := λ 0 ψ(σ)dσ.
(Ψ * f )(1) ⩽ 1 0 W (τ )e pR 2 τ e pR 2 (1-τ ) g(1 -τ )dτ ⩽ e pR 2 g(1) 1 0 W (τ )dτ = Cf (1),
(where the second inequality uses the fact that g is nondecreasing), and therefore |δ|(1) ⩽ ∆(1) = e 0 f (1) + γ (Ψ * f )(1) ⩽ (e 0 + γC)f (1).

1 Maple worksheet used for these computations given in Appendix A.3.

A rigorous bound

For a fully rigorous bound on δ(λ), we consider again the dierential inequalities ( 12) and ( 13). Reasoning as in the previous section, we have

| δ| ′ (λ) ≪ γ Q ♯ (λ) Q(λ) | δ| ′ (λ) + φ(λ) + γ P ♯ (λ) Q(λ) | δ|(λ) + γ Q ♯ (λ)φ(λ) + P ♯ (λ) Q(λ) f (λ),
where the coecients of |δ| and |δ| ′ as well as the inhomogeneous term are series with nonnegative coecients. Since Q ♯ (0) = 0, Lemma 6.5 in [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] applies again and shows that δ(λ) ≪ ∆(λ), with ∆(λ) ≫ 0 satisfying the LDE

∆′ (λ) = φ(λ) ∆(λ) + γ Q(λ) Q ♯ (λ) f ′ (λ) + ∆′ (λ) + P ♯ (λ) f (λ) + ∆(λ) .
Let us write ∆(λ) = ê(λ) f (λ) and use ( 3) to obtain a LDE satised by ê,

where the right-hand side is a positive series:

1 -γ Q ♯ (λ) Q(λ) ê′ (λ) = γ Q ♯ (λ) Q(λ) φ(λ) + P ♯ (λ) Q(λ) (1 + ê(λ)).
Since a(λ) := γQ ♯ (λ)/Q(λ) satises a(0) = 0 and a(λ) ≫ 0,

1 1 -a(λ) = Q(λ) Q(λ) -γQ ♯ (λ) ≫ 0,
by composition of two series with nonnegative coecients. After multiplication by this series, we obtain

ê′ (λ) = γ Q ♯ (λ)φ(λ) + P ♯ (λ) Q(λ) -γQ ♯ (λ) (1 + ê(λ)). ( 17
)
This LDE has several poles due to the perturbation γQ ♯ (λ) of the denominator in the right-hand side. To overcome this additional diculty, we use the following lemma, proved in Appendix A.2, to obtain a unique pole, at the price of a slight increase in the parameter p. This is a key point for adapting the proof of Proposition 4 to the current setting.

Lemma 2. Assuming 7γ < 1, we have

1 Q(λ) -γQ ♯ (λ) ≪ 1 (1 -p + R 2 λ) 3 with p + := p 1 -3 √ 7γ .
This result allows us to bound the solution of LDE ( 17) by the solution of the simpler LDE

ê′ (λ) = γ ψ+ (λ)(1 + ê(λ)), ê(0) = e 0 , (18)
where ψ+ (λ) := ψ(λ){p ← p + } with ψ(λ) dened in [START_REF] Thomassin | Asteria : Autonomous collision risks management[END_REF].

Proposition 5. The terms cn computed in FP arithmetic in the main loop of Algorithm 1 satisfy:

+∞ n=0 |c n -c n | = |δ|(1) ⩽ e 0 + (1+e 0 )e ηpR 2 (e γC + -1) f (1),
with γ = γ 40 , e 0 as in [START_REF] Mathai | Quadratic forms in random variable: Theory and Applications[END_REF], η :

= 3 √ 7γ 1-3 √ 7γ and C + := C{p ← p + }. Proof. Denoting Ψ+ (λ) := λ 0 ψ+ (σ)dσ, LDE (18) gives ê(λ) = (1 + e 0 )e γ Ψ+ (λ) -1 = e 0 + (1 + e 0 ) +∞ k=1 γ k Ψ+ (λ) k k! .
This gives an explicit expression for ∆(λ) = ê(λ) f (λ).

To obtain ∆(ξ) s.t. ∆(λ) = λ -2 L ∆ (λ -1 ), let Ψ + (ξ) be the series such that Ψ+ (λ) = L Ψ + (λ -1 ). Then Laplace transform rules give the following identity of formal power series:

∆(ξ) = e 0 f (ξ) + (1 + e 0 ) +∞ k=1 γ k k! Ψ + * k * f (ξ), (19)
where

Ψ + * k = Ψ + * • • • * Ψ + (k times).
In Lemma 4, Appendix A.2, we prove the remaining inequality:

(Ψ + * k * f )(1) ⩽ e ηpR 2 (C + ) k f (1).

The nal rounding error bound

The truncated series approximation P 0:N = e -pR 2 s = e -pR 2 N -1 n=0 c n of P is obtained by evaluating the sum N -1 n=0 cn in FP arithmetic and by multiplying the result s with e -pR 2

. We call P0:N the FP number returned by Algorithm 1. The following theorem provides a relative rounding error bound w.r.t. P. Adding to this bound a relative truncation error bound on |P 0:N -P|/P derived from [12, III.C] would yield a total error bound on | P0:N -P|/P. Theorem 1. The total rounding error is bounded by

P0:N -P 0:N P ⩽ (1 + γ N )(1 + τ )(1 + e 0 ) 1 + e ηpR 2 (e γC + -1) -1,
with γ = γ 40 , and where the quantities e 0 , τ are dened in Proposition 1, and η, C + are dened in Proposition 5.

The rst-order error approximation in the roundo unit u for this bound on | P0:N -P 0:N |/P is

N + 8 + 2pR 2 + 2x 2 m σ 2 x + 2y 2 m σ 2 y + 40C u, ( 20
)
where 40Cu is the dominant term for large p, R, x m and y m .

Proof. Denote s =

N -1 n=0 cn and s its FP evaluation using N -1 additions.

By Proposition 5, |s -s| ⩽

N -1 n=0 |c n -c n | ⩽ νf (1) where ν := e 0 + (1 + e 0 )e ηpR 2 (e γC + -1). Then |s -s| ⩽ γ N -1 N -1 n=0 |c n | ⩽ γ N -1 (1 + ν)f (1)
| = |e -pR 2 s(1 + τ ′ )(1 + θ 1 ) -e -pR 2 s| ⩽ e -pR 2 s(τ + u + τ u) + |s -s|(1 + τ )(1 + u) ⩽ P τ + u + τ u + (γ N -1 + ν + γ N -1 ν)(1 + τ )(1 + u) ⩽ P (1 + γ N )(1 + τ )(1 + ν) -1 ,
which is exactly the bound claimed by Theorem 1.

Preventing Overows and Underows

Algorithm 1 may be subject to overows and underows, depending on the problem parameters and the number N of terms to be computed. First, according to [12, III.C, Prop. 4], this number N has to be at least 2e(p + ω x + ω y )R 2 to obtain a reasonable approximation of P. Since the c n sum to f (1) = e pR 2 P and P may be close to 1, the use of plain IEEE 754-1985 binary64 FP arithmetic with maximum exponent 1023 may cause overows for examples requiring more than 2e ln (2 1023 ) ≈ 4000 terms.

In the C implementation of [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF], a rescaling strategy is used to prevent overows and underows. At the end of each iteration, if the absolute value of the computed term cn is above A or below A -1 , for A = 2 800 , then the values of cn , cn-1 , cn-2 , cn-3 are rescaled by 2 k for some k so that their absolute values belong to [A -1 , A]. This number k is added to a signed 64-bit integer used to store the current exponent, and the summation of the terms cn keeps track of these intermediate rescalings. Two additional rescalings are also used for c 0 and the nal factor e -pR 2 to prevent underows. Note that this rescaling strategy, where the exponent is stored separately in a 64-bit integer, does not modify the relative rounding error model used in the previous sections.

The following theorem guarantees the absence of overows under reasonable assumptions on the size of input parameters. Theorem 2. We assume that the number N of terms required to approximate P is bounded by N * = 10 8 , and that this bound also holds for pR 2 , ω x R 2 and ω y R 2 . In addition, considering the size of a satellite and its distance to the space debris, we assume 1 ⩽ R ⩽ 10 3 and σ x , σ y , |x m |, |y m | ⩽ 10 6 (all these quantities are expressed in meters).

Then the execution of the C implementation of Algorithm 1 on σ x , σ y , x m , y m , R, N is not subject to overows.

Proof. We prove this property for all the steps of the algorithm, postponing to the end of the proof the additional eect of rounding errors.

Loop-independent parameters. By combining the inequalities assumed for the parameters in the theorem, it is straightforward that all the subexpressions involved in the computation of p, ϕ, ω x , ω y , Q i , P i and c 0 are much smaller than 2 1023 and do therefore not cause any overow.

Evaluation of c n for n ⩾ 1. At the beginning of iteration n, the preceding terms c n-1 , c n-2 , c n-3 , c n-4 are bounded by A = 2 800 in absolute value thanks to the rescaling strategy. A quick analysis shows that each of the four coecients in front of c n-i is bounded by 2N 4 * , and both their numerator and denominator are bounded by 2N 5 * . Hence, their evaluation cannot produce overows. Finally, c n ⩽

3 i=1 Q i + 3 i=0 P i A ⩽ 7N 4 * A ⩽ 2 910 < 2 1023 , so
that no overow can occur.

No overow of the 64-bit exponent. The total sum is bounded by f (1) = e pR 2 P ⩽ e N * ⩽ 2 2 20.5 < 2 2 63 -1 , so that the exponents of all c n and all partial sums t in the 64-bit integer. The nal multiplication by e -pR 2 cannot cause an overow either since the argument is negative.

The eect of rounding errors. The local rounding errors in the constants and in each iteration of the loop were bounded by small constants (Propositions 1 and 2), hence they do not modify signicantly the overow analysis of the rst two items above. The worst-case relative error bound given by Theorem 1 is (crudely) bounded by exp (γ N +8 + (2

• 10 12 + γ 2 + η)N * + γC + ) ⩽ 2 2 61.4
. Although huge, this bound is sucient to prove that the computed sum of the cn is smaller than 2 2 61.4 • 2 2 20.5 < 2 2 63 -1 .

The rigorous underow analysis is slightly more involved and postponed to future work. Roughly speaking, the rescaling strategy prevents underows in the exponentials in c 0 and e pR 2 (which would cause the output to be zero).

Underows can however occur when unrolling the recurrence, but then it means that the neglected terms c n are so small compared to the previous ones that this underow error is smaller than the relative error already computed for the partial sum.

Examples and conclusion

We exemplify the error bounds on the examples provided in [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF], together with additional numerically challenging examples that we custom made for illustration purposes.

The numerical behavior of the algorithm is illustrated on Test 1, given in the rst line of Table 3, for which 101 terms are computed. In Figure 1a, the magnitude of the coecients c 0 , . . . , c 100 is plotted on a log-scale. This is a higher precision 106-bits FP arithmetic computation, using the MPFR library [START_REF] Fousse | Mpfr: A multiple-precision binary oating-point library with correct rounding[END_REF], in order to accurately approximate their exact values. Their magnitude increases up to c 16 and then the convergent regime is observed. In Figure 1b the relative rounding error on each coecient is plotted, when the loop is evaluated with a 53-bit FP arithmetic (this rounding error is estimated by comparing with the shadow 106-bit higher precision computation). These errors are plotted in terms of the roundo unit u = 2 -53 . The corresponding evaluated sum s and probability P0:N are recalled in Figure 1d.

For comparison, we also tested an interval arithmetic implementation with a 53-bit precision interval format, using the MPFI library 2 . While these intervals provide enclosures of all the accumulated rounding errors, we observe in Figure 1c, where the radius of the intervals is plotted in terms of u, that they highly overestimate the actual rounding errors. This is conrmed in 1d: with interval arithmetic, the nal absolute enclosure radius is 2.9439e12 u ≃ 3.26e-4, while the shadowed absolute error bound is 2.5920e6 u≃ 2.87e-10. In Table 3, we computed both the fully rigorous relative error bound of Theorem 1 and its linearization (20) w.r.t. u. For easy examples requiring less than 50 terms, both MPFI and our bounds provide very sharp enclosures.

Both of our bounds are almost identical. A rapid increase of interval widths with MPFI is observed when N is larger than 100 in most of the cases, whereas our bounds continue to guarantee at least one correct digit in very hard cases requiring about N = 10 7 terms. However, when the ratio σ y /σ x becomes very small, like in Case Alfano 5 and Custom 8, the bounds provided are very loose and the rigorous bound deviates from its linearization. As mentioned also in [START_REF] Serra | Fast and accurate computation of orbital collision probability for short-term encounters[END_REF], such extremely degenerate cases are rarely occurring in practice (this roughly corresponds to the integration domain becoming uni-dimensional).

Therefore, we believe that the bound provided in this article can be of highly practical use. One can simply evaluate the provided closed-form bound and there is no need to overload the C code with shadowing computations or with an additional execution with interval arithmetic (which would encumber especially the on-board implementation). Our study also shows that the approach of [START_REF] Mezzarobba | Rounding error analysis of linear recurrences using generating series[END_REF] and the additional mathematical tools developed herein are applicable to an algorithm currently used in practice. While described by pen-and-paper, the formulas presented can be easily obtained by a computer algebra software (as shown in the jointed Maple code), so they can be at least partly automated. Future works, include possible renements concerning underow handling, an average case analysis in the framework of [START_REF] Higham | A new approach to probabilistic rounding error analysis[END_REF], a formal proof of these results and further generalizations to other implementations of similar mathematical functions like Chi-square densities (which have similar algebraic properties).

A Appendix

A.1 Complementary proofs for the rounding error analysis of the parameters of Algorithm 1

The following results complete the proof of Proposition 1. Note that we insist on bounding each P i , Q i in terms of the corresponding P ♯ i , Q ♯ i for simplicity reasons only, and slightly better bounds can be obtained without this restriction. Proposition 6. At step 2 of Algorithm 1, the computed value

Q2 of Q 2 satises |Q 2 -Q2 | ⩽ Q ♯ 2 γ 18 .
Proof. The value of Q 2 is computed as

Q2 = RN(RN(RN(RN( p2 ) • RN(RN(R 2 ) 2 )) • φ) • RN( φ + 2)).
By Equation ( 8), this implies

Q2 = (pR 2 (1 + θ 2 )) 2 φ( φ + 2)(1 + θ 4 ) = (p(1 + θ 2 )R 2 (1 + θ 2 )) 2 φ( φ + 2)(1 + θ 4 ) = p 2 R 4 φ( φ + 2)(1 + θ 12 ). It follows that Q2 -Q 2 = p 2 R 4 ((2 + φ + ϕ)( φ -ϕ) + φ( φ + 2)θ 12 )
and therefore

| Q2 -Q 2 | ⩽ p 2 R 4 (4γ 4 + 3γ 12 ) ⩽ Q ♯ 2 4 3 γ 4 + γ 12 ⩽ Q ♯ 2 γ 18 . Proposition 7. At step 2 of Algorithm 1, the computed value Q3 of Q 3 satises |Q 3 -Q3 | ⩽ Q ♯ 3 γ 24 .
Proof. One has with ψ dened in [START_REF] Thomassin | Asteria : Autonomous collision risks management[END_REF]. Then for all ξ ⩾ 0, the series Ψ(ξ) such that Ψ(λ) = L Ψ (λ

Q3 = p3 (R 2 (1 + θ 1 )) 3 φ2 (1 + θ 7 ) = (p(1 + θ 2 )) 3 (R 2 (1 + θ 1 )) 3 φ2 (1 + θ 7 ) = p 3 R 6 φ2 (1 + θ 16 ), hence Q3 -Q 3 = p 3 R 6 (( φ -ϕ)( φ + ϕ) + φ2 θ 16 ) and | Q3 -Q 3 | ⩽ p 3 R 6 [2γ 4 + γ 16 ] ⩽ p 3 R 6 • γ 24 = Q ♯ 3 • γ 24 .

3. 3 2 y.

 32 of the reference[START_REF] Higham | Accuracy and stability of numerical algorithms[END_REF] Chap. 3.4]. The case of the other parameters P i and Q i is similar; see Appendix A.1 for detailed proofs.For c 0 , denote z = -Then z = z(1 + θ 4 ), and c0 = R 2 2σxσy e z (1 + θ 2 ) (1 + θ 4 ) = c 0 e zθ 4 (1 + θ 6 ).

Proposition 3 (

 3 Given two series a(λ) = +∞ n=0 a n λ n and b(λ) = +∞ n=0 b n λ n , denote by a(λ) ≪ b(λ) the fact that |a n | ⩽ b n for all n ⩾ 0. In particular, this implies that the b n coecients are nonnegative real numbers. We denote by |a|(λ) = +∞ n=0 |a n |λ n the series of absolute values of coecients. Corollary of Proposition 2). The generating series of local errors satises this dierential inequality:

This is a bound

  on | δ|, in the Laplace plane. To go back to |δ| and obtain an inequality |δ|(λ) ≪ ∆(λ), consider the series ∆ and Ψ dened by ∆(λ) = λ -2 L ∆ (λ -1 ) and Ψ(λ) = L Ψ (λ -1 ) (with no λ -2 factor in the latter). Standard Laplace transform theory, see [17, Chap. 5, 8], gives ∆(ξ) = e 0 f (ξ) + γ (Ψ * f )(ξ), with(Ψ * f )(ξ)= ξ 0 Ψ(τ )f (ξ -τ)dτ the convolution of Ψ and f . A technical but straightforward computation 1 (see Lemma 3 in Appendix A.2)

  . Combining these two bounds yields |s -s| ⩽ (γ N -1 + ν + γ N -1 ν)f (1). Finally, the relative errors |τ ′ | ⩽ τ and |θ 1 | ⩽ u induced by the evaluation of exp(-pR 2 ) and the multiplication by s give | P0:N -P 0:N

  2 https://gitlab.inria.fr/mpfi/mpfi Coe. rel. err. (in terms of u) for 53-bit precision Interval radius (in terms of u)

Figure 1 :

 1 Figure 1: Loop evaluation results for Test 1.

Proposition 8 . 8 and|

 88 At step 3 of Algorithm 1, the computed value P0 of P 0 satises |P 0 -P0 | ⩽ P ♯ 0 γ 10 . Proof. One hasP0 = ((p( 1 2 φ + 1)(1 + θ 2 ) + ωx )(1 + θ 1 ) + ωy )R 2 (1 + θ 3 ) = {[p(1 + θ 2 )( 1 2 φ + 1)(1 + θ 2 ) + ω x (1 + θ 5 )](1 + θ 1 ) + ω y (1 + θ 5 )}R 2 (1 + θ 3 ) = {p( 1 2 φ + 1)(1 + θ 8 ) + ω x (1 + θ 9 ) + ω y (1 + θ 8 )}R 2 , hence P0 -P 0 = pR 2 ( 1 2 ( φ -ϕ) + ( 1 2 φ + 1)θ 8 ) + ω x R 2 θ 9 + ω y R 2 θ P0 -P 0 | ⩽ pR 2 ( 1 2 γ 4 + 3 2 γ 8 ) + (ω x + ω y )R 2 γ 9 ⩽ R 2 ( 3 2 p( 1 3 γ 4 + γ 8 ) + (ω x + ω y )γ 9 ) ⩽ P ♯ 0 γ 10 .Proposition 9. At step 4 of Algorithm 1, the computed valueP1 of P 1 satises |P 1 -P1 | ⩽ P ♯ 1 γ 18 . Proof. One has P1 = {[ 1 2 p φ( φ + 5)(1 + θ 3 ) + ωx ](1 + θ 1 ) + ωy (2 φ + 1)(1 + θ 2 )}p(R 2 (1 + θ 1 )) 2 (1 + θ 3 ) = {[ 1 2 p φ( φ + 5)(1 + θ 5 ) + ω x (1 + θ 5 )](1 + θ 1 ) + ω y (2 φ + 1)(1 + θ 7 )}pR 4 (1 + θ 7 ) = pR 4 ( 1 2 p φ( φ + 5)(1 + θ 13 ) + ω x (1 + θ 13 ) + ω y (2 φ + 1)(1 + θ 14 )) Proof. One has P3 = p3 ωy ((R 2 (1 + θ 1 )) 2 (1 + θ 1 )) 2 φ2 (1 + θ 7 ) = p 3 ω y ((R 2 (1 + θ 1 )) 2 (1 + θ 1 )) 2 φ2 (1 + θ 18 ) = p 3 ω y R 8 φ2 (1 + θ 24 ), hence P3 -P 3 = p 3 ω y R 8 ( φ2 -ϕ 2 + φ2 θ 24 ) and | P3 -P 3 | ⩽ p 3 ω y R 8 (2γ 4 + γ 24 ) ⩽ P ♯ 3 γ 32 .Proposition 12. At step 15 of Algorithm 1, the computed value ã of exp(-pR 2 ) satises |exp(-pR 2 ) -ã| ⩽ τ exp(-pR 2 ) with τ = exp(pR 2 γ 2 )(1 + γ 2 ) -1.Proof. Denote z = -pR 2 . Then z = z(1 + θ 2 ). This gives ã = exp(z)(1 + θ 2 ) = exp(z) exp(zθ 2 )(1 + θ 2 ).

  xi * invlaplace(psi_aux, lambda, xi): Psi := collect(expand(Psi), exp(p*R^2*xi)); # remove the negative term -2*omega[y]*R^2 Psi := Psi + 2*omega[y]*R^2; # use (exp(t)-1)/t << exp(t) with t=p*R^2*xi to bound Psi t

Table 1 :

 1 Absolute rounding error bounds for the parameters. Proposition 1. The absolute FP rounding error for the parameters appearing in lines 17 and the preconditionner exp(-pR 2 ) in line 15 in Algorithm 1 is bounded as indicated in Table

Table 2 :

 2 Main properties of real and Laplace plane series.

  , γ = γ 40 and

	C :=	7 96	p 3 ω x R 8 +	7 12	p +	1 2	ω x p 2 R 6 +	9 4	p
		+	5 4	ω x +	15 4	ω y pR 4 +	3 2	p + ω x + 3ω y R 2 .

[START_REF] Titolo | An abstract interpretation framework for the round-o error analysis of oating-point programs[END_REF]

Proposition 10. At step 5 of Algorithm 1, the computed value

Proof. One has

and take the inverse Laplace transform term-by-term, which yields

We conclude by using the inequality (e x -1)/x ≪ e x and dropping the nal negative term.

Proof of Lemma 2. The denominator Q(λ)-γQ ♯ (λ) factors as (1-

) with all β i ̸ = 0. We prove that

which implies 1

and nally the desired result.

Let β denote one of the β i . If |β| ⩽ pR 2 , then clearly (21) holds. Now suppose |β| > pR 2 . Since λ = β -1 is a root of Q(λ) -γQ ♯ (λ), we have Q(β -1 ) = γQ ♯ (β -1 ). We observe:

We thus obtain (1 -pR 2 /|β|) 3 ⩽ 7γ, and (21) holds. Lemma 4. For all k ⩾ 1, we have

Proof. Replace the value of p by p + in the proof of Lemma 3 to have Ψ + (ξ) ≪ W + (ξ) e p + R 2 ξ and C + = 1 0 W + (τ )dτ . Also bound f (ξ) by e p + R 2 g(ξ) since p + ⩾ p. Then we have

For k = 1, we have as in the proof of Proposition 4

Repeating this process k times to push the exponential e p + R 2 ξ out of the convolution product gives W + e p + R 2 ξ * k * e p + R 2 ξ g (ξ) = e p + R 2 ξ (W + * k * g)(ξ).

Evaluating the series at ξ = 1 gives:


Pierre Le Bras
email: pierre.lebras@univ-brest.fr

Florian Sévellec

Pierre Tandeo

Juan Ruiz

Pierre Ailliot

Selecting and weighting dynamical models using data-driven approaches

In geosciences, multi-model ensembles are helpful to explore the robustness of a range of results. To obtain a synthetic and improved representation of the studied dynamic system, the models are usually weighted. The simplest method, namely the model democracy, gives equal weights to all models, while more advanced approaches base weights on agreement with available observations. Here, we focus on determining weights for various versions of an idealized model of Atlantic Meridional Overturning Circulation. This is done by assessing their performance against synthetic observations (generated from one of the model versions) within a data assimilation framework using EnKF. In contrast to traditional data assimilation, we implement data-driven forecasts using the analog method based on catalogs of short-term trajectories. This approach allows us to efficiently emulate the model's dynamics while keeping computational costs low. For each model version, we compute a local performance metric, known as the contextual model evidence, to compare observations and model forecasts. This metric, based on the innovation likelihood, is sensitive to differences in model dynamics and considers forecast and observation uncertainties.

Finally, the weights are calculated using both model performance and model codependency, and then evaluated on climatologies of long-term simulations. Results show good performance in identifying numerical simulations that best replicate observed short-term variations. Additionally, it outperforms benchmark approaches such as model democracy or climatologies-based strategies when reconstructing missing distributions. These findings encourage the application of the proposed methodology to more complex datasets in the future, like climate simulations.

Introduction

In the geosciences, several numerical models are usually available to represent the same complex system. For example, the Earth's global climate system is implemented in a range of different numerical models that are, for some, gathered within the Couple Model Intercomparison Projet (CMIP, Eyring et al. 2016). Due to the different parameterizations (related to unresolved

processes and numerical implementations) and given particular structural equations, a model represents an imperfect image of the studied system. In other words, each model is characterized by its own strengths and weaknesses, which impact its predictive skills. The use of the multi-model ensemble (MME) to predict a specific climate variable is then considered more informative than any individual model. Here, the models are assumed to complement each other using MME which has the potential to improve climate system representation [START_REF] Abramowitz | ESD Reviews: Model dependence in multi-model climate ensembles: weighting, sub-selection and out-of-sample testing[END_REF].

Nevertheless it remains the question of how to gather the information produces by the range of models. In this context it is usual to follow a model democracy approach, giving equal weight to all models to enhance the representativeness. This is what is mostly done in the 6 th Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC 2021). However, this strategy has been debated in the community due to the strong assumptions it implies: equal probability of all models to represent all variables in the whole phase-space of the system and independence between them (e.g., [START_REF] Knutti | The end of model democracy?[END_REF][START_REF] Sanderson | A representative democracy to reduce interdependency in a multimodel ensemble[END_REF][START_REF] Knutti | Uncertainty quantification using multiple models-Prospects and challenges[END_REF].

A trade-off between MME-based model democracy and selecting a single model would be to weight-average the models. The weights can be found through model performance, which consists in evaluating the consistency between model simulations and available observations. This procedure should be adapted to the specificities of model simulations being processed [START_REF] Eyring | Taking climate model evaluation to the next level[END_REF]. Hence, the dynamic behaviour of the system needs to be properly taken into account. This can be done, for example, by assessing the models' ability to describe successive sequences of observations. In addition, the weights should account for the degree of dependency between the models, as they usually share some elements in their structure or parameterizations [START_REF] Knutti | Challenges in combining projections from multiple climate models[END_REF][START_REF] Abramowitz | ESD Reviews: Model dependence in multi-model climate ensembles: weighting, sub-selection and out-of-sample testing[END_REF]). Several approaches have been suggested to measure the model-observations consistency and thus obtain the individual model weights.

The first approaches are based on descriptive statistics. This is the case of "reliability ensemble averaging" [START_REF] Giorgi | Calculation of average, uncertainty range, and reliability of regional climate changes from AOGCM simulations via the "reliability ensemble averaging"(REA) method[END_REF] and the "ClimWIP method" (Climate model Weighting by Independence and Performance, [START_REF] Sanderson | A representative democracy to reduce interdependency in a multimodel ensemble[END_REF] and [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF]. In those two methods, weights are based on the ability of model simulations to reproduce observed statistical diagnostics for one or more variables (e.g., mean, variability, or trend fields), while including a component that measures inter-model dependence. However, the score formalism does not evaluate how the models reproduce sequences of observations (which reflect model dynamics). ClimWIP method was applied in numerous studies at global or regional spatialscales with a range of climate variables such as atmospheric temperatures or precipitations [START_REF] Sanderson | A representative democracy to reduce interdependency in a multimodel ensemble[END_REF][START_REF] Sanderson | Skill and independence weighting for multi-model assessments[END_REF][START_REF] Brunner | Quantifying uncertainty in European climate projections using combined performance-independence weighting[END_REF][START_REF] Brunner | Reduced global warming from CMIP6 projections when weighting models by performance and independence[END_REF][START_REF] Merrifield | An investigation of weighting schemes suitable for incorporating large ensembles into multi-model ensembles[END_REF], arctic sea-ice [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF], or antarctic ozone depletion [START_REF] Amos | Projecting ozone hole recovery using an ensemble of chemistry-climate models weighted by model performance and independence[END_REF].

Second approaches are included in a probabilistic Bayesian framework, such as "Bayesian model averaging" (BMA, e.g., [START_REF] Raftery | Using Bayesian model averaging to calibrate forecast ensembles[END_REF][START_REF] Min | Probabilistic climate change predictions applying Bayesian model averaging[END_REF][START_REF] Sexton | Multivariate probabilistic projections using imperfect climate models part I: outline of methodology[END_REF][START_REF] Olson | A simple method for Bayesian model averaging of regional climate model projections: Application to southeast Australian temperatures[END_REF][START_REF] Olson | Accounting for skill in trend, variability, and autocorrelation facilitates better multi-model projections: Application to the AMOC and temperature time series[END_REF] and "kriging methods" applied at global and local scale [START_REF] Ribes | Making climate projections conditional on historical observations[END_REF]Ribes et al. 2022, respectively). In BMA, the distributions of model outputs are updated using available observations to obtain a posterior weighted distribution. The weights for candidate models are determined using Bayes' theorem which combines prior beliefs (typically uniform weights) with the global likelihood of the observations estimated for each model. These weights reflect the confidence in each model considering both the data and uncertainty. However, evaluating model consistency over time is tricky because the distributions involve time aggregated quantities on large time periods, making it less suitable for exploiting information about model dynamics.

In this study, we evaluate the performance of competing models based on their ability to reproduce short-term dynamics of the system described by observed sequences. To achieve this, a data-driven approach is adopted within a data assimilation framework [START_REF] Ruiz | Analog Data Assimilation for the Selection of Suitable General Circulation Models[END_REF]. One of the main advantage is the use of already existing simulations to produce free-model probabilistic forecasts (i.e. without having to rerun the models, which greatly reduces the computational cost), while retaining the useful properties of data assimilation [START_REF] Lguensat | The analog data assimilation[END_REF]. The data-driven forecasts are estimated using regression based on the classic least squares, but applied locally in the phase space, to capture the nonlinearities in the evolution of the system [START_REF] Platzer | Using local dynamics to explain analog forecasting of chaotic systems[END_REF]. The proposed methodology is relevant for various reasons. Firstly, it is written in a Bayesian framework (as BMA and kriging strategies), where data-driven forecasts are considered as Gaussian a priori information which are sequentially updated thanks to observations (which are also considered as Gaussian). This provides accurate and reliable initialization for the next forecast. The difference between the data-driven forecasts and the observations allows the computation of the innovation likelihood (e.g., [START_REF] Carrassi | Estimating model evidence using data assimilation[END_REF], which is used locally in time as metric of model performance. Here also, the scores take into account the model dependency as in the formalism of ClimWIP [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF]. We then obtain weights that we test by applying them to long-term simulation climatologies.

This study aims to develop a weighting methodology. As a first step, we use numerical simulations of an idealized chaotic deterministic model representing the centennial-to millennial dynamic of Atlantic Meridional Overturning Circulation (AMOC).

The proposed alternative approach is compared to more classic approaches, such as model democracy, climatology, or single best model.

The study is organized as follows. Section 2 details the methodological framework for obtaining model weights and applying them. Section 3 describes the numerical model and the data used to implement the methodology. Section 4 is devoted to the numerical results and Sect. 5 concludes this study and presents some perspectives for future research.

Methodology

This section first describes the framework to measure the ability of a single dynamical model to fit a set of noisy observations. After applying it to an ensemble of competing models, the second part is devoted to the strategies for computing model weights and their application.

Evaluating the local performance of a dynamical model

Here, we evaluate the model performance based on its short-term dynamics. By synchronising model forecasts to available observations, Bayesian data assimilation (DA) is a suitable framework for addressing this. DA, like Kalman filtering strategies, is commonly used to improve the estimation of the latent unknown state of a system by sequentially including the available 3 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License. observations [START_REF] Carrassi | Data assimilation in the geosciences: An overview of methods, issues, and perspectives[END_REF]. The associated state-space model is expressed as

x k = M (x k-1 ) + η k , ( 1a
)
y k = H (x k ) + ϵ k , ( 1b
)
where k is the time index, x k is the true state, M is the model propagator, η k is the model error, y k is the noisy observations, H is the observational operator, and ϵ k is the observation error. Equation (1a) represents the model relation between the true state and the previous state given an additive model error. Equation (1b) makes the relationship between the observations to the true state with an additive observational error.

In the classic Kalman filter, an assimilation cycle aims to sequentially update a Gaussian forecast distribution (generated by the model equations) with the information provided by the observations (also Gaussian), when available. The posterior analysis distribution is a more accurate and reliable estimation of the latent state by satisfying robust statistical properties (i.e., best linear unbiased estimator). The ensemble Kalman filter (EnKF) is preferred when dealing with nonlinear systems [START_REF] Evensen | The ensemble Kalman filter: Theoretical formulation and practical implementation[END_REF]. At each time step, the two first moments of the Gaussian forecast distribution are approximated empirically using a Monte Carlo approach. To do this, N forecast members are generated from the same model equations, each one initialized from a sample of the previous analysis distribution. This is expressed as

x f (j),k+1 = M x a (j),k + η (j),k+1 (2)
where j is the member sample index, x f (j),k+1 is the j th forecast member at time k + 1, x a (j),k is the j th analysis member at the previous time k and η (j),k+1 is the model error assumed to be Gaussian. In Eq. (2), x a (j),k is calculated as

x a (j),k = x f (j),k + K k y k -H(x f (j),k ) + ϵ (j),k (3)
where ϵ (j),k is the observation error drawn from a random sample from the multivariate Gaussian with the 0-mean vector and an error covariance matrix R. The Kalman gain K k in Eq. ( 3) is defined as

K k = P f k H T HP f k H T + R -1 , (4)
where P f k is the empirical covariance matrix of the forecast distribution at time k calculated using the EnKF j members, H represents the tangent linear of the observation operator H.

DA allows the model performance evaluation by measuring the consistency between the model forecasts and the available observations. Directly computed in a DA cycle, the Contextual Model Evidence (CME) is defined as the log-likelihood of observations at a given time for a model, taking into account the forecast state as prior information [START_REF] Carrassi | Estimating model evidence using data assimilation[END_REF]).

Larger CME is obtained for more consistency between the model forecast and observations. Under Gaussian assumption included in the EnKF cycle, CME is approximated as follows where x f k is the empirical mean of the ensemble forecasts and r is the number of available observations. While RMSE is a L 2 metric between the prior forecast state and the observation, similarly the CME can be seen as a distance between the same two fields [START_REF] Metref | Estimating model evidence using ensemble-based data assimilation with localization -The model selection problem[END_REF]. The advantage of CME is that it takes into account reliability information from both forecasts and observations (via covariance matrices) in addition to accuracy information (i.e., the difference between mean states).

CME k (y; M) = - 1 2 y k -H(x f k ) HP f k H T + R -1 y k -H(x f k ) T - 1 2 ln det HP f k H T + R - r 2 ln(2π), (5
In principle, DA formalism requires access to the numerical model to propagate the state of the system in order to obtain the forecasts. When performing EnKF, a large number of model forecasts are generated, representing a significant computational cost. An efficient and flexible data-driven alternative aims at combining DA with a statistical forecasting strategy based on analogs [START_REF] Lorenz | Atmospheric predictability as revealed by naturally occurring analogues[END_REF]. In Analog Data Assimilation (AnDA), existing long-term simulations are used instead of performing new sequential models simulations [START_REF] Tandeo | Combining analog method and ensemble data assimilation: application to the Lorenz-63 chaotic system[END_REF][START_REF] Lguensat | The analog data assimilation[END_REF]. The simulations are decomposed as a catalog of short-term trajectories (i.e., pairs of analog states and their successors) sampling the short-term model dynamics. Hence, using an appropriate distance, the analog states closest to the assimilated state are selected from the catalog. Their successors are combined using a statistical operator to obtain probabilistic forecasts. Hence, AnDA corresponds to perform a classic DA process, but replacing in Eq. ( 1a) the model propagator by its analog-based approximation. This reads:

x k = M (x k-1 ) + η k , (6)
where M refers to the analog-based forecasting propagator.

In practice, it consists in fitting a local linear regression [START_REF] Cleveland | Locally weighted regression: an approach to regression analysis by local fitting[END_REF] relating the selected analogs to their successors, which has proven to be robust [START_REF] Zhen | An adaptive optimal interpolation based on analog forecasting: application to SSH in the Gulf of Mexico[END_REF][START_REF] Platzer | Using local dynamics to explain analog forecasting of chaotic systems[END_REF]). The regression is able to learn the timelocal dynamic of the model including nonlinearities since the linear adjustement is applied locally in the state space. By projecting the current state with the regression, statistical approximated forecasts are obtained. The forecast takes into account the regression uncertainty (assumed Gaussian) which is a robust estimate of the model error. When the catalog is large enough (to approximately cover the full state space), it has been shown that AnDA performs as well as DA [START_REF] Lguensat | The analog data assimilation[END_REF][START_REF] Ruiz | Analog Data Assimilation for the Selection of Suitable General Circulation Models[END_REF] have recently shown, using different idealized dynamical models, that CME can be robustly estimated using AnDA.

The combination of AnDA with CME is the basis of the current study.

Strategies for model weighting-average

The methodology described above is applied to an ensemble of models in order to individually measure their forecast skills given a set of observations. For each model, CME time series are obtained and used to derive individual weights. The weights are then applied to provide a model-average representative of the observations.

Calculating CME-based weights

The CME time series can be processed in various ways to obtain a single scalar (i.e., the score) used to derive the model weight. In this study, three CME scores are defined, together with three benchmark scores (not using the CME). Hence, we (i) . For simplicity, the scores are here described prior to their normalization, which leads to the model weights.

(a) The model democracy [START_REF] Knutti | The end of model democracy?[END_REF] gives the same weight to all models, regardless of model performance. It reads:

w (i) democracy = 1. (7a) (b)
The climatological score is based on the comparison of model and observation distributions. It measures the minimum cumulative value of the two histograms, highlighting the common area between both distributions. A score close to 0 denotes a model that poorly simulates the observed distribution, while a score of 1 is obtained by a model with a perfect climatological distribution [START_REF] Perkins | Evaluation of the AR4 climate models' simulated daily maximum temperature, minimum temperature, and precipitation over Australia using probability density functions[END_REF]). It reads:

w (i) climato = n j=1 minimum(Y (j), M (i) (j)), (7b)
where Y and M (i) are the normalized histograms (such that their respective sum equals 1) of the observations and the model simulations, respectively, j is the index of histogram equal-width bins, and n is their total number. Note that in the current study, single-variable observations are used but the score could be adapted using multivariate distributions.

(c) The single best model score assigns 1 to the model showing the best climatological score and 0 to the others. Opposite to the democracy score, it only takes into account the best-performing model over the whole set of observations [START_REF] Tebaldi | The use of the multi-model ensemble in probabilistic climate projections[END_REF]. It reads:

w (i) single =      1, if w (i) climato > w (j) climato ∀i ̸ = j; 0, otherwise. (7c) (d)
The CME-ClimWIP score is based on the ClimWIP approach [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF]. Our adaptation reads:

w (i) CME-ClimWIP = e D i σ D 1 + L j̸ =i e S ij σ S , (7d)
where D i measures the intra-model performance with its associated shape parameter σ D and S ij measures the intermodel dependency with its associated shape parameter σ S . D i and S ij are usually expressed with Euclidian distance [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF]). Here, we evaluate D i such as

D i = K k=0 CME k (y; M (i) )
, where K is the number of time step and S ij as the sum of the defined inter-model such that

S ij = K k=0 CME k (M (i) ; M (j) ) j̸ =i
w (i) punctual = 1 K K k=0 1 {CME k (y;M (i) )>CME k (y;M (j) )} i̸ =j , ( 7e
)
where 1 is the indicator function assigning 1 when the condition {CME k (y;

M (i) ) > CME k (y; M (j)
)} i̸ =j is satisfied and 0 otherwise. This score captures the good performance of model in specific regions of the phase space, despite not necessary being optimal in other regions. For example, extreme values are often only well captured by a few specific models. By retaining the information from only one model at each time k, the score bypasses the need to identify temporal inter-model similarities.

(f) The CME best persistent model score w

(i)
persistent is derived from the previous score by overweighting models that are the best on several consecutive states. The score is linearly proportional to the cumulative number of consecutive time steps where the model is the best (expression not shown here). This score emphasises consistently good local dynamics.

Assessing the performance of weighting-averages

Here, the methodology consists in two successive stages (Fig. 1). The first step corresponds to a training stage where AnDA is applied to compute the CME time series for each model. The second step is a testing stage in which the CME are used, to compute the weights (i.e, normalized scores) as described in the previous subsection. Applying the weights to individual model distributions leads to a reconstructed distribution. To test the skill of this reconstructed distribution, its overlap with the true distribution is computed (in the same way as the formalism of the climatological score w (i) climato , see Eq. ( 7b)). Comparison between the overlaps obtained using different scores allows us to compare and then rank the weighting strategies.

3 Experimental setup

Idealized AMOC model

The study is based on an idealized autonomous low-order deterministic model of the AMOC able of reproducing its millennial variability within a chaotic dynamics [START_REF] Sévellec | Millennial variability in an idealized ocean model: Predicting the AMOC regime shifts[END_REF]Sévellec and Fedorov 2015). The model derives from the long history of salinity loop-models (originally proposed by [START_REF] Welander | Note on the self-sustained oscillations of a simple thermal system[END_REF]Welander 1965;Welander 1967). Its formalism with three equations [START_REF] Dewar | Fluid flow in loops driven by freshwater and heat fluxes[END_REF][START_REF] Huang | Haline circulation: Bifurcation and chaos[END_REF][START_REF] Huang | Haline circulation: Bifurcation and chaos[END_REF] reads: where t is the time, ω is the variable component of the overturning circulation strength, S BT and S N S are the bottomtop salinity gradients and the north-south salinity gradients, respectively, Ω 0 is the steady part of the overturning circulation strength which acknowledges the impact of constant temperature and surface winds, ϵ is the buoyancy torque, λ is a linear friction, β is the haline contraction coefficient, K is a linear damping and F 0 S 0 /h is the surface salt flux (where F 0 is the freshwater flux intensity, S 0 is a reference salinity, and h is the loop thickness i.e., the depth of the level of no motion).

dω dt = -λω -ϵβS N S , (8a)
dS BT dt = (Ω 0 + ω) S N S -KS BT + F 0 S 0 h , (8b)
dS N S dt = -(Ω 0 + ω) S N S -KS N S , (8c)
Equation (8a) describes the momentum balance. Here, the rate of change of the ocean overturning strength is driven by the meridional salinity gradients accounting for the buoyancy torque and a linear friction. Equations ( 8b) and (8c) describe the evolution of the salinity gradients driven by advection, a linear damping representative of diffusion, and the surface salt flux.

Numerical integrations are done using a 4 th order Runge-Kutta method with a time step of 1 year. All reference parameters are set following Sévellec and Fedorov 2014 (which we refer interested reader to). Their values are h = 1000 m, Ω 0 = -2.5 × 10 -2 yr -1 , F 0 = 1m yr -1 , S 0 = 35 psu, λ = 10 -2 yr -1 , ϵ = 0.35 yr -2 , K = 10 -4 yr -1 and β = 7 × 10 -4 psu -1 .

Ensemble of AMOC model versions

Various versions of the AMOC model are obtained by perturbing some parameters in the equations: Ω 0 , K and λ ϵ , while other remain fixed. All possible combinations between the reference values of the parameters, twice of them, and half of them lead to twenty-seven distinct parameterizations. Eleven are selected by ignoring the versions describing singular dynamics (e.g., with strong instabilities, constant over time, or presenting periodic behaviour). Thus, for the same equations, the 11 sets 215 of parameters induces their own dynamics. These dynamics show similarities and differences that can be visualised by the computation of the trajectories in the phase space (Fig. 2).

Synthetic data

A model-as-truth experiment framework is set up (e.g., [START_REF] Herger | Selecting a climate model subset to optimise key ensemble properties[END_REF], in which one model is used to generate noisy pseudo-observations. Here, only ω is observed and will be used as the variable of interest. To efficiently select the analogs, the three variables are normalized (here rescaled). In particular, the normalized values of the original variable ω are associated with the new variable called "normalized AMOC" (Fig. 3). In all the experiments, the covariance of the observations for the normalized AMOC is fixed to R = 0.5I 3 . The AnDA framework is applied over a period of K=8000 years with an assimilation step of 20 years, thus resulting in 400 assimilation cycles. In each AnDA cycle, the forecast distributions are estimated by selecting 10000 analogs (which avoids computational divergence for certain attractors with highly nonlinear dynamics in the set) from the associated catalog of size 400 000 years, for each of the 200 EnKF-members. The choice of the assimilation step allows the forecasts to be sufficiently differentiated, so that the competing models can be distinguished locally over time when using the CME (sensitivity experiments not shown). The different methods are evaluated using two protocols. In the first one, we consider a "perfect experiment", where the model used to generate the pseudo-observations is also included in the catalog. In the second one, the model is excluded from the catalog and we talk about "imperfect experiment". For both perfect and imperfect cases, we first describe an illustrative experiment when model 8 is used as pseudo-observations, and then summarized the results for each model alternately used to generate the pseudo-observations following a Leave-One-Out experiments design.

Perfect model experiment

The goal of the perfect model approach is to measure the ability of the three CME-based scores to retrieve the correct model (by giving it a predominant weight) from a pool of catalogs, including the true one. The skills of these scores are then compared with the skills of the three benchmarks. In all the perfect model experiments of this study, the CME-ClimWIP score is computed using the optimal values of σ D =0.3 and σ S =0.4 which maximize the overall overlap scores in all the experiments combined (following [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF][START_REF] Lorenz | Prospects and caveats of weighting climate models for summer maximum temperature projections over North America[END_REF].

Experiment with model 8

AnDA is performed on the eleven models with the same pseudo-observations from model 8 and the CME is computed for each assimilation cycle. The resulting CME series for each individual model varies significantly over time and differs locally from model to model (Fig. 4). CME of models 0, 1, 2, 6, and 10 is more than 75% lower than -1 (Fig. 5), showing their inefficiency to replicate model 8 system dynamics. The CME distributions of models 3, 7, and 9 show a narrower distribution with 75% values greater than -1. This means they are both accurate and reliable. As expected, the model with the highest and more consistent CME is model 8. This demonstrates the ability of the CME to retrieve the correct model from a range of models.

The climatologies-based score produces almost uniform weights close to model democracy (Table 1). Indeed, the model climatological distributions are very comparable. This impacts the reliability of selecting a single model following the best model score. Beyond that, the observation error leads the climatologies-based score to misidentify model 9 instead of model 8. The CME-ClimWIP score assigns its highest weight to the correct model 8, although models 7 and 9 have close weights.

Regarding the two other CME-based scores, although the weight associated to model 8 is high, these scores give the highest weight to model 9. We conclude that the time selection characteristics of them make both scores more sensitive to observation errors.

Overall eleven experiments

By varying the model used to generate the pseudo-observations in the eleven AnDA experiments, we can robustly assess the ability of the six scores to efficiently retrieve the correct model (Table 2). As emphasized in the specific experiment for model 8, the scores based on climatology do not allow a clear discrimination between models. Hence the climatological and the bestsingle-model only retrieved the correct model twice over the eleven experiments (i.e., experiments with models 1 and 10). The CME-based scores perform better, attributing their highest weight to the correct model 8 times for CME-ClimWIP and both 6 times for CME best punctual and CME best persistent. As for model 8 experiment, in most of the unsuccessful experiments of the three CME scores, the correct model still has a weight close to the weight of the selected one (not shown here). 2. Summarised perfect model results of the Leave-One-Out model-as-truth experiments. For each column representing an experiment (i.e., the model index used to generate pseudo-observations), the index of the model with the highest weight is specified for each score in the row. The indices in bold show when the true model is recovered. The last column summarizes how frequently the correct model is identified among the candidate models across the 11 experiments for the five scores. Note that 'N/A' is indicated for model democracy since the score prevents differentiation between models.

Imperfect model experiment

The imperfect model approach aims to reconstruct the statistical properties of the distribution of a missing model using distinct models. In all the imperfect experiments of this study, the CME-ClimWIP score is computed using the optimal values of σ D = 0.5 and σ S = 0.7 which maximize the overlap scores in all the experiments combined (following [START_REF] Knutti | A climate model projection weighting scheme accounting for performance and interdependence[END_REF][START_REF] Lorenz | Prospects and caveats of weighting climate models for summer maximum temperature projections over North America[END_REF].

Experiment with model 8

The model democracy approach yields the worst reconstruction (Fig. 6, left panel), stressing that, in this case, the observations information should be included in the weight calculation. skewed, underestimating the true normalized AMOC (Fig. 6, left panel). The three CME-based scores outperform the three benchmark ones, with greater overlaps of their reconstruction with the truth. In particular, they all reconstruct better statistics for high positive values of normalized AMOC. For CME-ClimWIP, three models contribute the most to the reconstruction of the true missing model: model 3 by 18.1%, model 7 by 49.1%, and model 9 by 31.7%. To a lesser extent, the same three models are also selected for the best punctual and best persistent scores, but still leaves a greater contribution from other seven models. This suggests that the information provided by these seven models is negligible, as the two last CME strategies have lower overlapping scores than the one obtained using CME-ClimWIP.

Overall eleven experiments

When the eleven models are reconstructed independently in the imperfect framework, the three CME-based scores produce better reconstructions on average than the democracy and climatologies scores (Fig. 7). Similar to the specific model 8 reconstruction, the climatologies comparison score does not provide meaningful distinctions between the models. This results in high variations in the reconstruction score when using the best single model across all experiments, further highlighting its unreliability. For instance, while the score effectively captures most of the distribution when reconstructing models 1 and 2 (due to their highly similar dynamics), for other experiments, it yields lower reconstruction performance.

The reconstruction performances in the eleven experiments show that CME-ClimWIP greatly outperforms the model democracy for 7 of 11 model reconstructions. As for model 8 reconstruction, CME-ClimWIP gives higher weights to a small subset of models, making it well suited to reconstruct distributions that share similarities with some others in the ensemble.

On the other hand, CME-ClimWIP is close but less useful than democracy (and climatologies) for reconstructing models 3, 5, 9, and 10. Here, despite the six scores having lost performance, uniform weights are appropriate for good reconstruction.

Relative influence of the CME-ClimWIP score components

As expressed in (7d), the CME-ClimWIP score is composed of the performance and dependence of the models, both expressed using the CME. The CME at time k depends on the innovation (i.e., the difference between x f k and y k ) and the covariance error matrices P f k and R. The degree of contribution of the performance and independence is controlled by two tuned shape parameters: σ D and σ S . Hence, to assess the importance of each component in the reconstruction performance, various versions of the CME-ClimWIP score are tested, by taking into account or not the influence of each component (Table 3).

The degree of dependence between models does not contribute significantly to the reconstruction performance (Table 3).

Since the eleven model versions are based on the same set of equations and the three variables are normalized, the modelmodel redundancy is significant. This implies that the inter-model dependency (i.e., the denominator in Eq. (7d)), is too similar across the models to contribute notably to differentiate them. The shape parameter (σ D ) has the main impact on the CME-ClimWIP by driving the score performance. The reconstruction performance of the CME-ClimWIP version without tuning the shape parameters is more than 10% lower to the reference performance. The score computed without information from P f and R has lower performance than the reference, confirming the usefulness of CME (including the covariances in its formulation) over a metric only based on Euclidian distance. Finally, when both covariances and shape parameters are excluded from the score, the performance lost is even more important. 3. Sensitivity of the reconstruction performance (based on overlap of the reconstructed distribution with the truth) to various alternative/simplified versions of the CME-ClimWIP score in the imperfect model approach. The second column shows the overlap associated with model 8 experiment. The third column shows the results of all Leave-One-Out experiments combined

Conclusion and perspectives

This study aims at developing a data-driven methodology for weighting models, only based on their ability to represent the dynamics of observations. For this purpose, a set of dynamical models is compared to noisy observations, where model equations are replaced by available simulations. The method combines a machine learning approach (i.e., an analog forecasting method) to estimate the model forecasts in a cost-effective manner (i.e., without the need to re-run the model) and a sequential data assimilation algorithm (i.e., the ensemble Kalman filter). The time-varying performance of the models with respect to observations is evaluated using the contextual model evidence (also known as innovation likelihood), which benefits from DA properties.

To test this methodology, an ensemble of eleven models is obtained by perturbing parameters of an idealized chaotic model of the AMOC. For each model version, the equations are only used to generate large simulations of its three variables. Each version alternately plays the role of truth (i.e., Leave-One-Out experiments) used to construct pseudo-observations of the AMOC component. In the eleven experiments, model weights are extracted from the CME series using different strategies.

They take into account the performance of the models with respect to pseudo-observations and the degree of similarity with other models. The method is then assessed by applying the weights to the distributions of long-term model simulations. In this way, we test the extent to which the short-term dynamics of individual models can provide relevant information for reconstructing the statistics of a targeted distribution. Reconstruction performance is measured by the percentage of overlap between the reconstructed and the true distributions. The reconstruction performance associated with the three CME-based scores is compared to three benchmark approaches that do not include DA framework (i.e., model democracy, climatological distributions comparison, and best single model). The results of the perfect model approach highlight a better performance of CME-based scores in recovering the correct model, compared to benchmarks, suggesting the importance of using DA. The results of the benchmark strategies generally suffer from a lack of discrimination between models to correctly identify the right one. In the context of the imperfect models approach, the scores based on the CME are able to reconstruct the targeted distribution more suitably than the benchmarks using the same partial noisy observations. The results underline that CME-based scores can be seen as a compromise between the "democracy score" and the so-called "dictatorial score" which selects only a single model. Within the CME-based score, the CME-ClimWIP score is relatively closer to "dictatorship". It is suitable for reconstructing distributions sharing similarities with few models in the ensemble and when democracy is not appropriate. On the other hand, the CME best punctual and CME persistent scores are closer to "democracy". By exploiting temporal performance, their weights are more adaptive, which is advantageous for outperforming the already successful democracy in any experiment (when CME-ClimWIP does not succeed).

This study initiates the implementation of the methodological framework whose ultimate goal is to weight a larger ensemble of dynamic models representing more complex systems, such as those provided by the CMIP6 [START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF]. Here, the aim will be to find weights that will be calculated on the performance of climate models (represented by their catalog of simulations) in correctly reproducing current partial observations. Extrapolating the weighted model averages using future climate projection simulations could provide a more accurate and reliable response of the AMOC to anthropogenic forcing. However, some aspects need to be examined due to the differences between a CMIP study and the present one. Firstly, whereas here the simulations are stationary, the CMIP simulations describe a trend and variabilities in response to the time-varying forcing (e.g., solar activity, volcano eruption, and greenhouse gas emissions). Secondly, in this proof of concept study, the full phase space of the idealized model is known, whereas using CMIP simulations, additional variables would need to be investigated to increase the unknown number of dimensions. This could be valuable for obtaining accurate and reliable AnDA reconstructions needed for deriving relevant weights. The diversity of available observations types could be particularly beneficial in this case.

  doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

  as the score assigned to model M

  . The CME-ClimWIP score is higher if model forecasts match observations and if they are sufficiently different from other-model forecasts. The two parts of Eq. (7d) are balanced thanks to the use of the parameters σ D and σ S . When they are appropriately set within the score, a few top-performing models are selected. Notably, higher values for these parameters bring the score closer to the model democracy while when σ D =1 and σ S tends to ∞ the BMA expression is retrieved (see Appendix A for more details on the link between the CME-ClimWIP score and the BMA framework). Note that σ D and σ S are determined according to the specific experiments model tests following Knutti et al. 2017 and Lorenz et al. 2018. 6 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License. (e) The CME best punctual model score exploits the time-local performance provided by the CME. It assigns a local-value of 1 to the best model at time k and 0 to the other models. The score reads:

Figure 1 .

 1 Figure 1. Schematic of the two-stages methodology for one experiment using 11 models. Left panel: training stage performing AnDA with the 11 catalogs and a single set of observations leading to CME computation. Bottom center: from the CME-based score and benchmark score model weights are obtained (blue and red, respectively). Right panel: by applying the weights to individual model distributions (grey), a single reconstructed distribution associated with each score is obtained (blue an red) which is compared to the truth (dashed black). The quality of the reconstruction is assessed using a distribution overlap (blue and red shadings).

8

  https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 2 .

 2 Figure 2. Chaotic trajectories in the phase space of the eleven perturbed versions of the 3-variable (normalized here) AMOC model (Sévellec and Fedorov 2014). The values of the perturbed parameters λ, ϵ, Ω0 and K used to generate the eleven model versions are indicated in the title. The attractor labelled 0 (top left) corresponds to the model with the reference set of parameters.

Figure 3 .

 3 Figure 3. Distributions of normalized AMOC for the eleven model versions

  10 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

  11 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 4 .

 4 Figure 4. AnDA results on 400 cycles using model 8 to generate the pseudo-observations: examples of three assimilated models among the eleven. Top panel: time series of normalized AMOC with forecasts mean states and 90% prediction interval (in shades). The red series denotes the true model (i.e., model 8), the orange and blue ones refer to models 0 and 7, respectively. The black dots represent the noisy observations generated from the true states denoted by the dashed black line. Bottom panel: associated CME time series for each model. Close to zero values of the CME indicate better forecast distribution.

Figure 5 .

 5 Figure 5. CME distributions from AnDA results using model 8 to generate the pseudo-observations. CME distribution are on the y-axis and model number (from index 0 to 10) on the x-axis. The boxes span the 25 th , median, and 75 th percentiles of the distributions, while the whiskers show the minimum and maximum values (excluding the outliers denoted with the black diamonds). Mean is indicated with the small yellow triangle.

  The climatology-based score slighty improves the reconstruction provided by the model democracy (+0.8% compared to model democracy). By selecting model 9, the best single-model score improve by 6% the reconstruction of the true model 8 distribution, when compared to the model democracy. The reconstructed distributions associated with the three strategies (i.e., model democracy, climatology, best single model) are all positively 13 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

Figure 6 .

 6 Figure 6. Imperfect model results of model 8 experiment (by excluding it). Left panel: reconstructed distribution of normalized circulation strength for the six scores with different colors. They are compared with the true distribution in dash black. The overlap score (in %) for each score is expressed in the legend in brackets. Right panel: table of model weights (in %) for each score (the sum of the weights per row is 100). Values in bold highlight the weights greater than the model democracy weight (>10%) which correspond to the models contributing the most to the distribution reconstruction.

A

  typical example of such case is the reconstruction of model 10 whose distribution is quite symmetric, whereas all the other 14 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License. models have asymmetric distribution. For model 3, 5, 9, and 10 experiments, CME best punctual and persistent give better results than CME-ClimWIP. It is worth noting that these two scores also have greater reliability across all eleven experiments than CME-ClimWIP, since their reconstruction performance remains consistently superior to the model democracy score with minimal variation between all the experiments.

Figure 7 .

 7 Figure 7. Imperfect model results (i.e., excluding the true model) of the Leave-One-Out model-as-truth experiments. Left panel: the model index used as pseudo-observations is on x-axis, whereas the overlap scores (in %) between the reconstructed PDF and the true PDF read in y-axis. Each colored marker corresponds to a score strategy as denoted in the legend. Right panel: boxplots spanning for each score (keeping the same colour code) the 25 th , median, and 75 th percentiles of the distributions of overlap scores for all the model-as-truth experiments combined. The whiskers show the minimum and maximum values excluding the outliers denoted with the black diamonds. Small yellow triangle indicate the means.

  15 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

  16 https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

  https://doi.org/10.5194/egusphere-2023-2649 Preprint. Discussion started: 1 December 2023 c Author(s) 2023. CC BY 4.0 License.

Table 1 .

 1 Perfect model results for model 8 as the pseudo-observation. For each score (in row): weights (in %) associated with the eleven models. The sum of the weights per row is 100. Values in bold highlight the model with the highest weight for each score.

		0	1	2	3	4	5	6	7	8	9	10
	Model democracy	9.1 9.1 9.1	9.1	9.1 9.1 9.1	9.1	9.1	9.1	9.1
	Climatologies	6.2	9	9.1	9.8	9.4 9.5 8.7	9.3	9.7	10	9.2
	Best single model	0	0	0	0	0	0	0	0	0	100	0
	CME ClimWIP	0.1	0	0	16.2 0.3 0.5 0.1 30.4 31.2 21.2	0
	CME best punctual	4	2.3	2	10	8.3 6.5 4.8 12.3 16.8	24	9.3
	CME best persistent 3.3 1.9 1.7	8.5	7.3	6	5	12.5 17.5 28.1 8.1

Acknowledgements. The authors thank the financial support for the AMIGOS project from the Région Bretagne (ARED fellowship for PLB) and from the ISblue project (Interdisciplinary graduate school for the blue planet, ANR-17-EURE-0015) co-funded by a grant from the French government under the program Investissements d'Avenir. The authors also want to thank Maxime Beauchamp and Noémie Le Carrer for their helpful comments on the study.

Code and data availability. Python codes are available in the following GitHub open repository: https://github.com/pilebras/AnDA_weight_ idealAMOC/tree/main under the GNU license. They include the code for generating the data (using an idealized AMOC model developped in [START_REF] Sévellec | Millennial variability in an idealized ocean model: Predicting the AMOC regime shifts[END_REF], performing the experiments using AnDA (adapted from the code available at https://github.com/ptandeo/AnDA) and obtaining the figures of the paper.

Author contributions. All authors contributed equally to the experimental design. PLB wrote the article and performed the numerical experiments and the analyses. FS, PT, JR and PA helped on the redaction of the paper.

Competing interests. At least one of the (co-)authors is a member of the editorial board of Nonlinear Processes in Geophysics. 7d), represents the approximate marginal likelihood of the observations. This is estimated within the AnDA process using successive innovation likelihoods (i.e., the CME values, e.g., [START_REF] Carrassi | Estimating model evidence using data assimilation[END_REF][START_REF] Pulido | Stochastic parameterization identification using ensemble Kalman filtering combined with maximum likelihood methods[END_REF][START_REF] Metref | Estimating model evidence using ensemble-based data assimilation with localization -The model selection problem[END_REF][START_REF] Tandeo | A review of innovation-based methods to jointly estimate model and observation error covariance matrices in ensemble data assimilation[END_REF]). This 475 approach offers two key advantages. Firstly, it avoids the need for calculating a climatological distribution, which can be challenging to estimate accurately and computationally expensive. Secondly, it provides more informative insights into current conditions, including forecast states with their uncertainties. This is especially useful for dealing with nonlinearities in dynamic systems, instead of using globally defined climatological distributions [START_REF] Metref | Estimating model evidence using ensemble-based data assimilation with localization -The model selection problem[END_REF]. Furthermore, the denominator in Eq. (7d) serves as a prior that contains information about inter-model dependency. This provides more valuable insights compared to a uniform prior commonly used in standard BMA approaches (e.g., George 2010; Garthwaite and Mubwandarikwa 2010).