Eric Adjakossa
email: ericadjakossah@gmail.com

Yannig Goude
email: yannig.goude@edf.fr

Olivier Wintenberger
email: olivier.wintenberger@upmc.fr

Kalman Recursions Aggregated Online

Keywords: online aggregation, Kalman filter, experts ensemble

In this article, we aim at improving the prediction of expert aggregation by using the underlying properties of the models that provide expert predictions. We restrict ourselves to the case where expert predictions come from Kalman recursions, fitting state-space models. By using exponential weights, we construct different algorithms of Kalman recursions Aggregated Online (KAO) that compete with the best expert or the best convex combination of experts in a more or less adaptive way. We improve the existing results on expert aggregation literature when the experts are Kalman recursions by taking advantage of the second-order properties of the Kalman recursions. We apply our approach to Kalman recursions and extend it to the general adversarial expert setting by state-space modeling the errors of the experts. We apply these new algorithms to a real dataset of electricity consumption and show how it can improve forecast performances comparing to other exponentially weighted average procedures.

Introduction

The aim of this paper is to aggregate Kalman recursions in an online setting in order to increase the accuracy of the prediction. We observe (y t ) sequentially through time t ≥ 1 and predictions ŷ(m) t , 1 ≤ m ≤ M , issued from Kalman recursions at time t ≥ 1. Here M ≥ 0 denotes the number of different Kalman recursions used as experts. The Kalman recursions are imbedded into a statespace model (see Section 2.1 for a formal definition). We introduce different Kalman recursions Aggregated Online (KAO) procedures that compute recursively weights ρ We obtain bounds on the regret of KAO algorithms that are similar to the ones encountered in the literature. The book of reference on aggregation is undoubtedly the book from [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF] and we refer to it for classical regret aggregation bounds. The novelty of our approach is to derive regret bounds directly on the cumulative quadratic predictive risk as defined by [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF]. The predictive risk or risk of prediction of a predictor ŷ ∈ F t-1 is defined as

L t (ŷ) = E (ŷ -y t ) 2 | F t-1 , a.s., t ≥ 1 , (1)
where (F t ) is the natural filtration of the past response variables σ(y s ; 0 ≤ s ≤ t) = F t , t ≥ 0.

The risk of prediction arises naturally when dealing with Kalman recursions. Indeed, Kalman recursions are online algorithms that provide the best linear predictions in gaussian state-space models. We refer to the classical monograph from [START_REF] Durbin | Time series analysis by state space methods[END_REF] for details. The cumulative predictive risk of a recursive algorithm predicting ŷt at each time t ≥ 1 is the sum

T t=1
L t (ŷ t ) up to the horizon T ≥ 1. Our regret bounds on KAO algorithm predictions (ŷ t ) are a.s.

deterministic bound called respectively model selection regret or aggregation regret and defined as

R S T (m) ≥ T t=1 L t (ŷ t ) -L t (ŷ (m) t ) 1 ≤ m ≤ M , (2)
R A T (π) ≥ T t=1 L t (ŷ t ) -L t M m=1 π (m) ŷ(m) s , (3)
for any vector of weights π := (π (m) ) 1≤m≤M in the simplex. We suppress the dependence in m and π in R S t and R A t when the regret bounds are uniform in m and π, respectively.

Regret bounds on the cumulative predictive risk have attracted some interest since [START_REF] Audibert | Regret bounds and minimax policies under partial monitoring[END_REF] showed that the classical EWA algorithm from [START_REF] Vovk | Aggregating strategies[END_REF] does not achieve a fast rate model selection regret even in the most favorable iid case. The fast rate model selection regret was first proved by the BOA algorithm in the iid setting with strongly convex loss in [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF] and then extended to any stochastic setting and exp-concave risk in [START_REF] Gaillard | Sparse accelerated exponential weights[END_REF]. For adaptative procedures, we improve their optimal regret R S T = O(log M + log log T + x) , T ≥ 1 with probability 1 -e -x , x > 0, to the a.s. bound

R S T = O(log M + log log T ) T ≥ 1 .
This optimal regret bound holds when the observations satisfy some unbounded state-space model defined in the next section 2.1. That the observations satisfy such a model is very unlikely in prac-tice; It is the price to pay to guarantee a.s. optimal regret bounds on the cumulative predictive risk (1) for unbounded responses. Existing regret bounds such as the one of [START_REF] Gaillard | Sparse accelerated exponential weights[END_REF] requires the boundedness of the response.

We present simulations and applications in cases where our assumptions are certainly not satisfied. Our new aggregation procedure improves the state of the art methods in aggregation such as MLPoly of Gaillard et al. (2014).

Preliminaries

State-space models

Assume that we observe (y t , X t ) with y t ∈ R the variable of interests and X t ∈ R d is the predictable design, i.e. X t ∈ F t-1 , t ≥ 1. Notice that the design X t can be either deterministic or random. We consider a collection of M ≥ 1 experts ŷ(m) ∼ N (0, Q (m) ) constitute an iid sequence and θ

(m) 0 ∈ R d is deterministic. The sequence (θ (m) t
) is a Gaussian Markov chain and admits the representation

θ (m) t = K (m) t θ m 0 + t-1 k=0 K (m) k z (m) t-k , t ≥ 1 , 1 ≤ m ≤ M . ( 5
)
It converges weakly to a stationary solution if and only if ρ(K (m) ) < 1, the spectral radius of m) is smaller than one. This assumption is not required in this work. Actually one of the most popular state models is the dynamic setting where one considers random walk under K (m) = I d , the identity matrix of R d . We notice that as the state model is hidden (latent, not observed),

K (
any assumption on the state recursion, such as the Gaussian assumption, is not restrictive for the observations (y t , X t ).

Our main assumption is the following one.

(H) The vectors (y t , θ (m) t

) constitute a Gaussian sequence,

E[y t | θ (m) t ] = X t θ (m) t , t ≥ 1 , 1 ≤ m ≤ M ,
Algorithm 1 Kalman recursion in the m-th state-space model

Parameters:

The matrices Q (m) and K (m) . Initialization: The matrix P (m) 0

and the vector θ(m) 0 . Recursion: For each iteration t = 1, . . . , T do:

θ (m) t+1 = K (m) θ (m) t + 1 X t P (m) t X t + 1 P (m) t X t (y t - ŷ(m) t ) , P (m) t+1 = K (m) P (m) t - 1
X t P (m) t X t + 1 P (m) t X t X t P (m) t K (m) + Q (m) ŷ(m) t+1 = X t+1 θ (m) t+1 .
and the conditional variance

σ 2(m) := Var(y t | θ (m) t
) > 0 is constant through time and known.

Condition (H) have different consequences upon the observations (y t , X t ). The first obvious one is that (y t ) constitutes a Gaussian sequence. The second one is that y t satisfies the linear model

y t = X t θ (m) t + ε (m) t , t ≥ 1 , 1 ≤ m ≤ M .
The gaussian property of the couple (y t , θ

(m) t ) ensures that ε (m) t
is a gaussian random variable with mean zero and variance σ 2(m) = Var(y t | θ (m) t

) independent of t ≥ 1. A direct implication from the expression (5) is the following mean-variance identity.

Proposition 2.1. Under Condition (H) the following mean-variance identity holds for all 1 ≤ m ≤ M and t ≥ 1:

E[y t ] = E[X t θ (m) t ] = E[X t ] K (m) t θ m 0 , Var(y t ) = Var[X t θ (m) t ] + σ 2(m) = E X t t-1 k=0 K (m) k Q (m) K (m) k X t + σ 2(m) .
The static state-space model setting corresponds to the case where

Q (m) = 0 so that Var(y t ) = σ 2(m) = σ 2 , 1 ≤ m ≤ M , t ≥ 1.

The Kalman recursion

For the sake of completeness, we recall the Kalman recursion associated with the m-th statespace model in Algorithm 1. For details on the Kalman recursion we refer to the monograph from [START_REF] Durbin | Time series analysis by state space methods[END_REF]. We notice that the Kalman recursion does not require any inversion of matrices. Each iteration has thus a O(d 2 ) computational cost. Moreover it does not require the knowledge of the parameters σ 2(m) > 0. In addition, in many cases X t is in fact a vector of

size d = M m=1 d m that stacks M vectors X (m) t ∈ R dm . In this case one considers d m sparse vectors θ (m) t
and identifies them with their non-null components θ

(m) t ∈ R dm .
Then the space equation is written as In the static case when K (m) = I and Q (m) is the null matrix then, using the Shermann-

y t = X (m) t θ (m) t + ε (m) t , t ≥ 1 , 1 ≤ m ≤ M ,
Morrisson formula, we have the alternative recursion for R (m) t

, the inverse of

P (m) t , R (m) t+1 = R (m) t + X t X t .
When P (m) 0

is taken equals to 1/λ (m) I d , for some λ (m) > 0, then the estimator computed recursively using the Kalman recursion coincides with the Ridge estimator θ(m) t = arg min

θ∈R d t s=1 y s -X s θ 2 + λ (m) 2 θ - θ(m) 0 2 2 .
This equivalence has been first established by [START_REF] Diderrich | The kalman filter from the perspective of goldberger?theil estimators[END_REF].

Notice that there is no assumption on the dependence among the Kalman recursions. Otherwise, it was possible to consider a Kalman recursion over the stack of the models in a dM dimensional state-space model. However, this approach is not practical when the computational cost O((dM ) 2 ) of the complete Kalman recursion is prohibitive because M is too large. The aim of this work is to show that this ideal procedure, that is uncertain in practice when the dependence among the recursions has to be estimated, can be easily overcome by a simple aggregation procedure over M Kalman recursions.

Examples

We provide some classical examples of state-space models satisfying (H):

1. The static iid setting: this degenerate case coincides with the usual gaussian linear model for fixed or random (iid) design (X t ). We assume the relation = 0 for all 1 ≤ m ≤ M . Then the mean identity of Proposition 2.1 is automatically satisfied as E[y t ] = 0 for all t ≥ 1. The variance identity requires that m ) for any 1 ≤ m, m ≤ M and t ≥ 1. The Kalman recursion can be used for tracking the signal (y t ) on different explanatory variables

y t = X t θ (m) t + ε (m) t , t ≥ 1, associated with state equations θ (m) t = θ (m) t-1 = • • • = θ (m) 0 (Q (m) = 0 for 1 ≤ m ≤ M ).
E[X t Q (m) X t ] = E[X t Q (m ) X t ] and σ 2(m) = σ 2(
X (m) t stacked in X t .
3. The expert setting: We consider the case where we have M deterministic experts without any information about their generation process. This situation is very common in real-life applications as the forecast can come from different sources (physical models, different data sources, different machine learning models).

For each 1 ≤ m ≤ M expert we stack its prediction

f m,t ∈ R in X (m) t
together with the intercept and the past error e m,t-1 = (y t-1 -f m,t-1 ), i.e.,

X (m) t = (1, f m,t , e m,t-1 ), t ≥ 1 .
and each state-space model is defined by the state equation:

θ (m) t = K (m) θ (m) t-1 + z (m) t , t ≥ 1.

Kalman recursions Aggregated Online (KAO) algorithm

Consider the state-space models (coinciding with Equation (4) for the mth state equation,

1 ≤ 1 ≤ M )    y t = X t θ (m) t + ε (m) t θ (m) t = K (m) θ (m) t-1 + z (m) t , t ≥ 1, 1 ≤ m ≤ M.
Recall that under (H) we have

E[y t | m] := E[y t | z (m) t , . . . , z (m) 1 , F t-1 ] = X t θ (m) t , 1 ≤ m ≤ M.
We aggregate Kalman recursions using a version of the exponentially weighted average forecaster defined as

ŷt = M m=1 ρ (m) t ŷ(m) t (7) with ρ (m) t ≥ 0 for 1 ≤ m ≤ M , M m=1 ρ (m) t = 1 and ŷ(m) t = X t θ (m) t
is the mth Kalman forecaster of y t .

Convex properties

The ability to find rapidly the solution of an optimization problem depends heavily on the convex properties of the objective function. In our case, the objective function is the conditional risk defined in (1) as Thus one cannot expect that our procedure converges in general and we rather study its regrets R S t and R A t defined in Equations ( 2) and (3) as the model selection and aggregation regret, respectively. The objective function is

L t (ŷ) = E t-1 [(ŷ -y t ) 2 ].
t s=1 L s M m=1 π (m) ŷ(m) s for any (π (m) ) 1≤m≤M in the canonical basis or in the simplex, i.e. π (m) ≥ 0 such that M m=1 π (m) = 1.
The optimal rates of convergence in the model selection and the aggregation problems depend on the convex properties of the objective function and the observation of an approximation of the gradients. As the objective functions are convex, we will extensively use the gradient trick which consists to bounding the regret with the linearized risks L (m)

s = L s (ŷ s )(ŷ (m) s -ŷs ) as L t (ŷ t ) -L t M m=1 π (m) ŷ(m) t ≤ - M m=1 π (m) L (m) t . ( 8
)
Fast rates of convergence could be obtained easily if the objective function was strongly convex.

Despite we use the square loss, it is not the case since the Hessian matrix

2 t s=1 (ŷ (m) s ) 1≤m≤M (ŷ (m) s ) 1≤m≤M
of the objective function is a sum of rank-one matrices which are very unlikely to converge in any non-stationary settings. This issue is bypassed in online convex optimization thanks to the notion of exp-concavity extensively studied by [START_REF] Hazan | Introduction to online convex optimization[END_REF].

Definition 3.1. A loss function is η-exp-concave (with η > 0) on some convex set Y if the function F (y) = exp(-η (y)) is concave for all y ∈ Y.
We need to find out for which values of η the conditional risks (1) are exp-concave. Moreover, we can use the exp-concave property of the risk L t to refine the gradient trick.

Theorem 3.1. Assume that it exists D > 0 such that |ŷ

(m) t -µ t | ≤ D a.s. 1 ≤ m ≤ M , t ≥ 1 with µ t = E t-1 [y t ]. Then the conditional risk L t is a.s. (2D 2 ) -1 -exp-concave for any y = M m=1 π (m) ŷ(m) t , t ≥ 1, with (π (m) ) 1≤m≤M in the simplex. Moreover if the linearized risk satisfies |L (m) t | ≤ G (m)
for any t ≥ 1, 1 ≤ m ≤ M , then we have

L t (ŷ t ) -L t (ŷ (m) t ) ≤ -L (m) t -η (m) L (m) t 2 , with η (m) = 1 8(2G (m) ∨ D 2 ) , 1 ≤ m ≤ M . The refined linearized risk L (m) t +η (m) L (m) t 2 is called the surrogate risk. It is itself exp-concave
due to the quadratic term whereas the linearized risk cannot be exp-concave.

Proof. Let 0 < η ≤ 1 (2D 2 ) . Consider the function ϕ η (y) = e -ηLt(y) for y = M m=1 π (m) ŷ(m) t and
(π (m) ) 1≤m≤M in the simplex. The function ϕ is a.s. at least twice differentiable and we have

ϕ η (y) = -2ηϕ η (y) 1 -2η(y -µ t ) 2 .
using the derivation under the integral sign. For θ ∈ Θ, we get the concavity since ϕ η (y) ≤ 0 as

2η(y -µ t ) 2 ≤ 2η( M m=1 π (m) (ŷ (m) t -µ t )) 2 ≤ 2ηD 2 = 1
and the first assertion follows.

We proceed as in the proof of Lemma 4.2 of [START_REF] Hazan | Introduction to online convex optimization[END_REF] considering

γ (m) = 1 2(2G (m) ∨ D 2 ) ≤ η. One deduces from the concavity property of ϕ γ (m) that ϕ γ (m) (y) -ϕ γ (m) (z) ≤ ϕ γ (m) (z)(y -z) which, taking y = ŷ(m) t and z = ŷt , provides exp(-γ (m) L t (ŷ (m) t )) -exp(-γ (m) L t (ŷ t )) ≤ -γ (m) L t (ŷ t ) exp(-γ (m) L t (ŷ t ))(ŷ (m) t -ŷt ).
One deduces that

γ (m) (L t (ŷ t ) -L t (ŷ (m) t )) ≤ log(1 -γ (m) L t (ŷ t )(ŷ (m) t -ŷt )). Using the relation log(1-z) ≤ -z - 1 4 z 2 that holds for any |z| ≤ 1/4 applied on |γ (m) L t (ŷ t )(ŷ (
(m) t , 1 ≤ m ≤ M . Recursion: Do: ρ (m) t+1 = exp -η X t P (m) t X t + σ 2(m) ρ (m) t M m=1 exp -η X t P (m) t X t + σ 2(m) ρ (m) t ŷt+1 = M m=1 ρ (m) t+1 ŷ(m) t+1 . ŷt )| ≤ 1/4 one obtains γ (m) (L t (ŷ t ) -L t (ŷ (m) t )) ≤ γ (m) L t (ŷ t )(ŷ t - ŷ(m) t ) - 1 4 (γ (m) L t (ŷ t )(ŷ t - ŷ(m) t )) 2
and the second assertion follows.

KAO for model selection

In this Section we assume the exp-concavity of the conditional risks and we adapt the classical analysis of the Exponentially Weighted Average (EWA) algorithm of [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF] to our setting. The aggregation procedure, called KAO, is described in Algorithm 2.

KAO achieves the optimal rate for model selection.

Theorem 3.2. Under assumption (H) and if it exists D > 0 such that |ŷ

(m) t -µ t | ≤ D a.s. 1 ≤ m ≤ M , t ≥ 1 then KAO for model selection with η = 1 (2D 2 )
achieves the regret bound

R S t (m) ≤ -2D 2 log(ρ (m) 0 ) 1 ≤ t ≤ T , 1 ≤ m ≤ M .
We note that the classical EWA algorithm satisfies a similar regret bound under the stronger assumption

|ŷ (m) t -y t | ≤ D, 1 ≤ m ≤ M, t ≥ 1, a.s.
which never holds in our Gaussian setting. One usual way to bypass this well-known restriction of EWA is to use a doubling trick which deteriorates the regret bound, see [START_REF] Cesa-Bianchi | Prediction, learning, and games[END_REF] for more details.

Proof. The proof is standard and follows the line of the proof of the EWA regret in Cesa-Bianchi and Lugosi ( 2006). The crucial step consists in identifying the conditional risk of any Kalman prediction ŷ(m) t under (H). We have the following Lemma Lemma 3.1. Under (H) we have the identity ŷ(m)

t = E t-1 [E[y t | m]].
Proof. The Kalman recursion produces the best linear prediction which is equal to the conditional expectation in the gaussian case. Then we have ŷ(m)

t = E t-1 [X t θ (m) t ] = E t-1 [E[y t |m]] by definition.
We have explicitely

L s (ŷ (m) s ) = E s-1 [(y s -E s-1 [E[y s | m]]) 2 ] = E s-1 [(y s -E[y s | m]) 2 ] + E s-1 [(E[y s | m] -E s-1 [E[y s | m]]) 2 ] +2 E s-1 [(y s -E[y s | m])(E[y s | m] -E s-1 [E[y s | m]])] = E s-1 [(y s -X t θ (m) t ) 2 ] + E s-1 [(X t (θ (m) t - θ(m) t ) 2 ] = σ 2(m) + X s P (m) s X s ,
since the third term of the sum is zero and since

E s-1 [(X t (θ (m) t - θ(m) t ) 2 ] = X s P (m) s X s thanks to
the Kalman recursion properties in the gaussian case. We also have, using the exp-concavity of L t and Jensen inequality,

e -ηLt(ŷt) = e -ηLt( M m=1 ρ (m) t ŷ(m) t ) ≥ M m=1 ρ (m) t e -ηLt(ŷ (m) t ) , ≥ M m=1 ρ (m) 0 e -η t-1 s=1 Ls(ŷs) e -ηLt(ŷ (m) t ) M m=1 ρ (m) 0 e -η t-1 s=1 Ls(ŷs) ≥ M m=1 ρ (m) 0 e -ηR S t-1 (m) e -ηLt(ŷ (m) t ) M m=1 ρ (m) 0 e -ηR S t-1 (m) ,
multiplying by e η t-1 s=1

Ls(ŷ (m) s

) above and below the fraction. We get the recursive relation

1 = M m=1 ρ (m) 0 ≥ M m=1 ρ (m) 0 e -ηR S t-1 (m) ≥ M m=1 ρ (m) 0 e -ηR S t (m)
and the desired result follows.

KAO for aggregation

In the case where the best expert is not worthy of confidence, it is generally much more interesting to compete with the best convex combination of the experts at hand. In this context, the aim is to provide a bound on the regret for aggregation R A t (π) where π := (π (m) ) 1≤m≤M

L (m) t = (9) ρ (m) t+1 = exp -ηL (m) t ρ (m) t M m =1 exp -ηL (m) t ρ (m ) t ŷt+1 = M m=1 ρ (m) t+1 ŷ(m) t+1 .
belongs to the simplex. As the conditional risk L s is a convex function that is differentiable, one applies the gradient trick and we consider an explicit biased version of the linearized risk

L (m) t = X t P (m) t X t + σ 2(m) -(ŷ t - ŷ(m) t ) 2 - M m =1 ρ (m ) t X t P (m ) t X t + σ 2(m ) -(ŷ t - ŷ(m ) t ) 2 . (9)
By convention L

(m) t | ≤ G a.s. for 1 ≤ t ≤ T , 1 ≤ m ≤ M .
Then KAO for aggregation starting with ρ

(m) 0 = 1/M and η = 1 G 2 log M t satisfies the regret bound R A t ≤ G 2t log M , 1 ≤ t ≤ T. ( 10
)
The regret bound matches the optimal bound for M ≥ √ t. Note that the boundedness assumption on L (m) t involves only the predictions and does not require to bound (y t ).

Proof. Since L s is convex and differentiable, we apply the gradient trick

R A t (π) ≤ - t s=1 M m=1 π (m) L (m) s .
Moreover, the expression of L s (ŷ s )ŷ (m) s can be developed as

L s (ŷ s )ŷ (m) s = 2 E s-1 [(ŷ s -y s )ŷ (m) s ] = 2ŷ s ŷ(m) s + E s-1 [(y s -ŷ(m) s ) 2 ] -ŷ(m) 2 s -E s-1 [y 2 s ] = E s-1 [(y s -ŷ(m) s ) 2 ] -(ŷ s -ŷ(m) s ) 2 + ŷ2 s -E s-1 [y 2 s ] = X s P (m) s X s + σ 2(m) -(ŷ s -ŷ(m) s ) 2 + ŷ2 s -E s-1 [y 2 s ].
Since the two last summands of L s (ŷ s )ŷ (m) s do not depend on m we obtain the identity

L (m) t = L s (ŷ s )(ŷ (m) s -ŷs ) = (9). As it exists G > 0 satisfying |L (m) t | ≤ G,
ρ (m) s = exp(-η s-1 r=1 L (m) t )ρ (m) 0 M m =1 exp(-η s-1 r=1 L (m ) t )ρ (m ) 0 we get log   M m=1 exp -η s r=1 L (m) r ρ (m) 0 M m =1 exp -η s-1 r=1 L (m ) r ρ (m ) 0   ≤ η 2 2 G 2 .
Then, by summing over s, a telescoping sum appears and leads to

1 η log M m=1 exp -η t s=1 L (m) s ρ (m) 0 ≤ ηt G 2 2 . Moreover, exp -η t s=1 L (m) s ρ (m) 0 ≤ M m=1 exp -η t s=1 L (m) s ρ (m) 0 , which leads to M m=1 π (m) t s=1 -L (m) s + log ρ (m) 0 η ≤ 1 η log M m=1 exp -η t s=1 L (m) s ρ (m) 0 .
Combining those bounds, we obtain

R A t (π) ≤ t s=1 M m=1 -π (m) L (m) s ≤ M m=1 π (m) -log ρ (m) 0 η + ηt G 2 2 .
We get the desired result noticing that ρ (m) 0

= 1/M and the optimal choice of η = 1 G 2 log M t .

We notice that a unique learning rate yields a uniform regret bound, independent of π. We also notice that we can use the gradient trick despite we only observed a biased version of the Algorithm 4 KAO with multiple learning rates

Parameters:

The variances σ 2(m) , the weights ρ(m) 0 and the learning rates η (m) , 1 ≤ m ≤ M .

Initialization:

The initial weights ρ

(m) 0 = η (m) ρ(m) 0 /( M m =1 η (m ) ρ(m ) 0 ), 1 ≤ m ≤ M .
For each iteration t = 1, . . . , T : Inputs: The Kalman predictions ŷ(m) t+1 and the matrices P

(m) t , 1 ≤ m ≤ M . Recursion: Do: L (m) t = (9) ρ (m) t+1 = exp -η (m) L (m) t (1 + η (m) L (m) t ) ρ (m) t M m =1 exp -η (m) L (m ) t (1 + η (m ) L (m ) t ) ρ (m ) t ŷt+1 = M m=1 ρ (m) t+1 ŷ(m) t+1 .
linearized risk thanks to the exponential form of the weights that are not sensitive to the bias.

Online tuning of the learning rates

The theoretical guarantees on the regret for the model selection and the aggregation problems do not hold for the same algorithm. The gradient trick is a crucial step in the proof of the regret for the aggregation problem. However, the fast rate for the model selection does not hold for the gradient-based EWA since the linearized risk cannot be exp-concave. In order to bypass this issue, we adapt the approach of [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF] to our setting. The first step is to use a surrogate loss of the form L

(m) t 1 + ηL (m) t
where the quadratic part yields exp-concavity. The second step is to use a multiple learning rates version of KAO as described in Algorithm 4 in Section 4.1 where we show that multiple learning rates can be easily tuned online.

Multiple learning rates for KAO

In the context of expert aggregation, it is well known that using multiple learning rates help to increase the prediction accuracy, see Gaillard et al. (2014) and [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF]. Here we aim to provide a multiple learning rates version for KAO in a similar way than the multiple learning rate version of the BOA procedure (see [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF]). The following theorem provides regret bounds both for model selection and aggregation on the same algorithm.

η (m) = 1 G   - log ρ(m) 0 t ∧ 1 2   is bounded as R A t (π) ≤ 2G M m=1 π (m) -log ρ(m) 0 t -log ρ(m) 0 , 1 ≤ t ≤ T . ( 11
)
If moreover there exists D > 0 such that |ŷ

(m) t -µ t | ≤ D a.s. 1 ≤ m ≤ M , t ≥ 1 then the model selection regret of KAO with multiple learning rates η (m) = 1 8(2G ∨ D 2 ) for all 1 ≤ m ≤ M is bounded as R S t (m) ≤ -8(2G ∨ D 2 ) log ρ(m) . ( 12
)
Proof. We start by applying the gradient trick as in (8) inferring that

R A t (π) ≤ - t s=1 M m=1 π (m) L (m) s . Moreover, as x -x 2 is 1-exp-concave for x > 1/2, Jensen's inequality implies that E exp X -X 2 ≤ exp E[X] -E[X] 2 = 1 (13)
for any centered random variable X such that X ≥ -1/2 a.s. We notice that η

(m) = 1 G - log ρm 0 t ∧ 1 2
satisfies the relation

η (m) L (m) t ≤ 1 2 . Denoting ρ(m) t = exp - t-1 s=1 η (m) L (m) s 1 + η (m) L (m) s ρ(m) 0 M m =1 exp - t-1 s=1 η (m ) L (m ) s 1 + η (m ) L (m ) s ρ(m ) 0 , we have the identity ρ(m) t = ρ (m) t η (m) × 1 M m =1 ρ (m ) t /η (m )
, which leads to

M m=1 ρ(m) t η (m) L (m) t = 0.
Thus the random variable

η (m) L (m) t 1≤m≤M
is therefore centered for the distribution ρ(m) t 1≤m≤M

, and we have

M m=1 ρ(m) t exp -η (m) L (m) t 1 + η (m) L (m) t ≤ 1, ( 14
) since -η (m) L (m) t
> -1/2 a.s. for all 1 ≤ m ≤ M and 1 ≤ t ≤ T , using (13). By putting the expression of ρ(m) t into Equation ( 14), we have

M m=1 ρ(m) 0 exp - t s=1 η (m) L (m) s 1 + η (m) L (m) s ≤ M m=1 ρ(m) 0 exp - t-1 s=1 η (m) L (m) s 1 + η (m) L (m) s , which implies for 1 ≤ t ≤ T that M m=1 ρ(m) 0 exp - t s=1 η (m) L (m) s 1 + η (m) L (m) s ≤ 1,
since by convention that L (ŷ

(m) 0 ) = 0. Thus, for 1 ≤ m ≤ M , - t s=1 η (m) L (m) s 1 + η (m) L (m) s ≤ -log ρ(m) 0 , (15)
by applying the logarithm function using the previous inequality. We have

- M m=1 π(m) t s=1 η (m) L (m) s ≤ M m=1 π(m) η (m) 2 t s=1 L (m) s 2 -log ρ(m) 0 for π(m) = π (m) /η (m) M m =1 π (m ) /η (m ) , 1 ≤ m ≤ M .
Multiplying with

M m =1 π (m ) /η (m ) we obtain - M m=1 π (m) t s=1 L (m) s ≤ M m=1 π (m) η (m) t s=1 L (m) s 2 - log ρ(m) 0 η (m) ≤ M m=1 π (m) η (m) G 2 t - log ρ(m) 0 η (m)
and the desired result on R A t (π) follows from the specific choice of η (m) . The regret bound on R S t (m) follows by an application of Theorem 3.1 on the last bound specified for π in the canonical basis

t s=1 -L (m) s -η (m) L (m) s 2 ≤ π (m) - log ρ(m) 0 η (m) for all 1 ≤ m ≤ M as η (m) = 1 8(2G ∨ D 2 )
.

Algorithm 5 KAO with adaptive multiple learning rates

Parameters:

The variances σ 2(m) , 1 ≤ m ≤ M .

Initialization: Any initial weights ρ(m)

0 > 0 such that M m=1 ρ(m) 0 = 1, 1 ≤ m ≤ M .
For each iteration t = 1, . . . , T : Inputs: The Kalman predictions ŷ(m) t+1 and the matrices P

(m) t , 1 ≤ m ≤ M . Recursion: Do: L (m) t = (9) η (m) t = -log ρ(m) 0 1 + t s=1 L (m) s 2 ρ (m) t+1 = η (m) t exp -η (m) t t s=1 L (m) s 1 + η (m) s-1 L (m) s ρ(m) 0 M m =1 η (m ) t exp -η (m ) t t s=1 L (m ) s 1 + η (m ) s-1 L (m ) s ρ(m ) 0 ŷt+1 = M m=1 ρ (m) t+1 ŷ(m) t+1 .
However the regrets bound do not apply on KAO with the same learning rates. The slow rate aggregation regret bound holds for a O(1/ √ t) learning rate whereas the fast rate model selection regret bound holds for a constant learning rate.

Adaptive multiple learning rates

Multiple learning rates are easily adaptable as in Algorithm 5. Moreover, a single algorithm with unique adaptive learning rates achieves optimal regret bounds for both model selection and aggregation problems as for the BOA algorithm developed by [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF] and refined by [START_REF] Gaillard | Efficient online algorithms for fast-rate regret bounds under sparsity[END_REF].

Theorem 4.2. Under assumption (H) suppose there exist G (m) > 0 and D > 0 such that

|L (m) t | ≤ G (m) and |ŷ (m) t -µ t | ≤ D a.s. for 1 ≤ t ≤ T , 1 ≤ m ≤ M .
Then the regret of KAO with adaptive multiple learning rates such that η

(m) t-1 L (m) t < 1/2 for any 1 ≤ t ≤ T and 1 ≤ m ≤ M is bounded as R A t (π) ≤ M m=1 π (m) G (m) (3 + G (m) ) √ t + 1 -log ρ(m) 0 + r (m) t , R S t (m) ≤ 8(2G (m) ∨ D 2 )(3 + G (m) ) -log ρ(m) 0 + r (m) t 2 + -log ρ(m) 0 + r (m) t .
where r

(m) t = log log e 1/4 + G (m) √ t + 1 -log ρ(m) 0 .
Remark 4.1. The leading constant is proportional to G (m) 2 . It is not optimal and can be reduced to G (m) by refining the adaptive learning rates as in [START_REF] Cesa-Bianchi | Improved second-order bounds for prediction with expert advice[END_REF].

Proof. By adapting the inequality ( 14) as (η

(m) t-1 L (m) t ) 1≤m≤M is centered for ρ(m) t 1≤m≤M
, where ρ(m)

t = exp -η (m) t-1 t-1 s=1 L (m) s 1 + η (m) s-1 L (m) s ρ(m) 0 M m =1 exp -η (m ) t-1 t-1 s=1 L (m ) s 1 + η (m ) s-1 L (m ) s ρ(m ) 0 , for any t ≥ 2 we have M m=1 ρ(m) 0 exp -η (m) t-1 t s=1 L (m) s 1 + η (m) s-1 L (m) s ≤ M m=1 ρ(m) 0 exp -η (m) t-1 t-1 s=1 L (m) s 1 + η (m) s-1 L (m) s . (16) Since x ≤ x α + α -1 (α -1) for x ≥ 0 and α ≥ 1, by setting α = η (m) t-2 η (m) t-1 and x = exp -η (m) t-1 t-1 s=1 L (m) s 1 + η (m) s-1 L (m) s ,
we have for any t ≥ 2 the relation exp -η

(m) t-1 t-1 s=1 L (m) s 1 + η (m) s-1 L (m) s ≤ exp -η (m) t-2 t-1 s=1 L (m) s 1 + η (m) s-1 L (m) s + η (m) t-2 -η (m) t-1 η (m) t-2 , which leads to M m=1 ρ(m) 0 exp -η (m) t-1 t s=1 L (m) s 1 + η (m) s-1 L (m) s ≤ M m=1 ρ(m) 0 exp -η (m) t-2 t-1 s=1 L (m) s 1 + η (m) s-1 L (m) s + M m=1 ρ(m) 0 η (m) t-2 -η (m) t-1 η (m) t-2 . ( 17
)
Using a recursion argument on t ≥ 2 on Equation ( 17) yields

M m=1 ρ(m) 0 exp -η (m) t-1 t s=1 L (m) s 1 + η (m) s-1 L (m) s ≤ M m=1 ρ(m) 0 exp -η (m) 0 L (m) 1 1 + η (m) 0 L (m) 1 + t-1 s=1 M m=1 ρ(m) 0 η (m) s-1 -η (m) s η (m) s-1 . (18)
Moreover, we have

t-1 s=1 η (m) s-1 -η (m) s η (m) s-1 ≤ t-1 s=1 η (m) s-1 η (m) s dx x ≤ η (m) 0 η (m) t-1 dx x ≤ log η (m) 0 η (m) t-1
the estimate of the ratio

η (m) 0 η (m) t-1 = 1 + t-1 s=1 L (m) s 2 ≤ G (m) √ t + 1, for G (m) ≥ 1. and -η (m) 0 L (m) 1 1 + η (m) 0 L (m) 1 ≤ 1/4 . Equation (18) implies -η (m) t-1 t s=1 L (m) s 1 + η (m) s-1 L (m) s ≤ -log ρ(m) 0 + r (m) t ,
and we obtain similarly than above that for any π we have

- M m=1 π (m) t s=1 L (m) s ≤ M m=1 π (m) t s=1 η (m) s-1 L (m) s 2 + -log ρ(m) 0 + r (m) t η (m) t-1
.

In order to bound the second order term

t s=1 η (m) s-1 L (m) s 2 we denote V t = 1 + t s=1 L (m) s 2 so that η (m) s-1 L (m) s 2 = -log ρ(m) 0 V s -V s-1 V s-1 = -log ρ(m) 0 √ V s + √ V s-1 V s-1 ( √ V s - √ V s-1 ) = -log ρ(m) 0 V s /V s-1 + 1 ( √ V s - √ V s-1 ) ≤ -log ρ(m) 0 1 + G (m) 2 + 1 ( √ V s - √ V s-1 ) .
A telescoping sum argument yields

t s=1 η (m) s-1 L (m) s 2 ≤ (2 + G (m) ) -log ρ(m) 0   1 + t s=1 L (m) s 2 -1   ≤ (2 + G (m) ) -log ρ(m) 0 t s=1 L (m) s 2 .
Finally we get

- M m=1 π (m) t s=1 L (m) s ≤ M m=1 π (m)   t s=1 L (m) s 2 (3 + G (m) ) -log ρ(m) 0 +r (m) t + -log ρ(m) 0 + r (m) t
and the desired bound on R A t follows. In order to obtain the regret bound on R S t we use the Young inequality 2

√ ab ≤ γa + b/γ with γ = 4(2G (m) ∨ D 2 ) so that - M m=1 π (m) t s=1 L (m) s ≤ M m=1 π (m) 1 8(2G (m) ∨ D 2 ) t s=1 L (m) s 2 + 8(2G (m) ∨ D 2 )(3 + G (m) ) -log ρ(m) 0 + r (m) t 2 + -log ρ(m) 0 + r (m) t
and the desired result follows from an application of Theorem 3.1.

Discussion and examples

KAO algorithms require the knowledge of the variances σ 2(m) > 0. A natural estimator of this quantity is the mean square residuals

σ 2(m) t = 1 t t s=1 y s -ŷ(m) s 2 .
It can be tuned online but without any guarantee on the regret of the corresponding algorithm.

In our applications, we prefer to estimate σ 2(m) t on a burn-in period and use this fixed value in KAO.

Comparison with BOA

We start this Section with a short comparison with the BOA algorithm of [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF] that achieves similar regret bounds than the one obtained here. Theorem 4.2 in [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF] shows that BOA, which is an algorithm based on surrogate losses, has nice generalization properties that extend the regret bounds in the adversarial setting into similar regret bounds in the stochastic adversarial setting. The price to pay for the generalization is a factor 2 in the regret bounds. We show that this factor 2 is avoidable under assumption (H) with an algorithm such as KAO which uses the surrogate risks rather than the surrogate losses. Finally, notice that the use of the risk allows getting a.s. regret bounds in the well-specified stochastic unbounded setting rather than high-probability regret bounds only in bounded settings.

The static iid setting

In the iid setting, we consider an aggregation of static Kalman recursions with P (m) 0

= 1/λ (m) I, λ (m) > 0 which coincides with online ridge regression starting at θ(m) 0 . A natural estimator for σ 2 is the mean of the mean square residuals

M -1 M m=1 σ(m) t
. This setting is very specific

since µ t = E t-1 [y t ] = E[y t ] = X t θ (m) t = X t θ (m) = X t θ * for any 1 ≤ m ≤ M under (
L s (ŷ s ) ≤ min 1≤m≤M t s=1 L s (ŷ (m) s ) + O(-D 2 log ρ(m) 0 )) .
It is interesting to combine this bound with the regret bounds on the ridge regression when the design (X t ) is iid, bounded by X and such that E[X t X t ] has a positive lowest eigenvalue Λ min . Applying Theorem 14 of de Vilmarest and Wintenberger (2020), the mth Kalman recursion achieves for any θ ∈ R d and any t ≥ 1

t s=1 L s (ŷ (m) s ) ≤ t s=1 L s (X s θ) + O λ (m) 3 θ - θ(m) 0 6 2 + d log t λ (m) + log(δ -1 ) 3 ,
with probability at least 1 -δ. Moreover, the localization strategy of de Vilmarest and [START_REF] De Vilmarest | Stochastic online optimization using kalman recursion[END_REF] shows that under the same probability D can be considered as a constant. Then KAO achieves the regret bound in expectation, valid for any θ ∈ R d and any 1

≤ m ≤ M , t s=1 L s (ŷ s ) ≤ t s=1 L s (X s θ) + O λ (m) 3 θ - θ(m) 0 6 2 + d log t λ (m) -D 2 log ρ(m) 0 + log(δ -1 ) 3 ,
with probability 1 -M δ. Aggregation can be seen as an online alternative of cross-validation for tuning the starting point of the ridge regression algorithm and the regularization parameter.

As an illustration one should consider θ(m) 0 may be taken equal to α(e i ) 1≤i≤d where (e i ) 1≤i≤d

is the canonical basis and α takes value on [-d, d] ∩ Z. Moreover λ (m) should be taken on an exponential d finite grid of (0, ∞). The number of Kalman recursions is M = O(d) and choosing uniform weights yields to a regret for any λ > 0 on the grid, any 1

≤ i ≤ d and any α ∈ [-d, d] ∩ Z as t s=1 L s (ŷ s ) ≤ t s=1 L s (X s θ) + O λ 3 θ -αe i 6 2 + d log t λ + D 2 log d + log(δ -1 ) 3 ,
with probability 1 -dδ. Other aggregation strategies on least-squares estimators are described in [START_REF] Leung | Information theory and mixing least-squares regressions[END_REF]. Restrictions of our framework are the well-specification condition (H)

and the presence of the large constant D 2 in the model selection bound. One clear advantage is an explicit online procedure whereas least square estimators require the inversion of inverse matrices at each batch step.

The dynamic setting

In the dynamic setting, we consider that (y t ) behaves as a centered random walk conditionally on the design. The Kalman recursions track the trajectory of the linear coefficients associated to the explanatory variables. Assume that the design is standardized such that E[X

(m) t 2 ] = E[X (m ) t 2 ] for any 1 ≤ m, m ≤ d. Consider M = d univariate Kalman recursions d m = 1 with K (m) = Q (m) = 1. Then the random coefficients θ (m) t
satisfies the relation ( 6) and constitutes a random walk. If there exist D, X > 0 satisfying D = max

t≥1 |ŷ t -µ t | and |X t,m | ≤ X then one can bound, with high probability max 1≤t≤T max 1≤m≤d |L (m) t | ≤ 2DX max 1≤t≤T d m=1 | θ(m) t | ≤ CDX   T t=1 d m=1 E[( θ(m) t ) 2 ] 1/2 2   1/2
for some high constant C > 0. Then we can apply the result of [START_REF] Guo | Stability of recursive stochastic tracking algorithms[END_REF] asserting that

E[( θ(m) t -θ (m) t
) 2 ] 1/2 ≤ E for some E > 0. Together with the fact that Var(θ

(m) t ) = t by definition we obtain max 1≤t≤T max 1≤m≤d |L (m) t | ≤ CDXdT .
Then applying KAO with adaptive learning rate and doubling trick as in Remark 4.1 with G = CDXdT , we obtain with high probability the aggregation regret bound

T s=1 L s (ŷ s ) ≤ T s=1 L s M i=1 π (m) ŷ(m) s + O(DXT 3/2 d log d) .
This super-linear rate is due to the high fluctuations of the Kalman recursions when they track random walks (θ

(m) t
). The Kalman recursions inherit the high variability of the random walks which is responsible for the high variability of the gradient and large G = O(T ). However, due to the unboundedness of the response, none of the existing regret bounds seem to apply in this setting.

The expert aggregation setting

The setting is similar to the previous one as K (m) is a diagonal matrix with non-null coefficients equals to 1. Thus one has to assume the boundedness of the gradients to get a √ T regret for the aggregation problem. It is the usual assumption in the setting of aggregation of experts and then the regret is essentially divided by a factor 2 compared with the regret bound obtained for BOA in [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF] under the boundedness of the response. It is worth mentioning again that the boundedness of the gradients of the conditional risk does not imply the boundedness of the response.

Simulation study

In this simulation study, we use some of the variables contained in the downloadable data set on the website of the RTE company (french TSO) that describes the hourly electricity consumption and production per type of production units in France from 2013 to 2017. We chose to simulate synthetic data from these true ones to be closer to a real application but controlling the true model at the same time. We generate synthetic data from a subset of these variables: the temperature, the gas production, the fuel production, the charcoal production, and the nebulosity. The square of the temperature and the cubic of the gas are jointly utilized as predictors in X t to simulate, under a state-space model, the signal y t that represents the electricity consumption. All the covariates are normalized to be in [0, 1] by dividing each of them by their maximum value. The true model (that generates the true or the best expert) is a state-space model using the square of the temperature and the cubic of the gas as covariates in X t and Gaussian noise. Regarding the parameters of this state-space model, σ = 1.5, Q is of values 1 on the diagonal and 0.9 otherwise, θ 0 is generated according to a gaussian law of mean 500 and covariance matrix identity, and K is the identity matrix. We also compute 27 other Kalman experts using other combinations of covariates that are different from those used for getting the true (or best) expert. Each Kalman expert is computed in the sequential way as follows. We begin by fitting the model using the first observations (y 1 , . . . , y window ) contained in a window. Then the fitted model is utilized to predict the observations contained in a window ahead (i.e., y window+1 , . . . , y 2window ).

At the pth step we use the observations y 1 , . . . , y pwindow to fit the model that is utilized to predict y pwindow+1 , . . . , y (p+1)window . We chose window = 500 as a good trade-off between a correct number of observations to estimate the state-space models and a good adaptation to changes. The prediction resulting from this procedure is called the Kalman expert and we, therefore, have 28 Kalman experts.

Simulations are done under the R software (R Core Team, 2019) and the predictive performance of the Kalman experts aggregation using KAO is compared with the aggregation performed using the R package opera [START_REF] Gaillard | opera: Online Prediction by Expert Aggregation[END_REF] and the aggregation procedures therein. The aggregation obtained from the package opera is named OPERA when we are competing with the best expert and do not want to mention any specific aggregation procedure. We make one hour ahead prediction using KAO on the 28 Kalman experts. In the case where the oracle is the best Kalman expert, the resulting prediction is plotted by a red curve in Figure 1 at the left panel, where the signal y t is plotted by a black curve, and the experts are plotted using the gray color.

We can see that the red line tracks well the black one, meaning that the aggregation from KAO performs well its prediction. More precisely, the MSE of KAO is 66.507 which is approximately equal to the MSE of the best Kalman expert (66.503), and the MSE of OPERA is 253.06. The right panel of Figure 1 shows the cumulated error of KAO (in blue color) and OPERA (in red color).

We can see that KAO performs better than OPERA. Though both KAO and OPERA (precisely, EWA or BOA procedure) are based on exponential weights, the difference seen in their respective cumulated errors can be explained by the fact that KAO takes into account the underlying models that provide the experts, and OPERA doesn't have this information.

In the case where the oracle is the best convex combination of the Kalman experts, the one hour ahead predictions of y t , using KAO, are plotted in the left panel of Figure 2 in red color and the Kalman experts are plotted in gray color. We can also see that KAO tracks well the signal y t that is plotted in black color. Here, the MSE of KAO is 65.02 against 223.37 for OPERA, using the procedure BOA [START_REF] Wintenberger | Optimal learning with bernstein online aggregation[END_REF]. The corresponding cumulated errors are plotted in the right panel for KAO (in blue color) and OPERA (in red color). The curves of the cumulated errors show that KAO has a better predictive performance than OPERA.

We simulate 100 collections of Kalman experts corresponding to 100 simulated datasets. Each We can see that KAO (blue curve) has not only the best MSE but also presents more stability than BOA (red curve) and MLpoly (dashed green curve). Here, reversely to the case where the aggregation competes with the best expert, BOA is better than MLpoly. This simulation study seems to point out that it may be worth of interest to take into account the underlying model that generates the experts when aggregating them. is used for computing the loss and the pseudo-loss that are needed to perform KAO. The aggregation performance of KAO (on these experts) is compared with that of both MLpoly and BOA that are two aggregation procedures available in the opera package. The results are contained in Table 1 where GT means Gradient Trick. GT, therefore, refers to the case where the oracle of the aggregation procedure is the experts' best convex combination. For confidentiality reasons, errors are expressed relatively to the RMSE of the best convex combination.

The uniform procedure is the experts mean, and the procedure best convex is indeed the experts' best convex combination. All of these procedures are performed on the corrected experts.

We clearly see that KAO performs slightly better (rmse = 1.05 with GT and rmse = 1.07 without GT) than both MLpoly (rmse = 1.06 with GT and rmse = 1.16 without GT) and BOA (rmse = 1.07 with GT and rmse = 1.11 without GT). In order to check if the Kalman correction is worth of interest, we make a direct autoregressive correction of the experts that are then aggregated using MLpoly and BOA (not KAO as we need an estimate of the risk for that). The results are contained in Table 2 and show that all the procedures are less accurate when the Kalman correction is not applied.

The weights that are assigned to the corrected experts (by KAO and MLpoly) are plotted in ). This smoothness of KAO weights can be explained by the fact that the procedure uses the underlying properties of the model that provide the experts. This information is used to anticipate the forthcoming performance of each expert.

Conclusion

In this paper, we show that the prediction obtained by aggregating the predictions coming from a finite set of experts can be improved by taking into account the properties of the underlying models that provide the experts' prediction. We place ourselves in the case where all the predictions provided by the experts come from fitting state-space models using Kalman recursions. By using exponential weights, two settings are considered: 1) the aggregation competes with the best expert (also considered as model selection), and 2) the aggregation competes with the best convex combination of the experts. We consider adaptive multiple learning rates in order to achieve the optimal rates in these two schemes for a unique procedure. The quality of the aggregation's prediction has been improved by taking advantage of the full knowledge of the Kalman experts, using their predictive risk in an unbounded well-specified setting. In the simulations studies, we notice a great recovery of stability of KAO (our aggregation procedure), where all other existing aggregation procedures may be sometime somewhat unstable, potentially due to the lack of boundedness of the responses. The aggregation procedure KAO is also applied to some existing

∈

  R dm . Doing so, each Kalman recursion holds in a state space of dimension d (m) < d, lowering the computational cost of each iteration to O(d 2 m ).

  The variances σ 2(m) , 1 ≤ m ≤ M and the learning rate η. Initialization: The initial weights ρ m ≤ M . For each iteration t = 1, . . . , T : Inputs: The Kalman predictions ŷ(m) t+1 and the matrices P

Algorithm 3

 3 KAO for aggregation Parameters: The variances σ 2(m) , 1 ≤ m ≤ M and the learning rate η. Initialization: The initial weights ρ (m) 0 , 1 ≤ m ≤ M . For each iteration t = 0, . . . , T : Inputs: The Kalman predictions ŷ(m) t+1 and the matrices P (m) t , 1 ≤ m ≤ M . Recursion: Do:

  In our setting, the adaptation of the gradient-based EWA of Cesa-Bianchi and Lugosi (2006) yields Algorithm 3. The following theorem derives an upper bound for the regret R A t . Theorem 3.3. Under Assumption (H), suppose it exists G > 0 such that |L

Theorem 4. 1 .

 1 Under assumption (H) suppose it exists G > 0 such that |L (m) t | ≤ G a.s. for 1 ≤ t ≤ T , 1 ≤ m ≤ M .Then the aggregation regret of KAO with multiple learning rates

  H)and some fixed θ * ∈ Θ corresponds to the well specified setting. Consider for a moment D random. Moreover one can estimate G (m) = 4D 2 such that, applying KAO with adaptive multiple learning rates we obtain the model selection regret bound t s=1

Figure 1 :Figure 2 :

 12 Figure1: One hour ahead prediction of yt using KAO, and cumulated prediction errors for KAO and OPERA. The left panel shows one hour ahead prediction of yt using KAO and η within a grid of values. The value of η that minimizes the MSE is utilized to perform the prediction. These predictions are done in the case where the oracle is the best expert. The right panel shows the cumulated prediction errors for KAO and OPERA using the Kalman experts. the blue line represents the Kalman aggregation error, and the red one represents the error of the aggregation coming from the opera package. These predictions are done in the case where the oracle is the best expert.

  collection of Kalman experts contains 28 different experts. We then perform the aggregation of each collection of Kalman experts using KAO and the procedures within OPERA for each type of oracle. The MSE of the aggregations are computed and plotted in Figure 3. The left panel (Figure 3(a)) shows the curves of the MSE of the aggregations performed in the case where the oracle is the best expert. KAO (dashed blue curve) presents the lowest MSE within all the aggregation procedures, followed by MLpoly and EWA. The right panel (Figure 3(b)) shows the aggregations' MSE in the case where the oracle is the best convex combination of the experts.

Figure 3 :

 3 Figure3: MSE of 100 aggregations of Kalman experts using KAO and the procedures contained in opera. The left panel shows the curves of computed mse, where each aggregation procedure competes with the best expert. KAO is plotted in blue color, BOA in red color, EWA (Exponentially Weighted Average ) in black and MLpoly(Gaillard et al., 2014) in green color. The right panel shows the curves of the mse computed when the aggregation competes with the best convex combination of experts. KAO is in blue color, BOA in red color and MLpoly in green color.

  Q(m) , 1 ≤ m ≤ M and the variance of the noise σ 2(m) are estimated using an EM algorithm on the first half of the data (t ∈ (1, 2, ..., T /2)) and we use the second half to evaluate KAO performances and compare it to other aggregation rules.

Figure 4 :

 4 Figure 4: Experts weight according to KAO and MLpoly, using the gradient trick. The left panel shows the weights assigned to the corrected experts by KAO where the oracle is the best experts combination. The right panel shows the weights assigned by MLpoly, using the gradient trick. The experts are denoted by E1, . . . , E10.

Figure 4 .

 4 Figure 4. The weights coming from KAO are more smooth (see Figure 4(a)) than those provided by MLpoly (see Figure 4(b)). This smoothness of KAO weights can be explained by the fact

Table 1 : Root Mean Square Error square of different aggregation procedures (relative to RMSE of the best convex combination). Kalman Experts.

 1

	Procedure	rmse (with GT)	rmse (without GT)
	Best expert	1.15	1.15
	Uniform	1.11	1.11
	MLpoly	1.06	1.16
	BOA	1.07	1.11
	KAO	1.05	1.07
	Best convex	1	1
	7. Application

In this section we apply the KAO algorithm to aggregate ten experts f m,t , 1 ≤ m ≤ M that are meant to predict the daily electricity consumption in France (see

[START_REF] Ba | Adaptive learning of smoothing functions: Application to electricity load forecasting[END_REF]

,

[START_REF] Gaillard | Forecasting electricity consumption by aggregating experts; how to design a good set of experts[END_REF]

for previous work on french load data) at times t ∈ (1, 2, ..., T ). These experts are

Table 2 : Root Mean Square Error square of different aggregation procedures (relative to RMSE of the best convex combination). AR experts.

 2 provided by different models that are black boxes. Thus, we consider the expert setting previously defined in 2.2. For each expert we stack their predictions f m,t ∈ R in X

	Procedure	rmse (with GT)	rmse (without GT)
	Best expert	1.18	1.18
	Uniform	1.11	1.11
	MLpoly	1.07	1.19
	BOA	1.07	1.09
	Best convex	1	1
			(m) t	together with the

intercept and the past error e m,t-1 = (y t-1 -f m,t-1 ), i.e.,

X (m) t = (1, f m,t , e m,t-1 ), t ≥ 1 .

and each state-space model m is defined by the state equation:

experts coming from unknown models, where we suggest correcting the errors of the experts using Kalman recursions. This strategy allows for approximating the theoretical weights needed for KAO and shows a quite important increase in the accuracy of the aggregation. In the case where the errors of the experts show no stationary behavior (for example, when there exist some cluster of variance), it should be interesting to adapt the fitting of the underlying state-space model in order to remain accurate.


E Ferrero

A Longhetto

L Montabone

L Mortarini

) M Manfrin

Physical simulations of neutral boundary layer in rotating tank

Keywords: 92.60.Fm -Boundary layer structure and processes. 92.10.Lq -Turbulence and diffusion. 92.60.Ek -Convection, turbulence, and diffusion

Laboratory measurements of neutral atmospheric boundary layer (ABL) are presented. The experiment was carried out in the large rotating tank of the Coriolis-LEGI laboratory in Grenoble. An ABL was created at reduced scale and measured. The mean flow was generated, both increasing (spin-up) and decreasing (spin-down) the rotation speed of the platform. Preliminary measurements by acoustic probes were used to assess the decay with time of the flow velocity at the position where turbulence measurements were subsequently performed. The mean velocity and the turbulence fields were then measured by using PIV (Particle Image Velocimetry) technique, which allows to obtain high-resolution measurements in the simulated ABL. For each rotation period and velocity conditions, two vertical cross-sections of the flow were measured at different times, then 3D velocity fields were reconstructed from the two planar fields under convenient geometric (orthogonal) and physical assumptions (reproducibility of the flow). The aim of this work is to collect a useful data set for testing and comparing turbulence models and parameterisations. For this reason, particular attention was paid to the turbulence statistics, turbulent fluxes and scales. The results of the data analysis are presented and discussed.

-Introduction

It is well known today that in order to understand how irregularities of the ground surface can distort the mean and turbulent structure of the flow and thus to obtain a com-plete picture of its structure, it is necessary to solve a full set of fluid dynamics equations using finite differences or more sophisticated numerical methods (see for instance [1,2]).

The fundamental role played by the higher order turbulence moments in reproducing the turbulent transport not only in complex terrain and convective conditions but also in the neutral horizontally homogeneous conditions, has been recently demonstrated by [2].

On the other hand, the characteristics of the simulated turbulence depend very sensitively on the closure scheme which is adopted. Simple schemes, such as mixing-length model, which are essentially local in nature, are adequate only if there is a near equilibrium between local production and dissipation rates of turbulent kinetic energy. Otherwise non-equilibrium effects become significant and thus a realistic description of the turbulence structure requires considerably more sophisticated closure [3]. A simple model in which the eddy exchange coefficient does not take into account horizontal variations would not be able, in principle, to simulate flows over irregular terrain. Only the introduction of transport equations allows to account for such inhomogeneities.

For these reasons, the assessment of turbulent fluxes of momentum and wind field in the atmospheric boundary layer is a topic of great interest from both a theoretical and practical point of view and still an open problem.

One of the main problems in developing turbulence models is to find a complete dataset for determining parameterisations and empirical constants and for comparing the model results. These data-sets are difficult to obtain in open-field experiments because of the uncontrolled flow conditions. Useful data can be gathered in wind tunnel experiments, but these cannot take into account the Earth rotation [4,5]. Recently LES was adopted in literature as a base for comparison of other models and parameterisations [6,7].

In this work a method for obtaining detailed observations of a rotating neutral boundary layer from laboratory experiments carried out under controlled conditions is presented. In order to assess the reliability of the method, a few simple schematic examples of ABLs over both smooth and rough flat surfaces have been selected. The laboratory facilities and measurement methodologies are illustrated and the results in terms of turbulence moments and diffusion coefficient are shown.

-The laboratory experiment

This work consists in a series of laboratory simulations of the dynamics of a rotating neutral turbulent ABL and their influence on the diffusive conditions of the atmosphere.

The laboratory experiments took place on the hydrodynamic rotating tank of "Coriolis/LEGI" (Grenoble). The rotating speed of the tank can be continuously changed up to 5 revolutions per minute, allowing to achieve Rossby numbers typical of mesoscale circulations (10 -1 < Ro < 1). The large platform of "Coriolis/LEGI" laboratory is essential for getting the above-mentioned results not only on account of the horizontal spatial scale, but mainly because it allows the necessary resolution for reproducing the surface roughness elements and for determining the velocity field inside the topographical features of the simulated complex terrain.

The LASER-PIV (Particle Image Velocimetry) system, annexed to this hydrodynamic simulation laboratory, allows the quantitative measurement of fluctuating velocity fields. This can be accomplished by seeding the water with particles (60 µm in diameter), illuminating it with a light sheet produced by a continuous Argon laser equipped with a beam expanding optical device. A high-resolution numerical CCD camera records sequences of images shot at adjustable time rates and a dedicated software identifies the displacement of clusters of particles by means of cross-correlation technique between pairs of images. The velocity of the flow at any point of the measured region [8] is then calculated at each time.

A schematic picture of the experiment layout is shown in fig. 1. Two different configurations of the camera and the light sheet were used to get two horizontal components of the flow. These components were subsequently composed to reconstract the full 3D velocity field.

The experiments were carried out with different values of the roughness length. In the first case no roughness elements were put in the tank, while in the second series of experiments cubic-shaped elements were glued to the tank bottom.

The mean flow U (r) at the distance r from the center of the tank was created either by increasing (spin-up) or by decreasing (spin-down) the tank rotation speed. The initial velocity of the flow was controlled by choosing different values for the initial (T 0 ) and final (T 1 ) rotation periods.

The turbulent equilibrium boundary layer depth can be evaluated as

(1)

δ t = 0.4 u * f ,
where f = 4π/T 1 .

The frictional velocity was evaluated by the following empirical formula [START_REF] Ferrero | Proceedings of Physmod2003[END_REF]:

(2) u * = αu ∞ ,
where u ∞ is the free stream velocity outside the turbulent layer and the coefficient α = 0.07 was estimated, in a previous work [START_REF] Ferrero | Proceedings of Physmod2003[END_REF], from the velocity measurements performed with a total water depth of 30 cm and different spin-up and spin-down conditions. It can be noted that the free stream velocity in the tank plays the role of the geostrophic wind in the actual atmosphere.

Considering that U = 2π(1/T 1 -1/T 0 )r, we have

(3) δ t = 0.2α 1 - T 1 T 0 r
which does not depend on the rotation period. Since in all the experiments we have T 1 /T 0 = 1/2 and r = 5.1 m, we expect (4) δ t = 3.6 cm .

On the other hand, the Reynolds number based on the laminar depth is

(5) Re l = Uδ ν = 1 - T 1 T 0 r (2π) 1/2 (νT 1 ) -1/2 = 6390 T -1/2 1
, where

(6) δ = νT 1 2π 1/2 .
It may be useful to recall that the Ekman layer instability is attained for Re l of the order of 55. Hence, turbulence is fully developed for Reynolds numbers higher than this value.

In the experiments with values of T 1 ranging from 30 s to 240 s we obtained Reynolds numbers ranging from 412 to 1166. These values refer to the initial velocity; during the experiment the velocity decreases, due to the viscous forces acting on the system, and the Reynolds numbers lower. The flow, initially turbulent, after about 10 2 normalised times becomes laminar, as can be seen in fig. 2, where the normalised velocity U/U 0 (U 0 being the initial velocity after the spin-up or spin-down), measured by acoustic probes in the bulk of the flow, is depicted as a function of the normalised time, tU 0 /H (where H = 30 cm is the total depth of the fluid in the preliminary tests). It can be observed that the slopes of the curves change at this time indicating the different regimes of the two intervals. However, a turbulent boundary layer was generated close to the wall also for normalised times greater than 10 2 .

-Data preprocessing and processing

For each view of the flow, as explained above, the experiments have produced a time series of 714 raw images of the fluid uniformly seeded with particles. Cross-correlation between a pair of images is used to obtain the displacements of patches of particles over a sufficiently small time interval, a technique known as Particle Image Velocimetry (PIV) or Digital Particle Image Velocimetry (DPIV). See [START_REF] Raffel | Particle Image Velocimetry, a Practical Guide[END_REF] for a complete review of the bases of this technique, although it does not include recent developments and improvements. Before using PIV, raw images have been preprocessed to eliminate a very bright line which corresponds to the illumination of the particles settled down on the tank floor. We performed an average (pixel by pixel) over the whole series of 714 images and then removed, from each image, only the part close to the floor with a weighed subtraction.

The Correlation Image Velocimetry (CIV) algorithm by [START_REF] Fincham | [END_REF] is used to obtain the 2D velocity vectors on a planar grid. This state-of-the-art algorithm incorporates advanced features (such as the multi-passage option) which improve the performance and allow for obtaining accurate values for the velocities; see [8,[START_REF] Fincham | [END_REF] for details. Two versions of CIV (CIV1 and CIV2) have been developed. The more advanced algorithm corresponds to CIV2, while CIV1 is rather standard. For the purpose of these experiments, we have mainly used the on-line correlation peak averaging of Hart for reducing the number of wrong vectors, together with a post-processing outliers detection and filtering algorithm based on the low correlation value of an interrogate vector. The quality of most of the raw images after the preprocessing was such to use a rectangular interrogation window of 30 × 20 pixels, obtaining a final resolution of about 60 × 80 vectors on the grid (the exact numbers vary slightly in the different experiments, due to different camera view fields).

Each experiment includes two time series of 714 raw images of the fluid taken at an angle of 45 degrees with respect to the direction of the flow, as explained above. These time series span about 12 seconds at a frequency of 60 Hz, and for most experiments consecutive images have been used for cross-correlation, so that velocity fields are available at the same frequency.

-PIV measurements analysis

As said before, for each flow condition we measured two planar fields on two vertical orthogonal cross-sections. From them, we reconstruct 3D flow fields, assuming statistical reproductibility of the flow. The streamwise (azimuthal) and crosswise (radial) velocity components, u and v, respectively, are obtained from the horizontal velocity components u 1 and u 2 , measured with each camera

u = ( √ 2/2)(u 1 + u 2 ) , (7) v = ( √ 2/2)(u 1 -u 2 ) .
The wall is visible in each image, and the apparent vertical is defined as the direction perpendicular to this wall (it is nearly aligned with the image y-axis). In reality the apparent vertical particle displacement is a projection in the field of view, which is slightly tilted with respect to the true vertical (see fig. 3). The local tilt angle at height z from the wall is (zz c )/L, where z c = 15 cm is the camera height above the tank floor and L = 200 cm the distance between the camera and the field of view. Using the approximation of a small angle, the true vertical velocity w is related to the apparent vertical velocities w 1 and w 2 seen by the camera 1 and 2, respectively, by the geometric projections:

w = w 1 -(z -z c )u 2 /L , ( 8
)
w = w 2 + (z -z c )u 1 /L .
The time-averaged velocity components u , v , w are obtained by applying the time averaging operator to (7) and ( 8) taken at the same height z. The mean vertical velocity w is then obtained twice by each of the cameras, and we take as the final result the average of the two values.

Multiplying the two relations (8) by u 1 and u 2 respectively before averaging, we get Adding these two relations, we get the streamwise (azimuthal) Reynolds stress (the parallax corrections in u 1 u 2 cancel each other) [START_REF] Raffel | Particle Image Velocimetry, a Practical Guide[END_REF] wu = (

wu 1 = w 1 u 1 -(z -z c ) u 1 u 2 /L , ( 9
)
wu 2 = w 2 u 2 + (z -z c ) u 1 u 2 /L .
√ 2/2)( w 1 u 1 + w 2 u 2 ) .
Note that this method does not apply for the transverse Reynolds stress wv : the parallax corrections in u 1 u 2 would then add up, and the correlation u 1 u 2 cannot be obtained by our method, since the two components u 1 and u 2 are not measured simultaneously.

In fig. 4 an example of the power spectrum obtained in the case of the experiment with U = 4.6 cm/s and no roughness, at the level n = 10 is shown (inside the boundary layer). The turbulence layer appears to be well developed and a two-decades inertial sub-range regime can be observed. This result allows to analyse the boundary layer and to compare it with the case of the atmosphere.

In order to plot the vertical profiles in non-dimensional coordinates and comparing the different experiments, we should evaluate the boundary layer height and the friction velocity.

The boundary layer height z i can be deduced as the level where the turbulent kinetic energy (TKE) vanishes. The TKE is defined as (here

TKE = 1 2 u 2 + v 2 + w 2 (11) 1 2 u 2 1 + u 2 1 + 1 2 ( w 2 1 + w 2 2 )
u i = u i -u i ).
For this quantity, we can neglect the parallax corrections ( 8) to an excellent approximation and for w 2 we just take the average from the two cameras.

In fig. 5 the vertical TKE profiles corresponding to the six experiments considered are depicted. The boundary layer height z i was evaluated from these profiles considering the level at which the profiles attain zero. From the vertical profiles the free stream velocity was also estimated for the six experiments and the friction velocities were finally determined from (2). In table I the estimated parameters for each experiment are sumarised.

Figures 6 and7 show the vertical profiles of u and v velocity components made non-dimensional by using the above-estimated parameters. It can be observed that, as expected, the curves corresponding to different runs tend to collapse on the same profile, demonstrating the correctness of the estimation of the surface layer parameters.

We can conclude that the non-dimensional parameters, u * and z i , calculated as described above, are the characteristic parameters for the simulated flows.

Concerning the w components (fig. 8), it can be noted that they were not equal to zero, as it should be over a flat boundary; this could be due to the Ekman pumping generation by the tank rotation.

For the sake of comparison, together with the experimental results we also reported the profiles predicted by two Large Eddy Simulations (LESs).

LES model [7] utilises the incompressible Boussinesq form of the Navier-Stokes equations and considers a horizontally homogeneous boundary layer. The ABL variables (u, v, w, pressure, temperature, etc. ...) are spatially filtered to define resolved component and sub-grid scale (SGS) component. The filter width is of order of the numerical grid dimension and lies in the inertial subrange of scales in the turbulence spectrum. In the model [6] the SGS energy is found from a prognostic equation [12], while the SGS fluxes are obtained from a first-order closure or flux-gradient relationship. The stress tensor is given by ( 12)

τ ij = -2ν t S ij ,
where the deformation tensor S ij is defined as

(13) S ij = 1 2 ∂u i ∂x j + ∂u j ∂x i .
The eddy viscosity ν t is parameterized as ( 14)

ν t = C k le 1/2 s ,
where C k is a constant ( 0.1), l is a mixing length, e s is the SGS TKE and its dissipation rate is given by In this work we reproduce by the LES two neutral ABL utilising the same procedure already described by [6], that is first generating a mixed (sheared and convective) ABL, then setting surface heat flux to zero. They are characterised by (LESFE1) U ∞ = 15 m/s, u = 0.5 m/s, z i = 478 m and (LESSQ2) U ∞ = 15 m/s, u = 0.58 m/s, z i = 1574 m, respectively. Figure 9 shows that the experimental results underestimate the simulated values of about a factor two.

In fig. 10 the vertical profile of the u w momentum fluxes is shown. Also in this case the experimental results underestimate the LES data, with a different degree as greater as the velocity is smaller.

-Diffusion coefficients

When a dispersion process is considered over large vertical scales covering not only the inner and outer ABL but also the lower free atmosphere over an extended horizontal fetch, numerical atmospheric models incorporating turbulent diffusion processes are required in order to describe the diffusion and transport of airborne tracers ( [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF], pp. 386-389). This means that it is often necessary to be able to approximate the surface stress and fluxes in terms of mean variables at the grid points in the above-mentioned numerical models. Unfortunately, in some models the lowest grid point is well above the surface layer, making it impossible to use flux-profile relationships based on the non-local, Surface Layer Monin-Obukov similarity theory. However, a different similarity theory, the Surface Rossby Number Similarity, built on the matching of surface layer profiles to the external forcing given by flows in the mid boundary layer, can permit that surface stress and fluxes be related to conditions higher in the ABL in terms of mean variables at their grid points ( [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF], p. 359), [START_REF] Tennekes | Workshop on Micrometeorology[END_REF]15].

It is then possible to match the flow at the bottom of the mid-ABL with that at the top of the surface layer and relate the surface stress to the mean flow higher in the interior of the ABL, giving the bulk transfer, or drag coefficients, which allow one to find the surface flux given knowledge of mean variables aloft and at the surface ( [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF], p. 387). In fact, when we have the bulk transfer, or drag coefficient [CD = V 2 /u 2 * ], we can use the following relationships, parameterised by many authors [START_REF] Clarke | The Wangara Experiment: Boundary Layer Data[END_REF][17][18][19] to obtain surface fluxes Φ of a property ( [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF], p. 270):

(16) Φ = -C g V g [ top -bottom ] ,
where C g is called the geostrophic drag coefficient and V g is the magnitude of the geostrophic wind. Here the geostrophic wind and then the geostrophic drag coefficient are used, because sometimes in numerical or theoretical models one does not know the winds anywhere near the surface, but one can calculate the geostrophic wind instead [START_REF] Stull | An Introduction to Boundary Layer Meteorology[END_REF], p. 270.

In the particular event of the geostrophic drag coefficient C GN for neutral atmospheric conditions, the following parameterisation has been suggested: (17) C GN = 0.0123 Ro -0.14 0 , where Ro 0 = V g /f z 0 represents the surface Rossby number. This number differs from the Rossby number Ro based on the mesoscale scales for the lengths and velocities.

The main problem to deal with now is to look for the set of scales that produce the correct non-dimensionalization of the ABL equations, that is guarantee the Surface Rossby Number Similarity, assuring that everything becomes independent of the external parameters, in particular the non-dimensional vertical profiles of the most relevant turbulent parameters, such as the turbulent diffusion coefficient, the dissipation rate of turbulent kinetic energy, the vertical temperature gradient, and so on. As the geostrophic wind can be calculated quite easily, equations expressing the matching between the bottom and higher ABLs through this scale velocity, also known as geostrophic drag laws, have attracted particular interest in the community of ABL physicists. It has been demonstrated by [20] that the Surface Rossby-number similarity condition in the ABL exists above a height Z R = Z 0 , where Z R and Z 0 are the normalised (standardised) vertical height and roughness length, respectively, if the relevant physical quantities of the ABL are made dimensionless with the following internal scales: friction velocity u * (18) internal scale height h = ku * /f . (

It can be noted that h is the turbulent boundary layer depth δ t already defined (eq. ( 1)) and, in general may not coincide with the experimental value z i listed in table I. With these internal scales, it is possible to define the following reduced (nondimensionalized) variables with internal parameters:

Z = z h , ( 20
)
K m = k m h 2 f = k m ku 2 * /f , ( 21
) (U, V ) = k u, v u * , ( 22
) (T x , T y ) = (τ x , τ y ) ρu 2 * , ( 23
)
where K m is the reduced turbulent diffusion coefficient, (U, V ) are the reduced x and y wind speed components, (T x , T y ) are the reduced turbulent stress component and k is the von Karman constant. By using Prandtl's mixing length relation

(24) k m = 2 dv dz ,
where is the mixing length, a new expression can be obtained for the reduced K m . We obtain first (25)

K m = L 2 dV dZ 1 k 2 ,
where L = /h = f /(ku * ) is the reduced mixing length non-dimensionalized with internal parameters. Then, by recalling the flux-gradient relationship (τ

x , τ y )=ρk m d(u, v)/dz, one can obtain k m = l(τ /ρ) 1/2 or, in reduced form, ( 26
) K m = L k (T 2 x + T 2 y ) 1/4 .
If the orientation of the x-axis of the coordinate system is chosen along the direction of the surface stress τ 0 (rather than the surface geostrophic wind V g0 ) and reduced variables non-dimensionalized only with internal parameters are used, in [20] it is demonstrated that the equations of the ABL, as well as their boundary conditions, are independent of any external parameters for Z ≥ Z R (Surface Rossby Number Similarity). More precisely, they also demonstrated that this property, that holds true for almost all of the quantities defining the structure of a rotating turbulent ABL, does not apply for the profiles of K m . In fact, in the chosen coordinate system (where T x = T x0 , T y = 0), the non-dimensional mixing-length hypothesis reads as L 0 = kZ 0 at z = z 0 , and eq. (25) will become (again at z = z 0 ): K m0 = Z 0 and hence

K m0 = L 0 /k or, in dimensional form, k m0 = ku * z 0
which is the well-known Prandtl mixing-length hypothesis for z = z 0 , assuring a perfect coincidence of the different K m profiles in the lowest part of the ABL. This is not the same if we use reduced variables non-dimensionalized with external parameters, because in this case the equations of the ABL depend on external parameters which contain internal as well as external parameters. When conditions of Surface Rossby Number Similarity apply, reduced vertical profiles of T x (Z), T y (Z) and K m (Z) in barotropic ABLs are universal functions, valid for arbitrary external parameters. The experiments carried out in the Coriolis tank offered a special opportunity to test the validity of the approximation of the above law.

In figs. In order to make the vertical coordinate non-dimensional, we use the experimental ABL height z i .

As has been explained before, the different experiments have been carried out with two different roughness conditions. The lowest one corresponded to the roughness of the floor of the tank, which was quite smooth; the second one corresponded to roughness created on the tank floor by gluing on it a large number of small cubes whose sides were l c = 4 × 10 -3 m, almost randomly distributed with a mean distance of 8 × 10 -2 m in order to avoid regular and stationary streamlines at the bottom of the water flow. This distribution of cubes produced an equivalent roughness z 0 ∼ 10 -1 l c , corresponding to a dimensionless roughness Z 0 = 5 × 10 -3 , considering a turbulent boundary layer depth equal to 8 × 10 -2 m.

According to the previous theoretical development, the surface Rossby number similarity should exist everywhere in the ABL above a height Z g ≥ 10Z 0 , where the vertical profiles of all physical quantities should collapse on a unique, universal profile, except for the reduced turbulent diffusion coefficient K m when a relationship, expressing the mixing-length hypothesis, is introduced. In which case, the surface Rossby number similarity should exist for Z g ≥ Z 0 .

However, our normalised K m profiles seem to violate this second theoretical expectation, and to comply rather with the first and more general property, according to which surface-Rossby number similarity appears only above the reduced height Z R ≥ 10 Z 0 . As a matter of fact, the dimensionless profiles collapse into a unique profile in the layer between Z = 0.05 and Z = 0.1 (fig. 11b).

Concerning the cases without roughness, it can be observed that the two dimensionless profiles match at a lowest level of the range where the Rossby number similarity applies (Z = 0.05).

One explanation of this departure from the theoretical expectation could be that the turbulence developed near the bottom of the water layer in the rotating tank was not similar to the homogeneous, small-scale turbulence assumed by the mixing-length assumption for the lower layers of the atmospheric ABL. This explanation is suggested by the fact that considerable doubt has been long since thrown on the fundamental assumption that the shearing stress in a boundary layer is directly related to the gradient of mean velocity. On the other hand, the Prandtl-mixing length theory appears to be a "small-eddy theory" ([13, p. 208]) for turbulence generated mechanically and does not take proper account of the influence of the bigger eddies; as these large eddies are absent in the turbulent flow of the lower surface atmospheric boundary layer, some success attended the introduction of mixing-length concepts on account of the fact that momentum seems to be transferred mainly by smaller isotropic eddies (the same does not occur for heat transfer, for which some extraneous effect such as buoyancy enters). However, across the whole depth of a real rotating atmospheric ABL, bigger nonisotropic vortices can be present (as well as wakes due to complexity of terrain, or large 2D coherent structures due to the Earth rotation). This fact can then support the presence of some transport processes in turbulent flow that is independent of the gradient of the entity being transferred, i.e. there might exist some "convective-like" processes as well as the usual diffusion mechanism.

Moreover, measurements of wind profiles near the ground indicate that is of the same order of magnitude as the height of observation, and therefore is not "small in the free-path sense". This fact too has the consequence that it may not be possible always to formulate a consistent theory of turbulent properties by using simply the down gradient approximation: there may be a term, or terms, depending on the transported variable and some mean velocity [2,21].

In conclusion, our experiments highlighted that only the theoretical prediction of surface Rossby-number similarity in the atmospheric boundary layer is confirmed by experimental observation conducted in a rotating tank, while the accompanying Prandtl mixing length assumption is untenable.

-Conclusions

The experiment carried out in the rotating Laboratory of LEGI-Coriolis is described. The experiment was aimed to simulate a neutral boundary layer developed over flat boundary with different roughness lengths. Advanced measurement techniques, as PIV, were used in order to gather with high resolution data both in time and space. The results of the data analysis are presented in terms of mean flow and turbulence moments vertical profiles. The surface layer parameter and the boundary layer height were estimated in order to make the dynamical variable non-dimensional and to compare the different experiments with LES results. The experimental data underestimate the simulated ones. Nevertheless they show that the main characteristics of an ABL are capturated and, at least qualitatively reproduced. Concerning the diffusion coefficients, it has been shown that a partial surface Rossby number similarity exists on for Z R ≥ 10 Z 0 due to the untenability of the Prandtl mixing length when non-isotropic big vortices exist in the rotating turbulent flow.

These conclusions suggested that the methodologies presented can be useful for comparing new turbulence parameterisation or more sophisticated models and for assessing their ability to reproduce the main features of the boundary layer. This is particularly interesting due to the possibility of the experimental facility to account for rotation.

Fig. 1 .

 1 Fig. 1. -Experimental set-up.

Fig. 2 .

 2 Fig. 2. -Normalised velocity as a function of the normalised time.

Fig. 3 .

 3 Fig. 3. -Parallax correction.

Fig. 4 .

 4 Fig. 4. -Power spectrum (σ 2 w /f ) as a function of the frequencies (f ). The line represents the f -5/3 behaviour in the inertial subrange.

Fig. 5 .

 5 Fig. 5. -Turbulent kinetic energy vertical profile.

Fig. 7 .

 7 Fig. 7. -Normalised v vertical profiles.

Fig. 8 .

 8 Fig. 8. -Normalised w vertical profiles.

  comparison of the normalised measured and simulated TKE has been shown.

Fig. 9 .

 9 Fig. 9. -Normalised turbulent kinetic energy vertical profiles compared with LES data.

Fig. 10 .

 10 Fig. 10. -Normalised vertical profiles of the u w compared with LES data.

  11 a and b the momentum diffusion coefficient K m vertical profiles, (27) K m = -u w dU/dz for the different experiments are depicted in dimensional and non-dimensional coordinates.

Fig. 11 .

 11 Fig. 11. -Diffusion coefficient (a), normalised diffusion coefficient (b).

Table I .

 I -Estimated parameters

	Experiments	U (cm/s)	u (cm/s)	zi (cm)	T0	T1
	norug92	10.67	0.74	9.	240	120
	norug46	5.95	0.42	8.	240	120
	rug92	11.20	0.78	9.	240	120
	rug65	7.80	0.55	8.	240	120
	rug46	5.65	0.4	8.	240	120
	rug325	4.50	0.31	7.	240	120

Fig. 6. -Normalised u vertical profiles.

The authors wish to thank D. Bertoni, R. Forza and S. Viboud for their skilful technical assistance before and during the experimental runs. This work has been carried out in the frame of the EU Program "Transnational Access to major Research Infrastructures".