Anant Raj
email: raj.anant12@gmail.com

Pierre Gaillard
email: pierre.gaillard@inria.fr

Christophe Saad

Regression Non-Stationary Online

Non-stationary Online Regression

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

We consider online linear regression in a nonstationary environment. More formally, at each round t = 1, . . . , n, the learner receives an input x t ∈ R d , makes a prediction y t ∈ R and receives a noisy output y t = x t θ t + Z t where θ t ∈ R d is some unknown parameter and Z t are i.i.d. sub-Gaussian noise. We are interested in minimizing the expected cumulative error

R n ( y 1:n , θ 1:n ) := n t=1 E ( y t -x t θ t ) 2 .
(

) 1
Of course, without further assumption, the cumulative error is doomed to grow linearly in n. Therefore, we assume there is regularity in the signal θ 1:n = θ 1 , . . . , θ n ∈ R d×n , measured by its total variation

T V (θ 1:n ) = n t=2 θ t -θ t-1 1 . (2)
We also assume that there exists B > 0 such that for all t ≥ 1, θ t 1 ≤ B. We emphasize that apart from boundedness in 1 -norm and in total variation, we do not make any assumption on the sequence θ 1:n . The latter is arbitrary and may be chosen by an adversary.

Related Works Online prediction of arbitrary time-series has already been well studied by the online learning and optimization communities and we refer to the monographs [START_REF] Cesa | Prediction, learning, and games[END_REF][START_REF] Hazan | Introduction to online convex optimization[END_REF] and references therein for detailed overviews. A very large part of the existing work only deals with stationary environment, in which the learner's performance is compared with respect to some fixed strategy that does not evolve over time. Thanks to many applications (e.g. web marketing or electricity forecasting), designing strategies that adapt to a changing environment has recently drawn considerable attention.

Online learning in a non-stationary environment was referred under different names or settings as "shifting regret", "adaptive regret", "dynamic regret", or "tracking the best predictor" but most of these notions are strongly related. Some relevant works are [START_REF] Besbes | Non-stationary stochastic optimization[END_REF][START_REF] Bousquet | Tracking a small set of experts by mixing past posteriors[END_REF][START_REF] Cesa-Bianchi | Mirror descent meets fixed share (and feels no regret)[END_REF][START_REF] Hazan | Efficient learning algorithms for changing environments[END_REF][START_REF] Herbster | Tracking the best linear predictor[END_REF][START_REF] Jadbabaie | Online optimization: Competing with dynamic comparators[END_REF][START_REF] Mokhtari | Online optimization in dynamic environments: Improved regret rates for strongly convex problems[END_REF][START_REF] Yang | Tracking slowly moving clairvoyant: Optimal dynamic regret of online learning with true and noisy gradient[END_REF][START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF]. [START_REF] Herbster | Tracking the best linear predictor[END_REF] first considered shifting bounds for linear regression using projected mirror descent. [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF] provides dynamic regret guarantees for any convex losses for projected online gradient descent. Most of these work considered however non noisy observations (or gradients), as we consider. [START_REF] Besbes | Non-stationary stochastic optimization[END_REF] proved matching upper and lower bounds for the dynamic regret with noisy observations. They provide dynamic regret bounds of order T V (θ 1:n ) 1/3 n 2/3 for convex losses and T V (θ 1:n )n for strongly convex losses. The latter was generalized to exp-concave losses by [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF].

Contributions Most of the above works consider the regret, while here we consider the cumulative error [START_REF] Katy | Relative loss bounds for on-line density estimation with the exponential family of distributions[END_REF]. In other words, in our case, the performance of the player is only compared with respect to the true underlying sequence θ 1:n which must have low total variation. This assumption allows us to prove stronger guarantees. Indeed, in the one-dimensional setting of online forecasting of a time-series with square loss, [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] could prove that the optimal rate of order T V (θ 1:n ) 2/3 n 1/3 instead of T V (θ 1:n )n for the cumulative error [START_REF] Katy | Relative loss bounds for on-line density estimation with the exponential family of distributions[END_REF]. Their technique is based on change point detection via wavelets and heavily relies on their simple setting (one dimension, no input x t ).

In this work, we generalize the result of [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] to online linear regression in dimension d and to reproducing kernel Hilbert spaces (RKHS). We ended up by using the meta-procedures of [START_REF] Hazan | Adaptive algorithms for online decision problems[END_REF] and [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] for exp-concave loss functions, combined with well-chosen subroutines. Carrying a careful regret analysis in our setting, we achieve the optimal error of [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF].

Finally, in Section 4, we corroborate our theoretical results on numerical simulations.

2 Warm-up: Online Prediction of Non-Stationary Time Series

In this section, we discuss the relevant background to our work and simple intuition for 1-dimensional problem. However, before going into the details of our approach, we first discuss the work of Baby and Wang [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] which considers one dimensional non-stationary online linear regression.

ARROWS [2]

ARROWS considers to solve the problem of online forecasting of sequences of length n whose totalvariation (TV) is at most n. The observed output is the noise contaminated version of original input sequence θ i for i in [n]. ARROWS considers to predict via the moving average of the output in an interval. If the total-variation within that time interval is small then the moving average in that time interval is reasonably good prediction to minimize the cumulative squared error. For that reason, the algorithm needs to detect intervals which has low total variation. This task of detection is accomplished by constructing a lower bound of TV which acts like a threshold to restart the averaging and hence acts like a non-linearity which can capture the non-linear variation in the sequence. The estimation of the lower bound is based on computing of Haar coefficients as it smooths the adjacent regions of a signal and then taking difference between them. A slightly modified version of the soft threshold estimator from from Donoho et al. [START_REF] David L Donoho | Minimax risk over hyperrectangles, and implications[END_REF] is considered for oracle estimator.

Overall, the restart strategy based on change point detection using Haar coefficients proposed in this work achieves the optimal error however, the approach is very hard to extend beyond one dimensional regression problem. Another drawback this work has is that ARROWS requires to know the noise level sigma to tun the algorithm even in one dimensional forecasting problem. To know the exact noise level is an unrealistic assumption in real life problems. We address here these two concerns.

One-Dimensional Intuition

In this section, we consider the simpler case with d = 1 and x t = 1 that was already considered by [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] as a warm up to understand the intuition behind our algorithm. Let us now define the formal problem. The problem formulation looks as: y t = θ t + Z t for t = 1, • • • , n and Z t be independent σ sub-Gaussian random variables. The goal of the problem is to recover θ t by minimizing the cumulative error R n ( y 1:n , θ

1:n ) = T t=1 E ( y t -θ t ) 2 .

Lower-bound and previous results

In [START_REF] Roy | Online and bandit algorithms for nonstationary stochastic saddle-point optimization[END_REF], the authors first prove that using online gradient descent with fixed restart (as considered by [START_REF] Besbes | Non-stationary stochastic optimization[END_REF]) is sub-optimal in this setting. Their theorem 2 shows a cumulative error for OGD with fixed restart of order Õ(B 2 + T V (θ

1:n ) 2 + σT V (θ 1:n ) √ n)
, where B is an upper-bound on θ 1 1 . Yet, they also prove the following lower-bound.

Proposition 1 ([2, Proposition 2]) Let n ≥ 2, σ > 0, and B, C n > 0 such that min{B, C n } > 2πσ.
Then, there is a universal constant c such that, for any forecaster, there exists a sequence θ 1 , . . . , θ n such that T V (θ

1:n ) ≤ C n and R n y 1:n , θ 1:n ≥ c(B 2 +C 2 n +σ 2 log n+n 1/3 C 2/3 n σ 4/3 ).
Our aim is to address the two major challenges of ARROWS discussed previously (address general d-dimensional problems and no need to know the exact noise level σ) while achieving an optimal error of order O(n

1/3 C 2/3 n ).
An hypothetical forecaster which achieves optimal error Let m ≥ 1. We first analyse the approximation error obtained by an hypothetical forecaster that produces moving average with at most m restarts. It first computes a sequence of restart times

1 = t 1 ≤ t 2 ≤ • • • ≤ t m+1 = n + 1 such that T V θ t i :(t i+1 -1) ≤ T V (θ 1:n ) m , (3)
for all 1 ≤ i ≤ m and then forms the prediction ỹt for t ∈ {t i + 1, . . . , t i+1 } ỹt := ȳt i :(t-1) where ȳt i :(t-1) :=

1 t -t i t-1 k=t i y k . (4)
We would assume the existence of similar hypothetical forecaster for non-stationary online linear regression (section 3.1) and non-stationary online kernel regression (section 3.2) with slight variation in the prediction function. Of course this forecaster is not practical since the restart times t i are unknown.

In Theorem 5 stated in Appendix A, we show that for m ≈ n 1/3 C 2/3 n , this hypothetical forecaster achieves the same optimal error of Proposition 1,

R n (ỹ 1:n , θ 1:n ) ≤ O(n 1/3 C 2/3 n ) , (5)
as was already obtained by [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF]. Of course, it remains to estimate the restart times t i .

A meta-aggregation algorithm to learn the restart times Contrary to [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF], which uses a change point detection method, we propose to do so by using meta-aggregation algorithms from non-stationary online learning such Follow the Leading History (FLH) [START_REF] Hazan | Adaptive algorithms for online decision problems[END_REF] based on exponential weights and presented in Algorithm 1.

Algorithm 1 Follow the Leading History (FLH) [START_REF] Hazan | Adaptive algorithms for online decision problems[END_REF] Input: black box algorithm A, learning parameter η > 0 1 Init:

S 0 = ∅ 2 for t = 1, . . . ,
p t+1 (i) = p t (i) exp -η(y t -ỹt (i)) 2 .
Basically, FLH is a meta-aggregation procedure that considers a subroutine algorithm, called A, producing a prediction based on past observations. A can be any online learning algorithm that aims at minimizing the static regret, that is the excess cumulative error compared to a fixed parameter.

The role of the meta-algorithm is to learn the restarts. To do so, at each round t ≥ 1, FLH builds a new expert (step 3 of Alg. 2) that applies A on the sequence of observations y t , . . . , y n (that is by not considering the past data before round t). This new expert is assigned a weight 1/t and the weights of previous experts are normalized so that they sum to 1 (step 4). All the experts are then combined using a standard exponentially weighted average algorithm (step 6 of Alg. 2). The prediction of FLH is finally obtained (step 8) by forming a convex combination of the expert predictions. The number of active experts grow linearly with time. In Alg. 2, we also present IFLH, introduced by [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF], which improves the computational complexity by removing experts over time.

In Theorem 1, we show that a cumulative error of optimal order Õ(C

2/3

n n 1/3 ) can be achieved by applying FLH with moving averaged (4) as subroutines.

Theorem 1 Let θ 1:n ∈ R n such that T V (θ 1:n ) ≤ C n . Assume that |θ t | ≤ B for all t ≥ 1.
If moving average predictions (4) are used as subroutine of Algorithm 2, the cumulative error is upper-bounded as

R n ( y 1:n , θ 1:n ) ≤ O(n 1/3 C 2/3 n log 2 n) ,
with high probability.

Proof First, with probability 1 -δ, all |y t | = |θ t + Z t | for 1 ≤ t ≤ n are bounded by C √ log n for some constant C depending on B and σ 2 . Thus, y → (y -y t ) 2 are α-exp-concave with α = C /(log(n/δ)) for some C > 0. Let m ≈ n 1/3 C 2/3 n and t 1 , . . . , t m be as defined in (3) and [START_REF] Calandriello | Second-order kernel online convex optimization with adaptive sketching[END_REF] (see also Thm. 5). From Claim 3.1 of [START_REF] Hazan | Adaptive algorithms for online decision problems[END_REF], we have for any i = 1, . . . , m

ti+1-1 t=ti ( y t -y t ) 2 -(ỹ t (t i ) -y t ) 2 ≤ 3 log n α ≤ O(log 2 n) .
Therefore, summing over i = 1, . . . , m and using that the subroutines are moving averages (i.e., ỹt (t i ) = ȳt i :(t-1) ) and the definition of ỹt in (4), we get

n t=1 ( y t -y t ) 2 -(ỹ t -y t ) 2 ≤ O(m log 2 n) . (6)
Thus, because

Z t = y t -θ t is independent of y t R n ( y 1:n , θ 1:n ) := n t=1 E ( y t -y t + y t -θ t ) 2 = n t=1 E ( y t -y t ) 2 -(y t -θ t ) 2 = n t=1 E ( y t -y t ) 2 -(ỹ t -y t ) 2 + (ỹ t -y t ) 2 -(y t -θ t ) 2 (6) ≤ O(m(log n) 2 ) + n t=1 E (ỹ t -y t ) 2 -(y t -θ t ) 2 .
It only remains to show that the last term corresponds to R n (ỹ 1:n , θ 1:n ) := n t=1 E (ỹ t -θ t ) 2 and apply Inequality [START_REF] Calandriello | Second-order kernel online convex optimization with adaptive sketching[END_REF]. Expending the squares, it indeed yields

E (ỹ t -y t ) 2 -(y t -θ t ) 2 = E ỹ2 t + 2(θ t -ỹt )y t -θ 2 t = E ỹ2 t + 2(θ t -ỹt )(θ t + Z t ) -θ 2 t = E (ỹ t -θ t ) 2 ,
where the last equality is because E[Z t ] = 0 and Z t is independent from ỹt and θ t .

Non-Stationary Online Regression

In this section, we discuss more general problem of non-stationary online regression. We consider the following problem :

y t = g t (x t ) + Z t (7)
where g t : R d → R is a non-linear function and Z t be independent σ-subGaussian random variables in one dimension with E[Z t ] = 0. Similar to the previous section, the goal in this section would be to track the sequence of g t with g t for all t such that y t = g t (x t ) to minimize the expected cumulative error R n ( y 1:n , θ 1:n ) with respect to the unobserved output g t after n time steps which we define as follow:

R n ( y 1:n , g 1:n ) = n t=1 E ( y t -g t (x t )) 2 = n t=1 E ( g t (x t ) -g t (x t )) 2 . (8)
However, we need to remember that we observe g t (x t ) only after perturbed through some noise variable Z t . Hence, we need to decompose our regret in terms of the observed response y t . Biasvariance decomposition directly provides the decomposition in terms of the observed variable y t .

Proof is given in Appendix B.

Lemma 1 For any sequence of functions gt : R d → R for t ∈ [n] independent of Z t for all t, the cumulative error (8) can be decomposed as follows:

R n ( y 1:n , g

1:n ) = n t=1 E ( y t -y t ) 2 -(g t (x t ) -y t ) 2 + n t=1 E (g t (x t ) -g t (x t )) 2 .

Non-stationary Linear Regression

In Lemma 1, we provided the general bias-variance decomposition result for squared loss while computing expected cumulative error. In this section, we will specifically discuss the result for linear predictor θ t for all t i.e. we assume that g t is linear function for all t. Hence, the problem can be formulated as follows. At each step t ≥ 1, the learner observes x t ∈ R d , predicts y t = x t θ t and observes

y t = x t θ t + Z t (9)
Algorithm 2 IFLH -Improved Following the Leading History (binary base) [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] Input: black box algorithm A, learning parameter η > 0 Init: S 0 = ∅ for t = 1, . . . , n do

Start a new instance of algorithm A denoted A t and assign weight p t (t) = 1 t .

Define its ending time as

τ t := t + 2 k where k := min{k ≥ 0 s.t c k > 0} and t := ∞ i=1 c k 2 k is the binary representation of t.
Define the set of active experts S t := {1 ≤ i ≤ t :

τ i > t}
Normalize the weight of each active expert j ∈ S t \ {t}

p t (j) := 1 - 1 t p t (j) j∈St\{t} p t (j)
Observe x t and get the prediction ỹt (i) from each black box algorithm A i , i ∈ S t .

Predict y t = i∈St p t (i)ỹ t (i) and observe y t ∈ R.

Update the weights for each

i ∈ S t p t+1 (i) = p t (i) exp -η(y t -ỹt (i)) 2 .
where Z t be independent σ-subGaussian zero mean random variable. We assume

θ 1 , . . . , θ n ∈ R d such that T V (θ 1:n ) = n t=2 θ t -θ t-1 1 ≤ C n and θ t ≤ B for all t ≥ 1.
The goal is to control the cumulative error with respect to the unobserved outputs ỹt = x t θ t = y t -Z t . We substitute g t (x t ) with x t θ t in Equation ( 8) and denote the prediction function g t (x t ) = x t θ t = y t . Hence, the expected cumulative error R n ( y 1:n , g 1:n ) can be written as

n t=1 E ( y t -ỹt ) 2 = n t=1 E ( θ t -θ t ) x t 2 .
Hypothetical forecaster We consider an hypothetical forecaster which similar to that of 1dimensional case. It computes a sequence of restart times 3) and then forms the prediction ỹt for t

1 = t 1 ≤ t 2 ≤ • • • ≤ t m+1 = n + 1 for all 1 ≤ i ≤ m as in equation (
∈ {t i + 1, . . . , t i+1 } ỹt := x t θt , (10)
where θt = θt j :(t j+1 -1) for t j ≤ t < t j+1 and θt j :(t j+1 -1) = 1 t j+1 -t j t j+1 -1 t=t j θ t . Below in Lemma 2, we show that the cumulative error can be controlled with respect to this hypothetical forecaster.

Lemma 2 (Adaptive Restart in d-dimension) Let X, B > 0. Assume that x t ≤ X and θ t ≤ B for all t ∈ [n].
Then, there exists a sequence of restarts

1 = t 1 < • • • < t m = n + 1 such that n t=1 (x t θt -x t θ t ) 2 = m j=1 tj+1-1 t=tj ( θtj:(tj+1-1) -θ t ) x t 2 ≤ X 2 n C n m 2 + 4X 2 B 2 m ,
where θt := θt j :(t j+1 -1 ) for t j ≤ t ≤ t j+1 -1 and

θt j :(t j+1 -1) = 1 t j+1 -t j t j+1 -1 t=t j θ t .
However, this forecaster cannot be computed and is only useful for the analysis since both the restart times t i and the parameters θ t are unknown. We use meta algorithm Improved Following the Leading History (IFLH, Algorithm 2) [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] to efficiently learn the restart time which is computationally more efficient than FLH presented in Algorithm 1.

To reduce the computation complexity, there is also an associated ending time for each expert in IFLH which tells that that particular expert will no longer active after its ending time. As we only have the access to the noisy gradient, we will utilize the result presented in [22, Theorem 1] with a probabilistic upper bound on the gradient to get the final upper bound on expected cumulative loss. We provide below an upper bound on the expected cumulative error.

Theorem 2 Let n, m ≥ 1, σ > 0, B > 0, X > 0, and C n > 0. Let θ 1 , . . . , θ n such that T V (θ 1:n ) ≤ C n and θ t ≤ B. Assume that x t ≤ X for all t ≥ 1. Then, Alg. 2 [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] with Online Newton

Step [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF] as subroutine and well-tuned learning rate η > 0 satisfy

R n ( y 1:n , θ 1:n ) d 1/3 n 1/3 C 2/3 n (X 2 σ 2 B+X 2 B 2 ) 1/3 ,
with high probability.

Discussion: The result presented in Theorem 2 provides an upper bound on the expected cumulative error of Alg. 2 for non-stationary online linear regression. This generalizes the result of Baby and Wang [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] which only works for one dimensional problem. Our algorithm is adaptive to the noise parameter σ which means we do not need to know the variance σ, which is not correct for the algorithm presented in Baby and Wang [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF]. While implementing the algorithm, all we need to know is the maximum value of y t observed so far.

On Lower Bound: The lower bound presented in Baby and Wang [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] can be extended easily for general d-dimension by considering the problem of d-independent variables. This will simply add an extra multiplicative factor of d in the lower bound (Proposition 1). Our upper-bound is thus optimal in n, d, and C n . However, the dependence in σ is worse than the one of Baby and Wang [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF]. This may be due to fact that our algorithm also adapt to the noise parameter and we do not need to know σ in our algorithm. It is an interesting question to know whether our dependence in σ is optimal in our case and we leave it for future work.

Non-stationary Kernel Regression

In this section, we consider the case of nonstationary online kernel regression. For the input space X and a positive definite kernel function K : X × X → R, we denote the RKHS associated with K as H. We further denote the associated feature map φ : X → H, such that K(x, x ) = φ(x), φ(x ) H . With slight abuse of notation, we write that K(x, x ) = φ(x) φ(x ). In this section, we assume that the functions g t lie in some RKHS H corresponding to the kernel K for all t. At each step t ≥ 1, the learner observes x t ∈ R d , predicts y t = φ(x t ) θ t and observes

y t = φ(x t ) θ t + Z t , (11)
where Z t be independent σ-subGaussian zero mean random variable. The case we consider comes under well specified case as the optimal functions θ 1 , . . . , θ n ∈ H lie in the same RKHS H corresponding to the kernel K where we consider our hypothesis space. We define K nn as (K nn ) i,j = φ(x i ), φ(x j ) and λ k (K nn ) denotes the k-th largest eigenvalue of K nn . Time dependent effective dimension d ef f (λ, s, r) is defined as follows,

d ef f (λ, s, r) = T r(K s-r,s-r (K s-r,s-r + λI) -1 ) .
We also assume that T V (θ

1:n ) = n t=2 θ t - θ t-1 H ≤ C n .
The goal is to control the cumulative error with respect to the unobserved outputs ỹt = φ(x t ) θ t = y t -Z t . We substitute g t (x t ) with φ(x t ) θ t in Equation ( 8) and denote the prediction function with θ 1 , • • • , θ n such that g t (x t ) = φ(x t ) θ t = y t . Hence, the expected cumulative error R n ( y 1:n , θ 1:n ) can be written as

R n ( y 1:n , θ 1:n ) = n t=1 E ( θ t -θ t ) φ(x t ) 2 .
For our analysis, we consider a similar hypothetical forecaster as in linear regression (see Equation ( 3)). The prediction ỹt for t ∈ {t i + 1, . . . , t i+1 } is simply given as ỹt := φ(x t ) θt where θt = θt j :(t j+1 -1) for t j ≤ t < t j+1 . In the result given below in Lemma 3, we show that the expected cumulative error can be controlled with respect to this hypothetical forecaster given the adaptive restart.

Lemma 3 (Adaptive Restart in RKHS)

Let B, κ > 0. Assume that φ(x t ) 2 ≤ κ 2 , and θ t H ≤ B for all t. Then, there exists a sequence of restarts

1 = t 1 < • • • < t m = n + 1 such that n t=1 (φ(x t ) θt -φ(x t ) θ t ) 2 = m j=1 tj+1-1 t=tj ( θtj:(tj+1-1) -θ t ) φ(x t ) 2 ≤ κ 2 n C n m 2 + 4κ 2 B 2 m ,
where θt := θt j :(t j+1 -1) for t j ≤ t < t j+1 , C n ≥ n t=2 θ t -θ t-1 H , and

θt j :(t j+1 -1) = 1 t j+1 -t j t j+1 -1 t=t j θ t .
As we have discussed previously, it is not possible to compute this forecaster and it will be only useful in the analysis of the algorithm. One simply has to use a meta algorithm as in [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] to learn these restart times. However, one cannot use online newton step as the black box subroutine in this meta algorithm like it was done for linear regression as the convergence of online newton step is not known for tracking prediction functions in RKHS. Hence, we use Kernel-AWV as the black box online learner [START_REF] Gammerman | On-line prediction with kernels and the complexity approximation principle[END_REF] (see also [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF]) as subroutine in Alg. 2 to estimate the prediction function. Kernel-AWV depends on a regularization parameter λ > 0. Note that other subroutines designed for Online Kernel Regression such as Pros-N-Kons [START_REF] Calandriello | Second-order kernel online convex optimization with adaptive sketching[END_REF] or PKAWV [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF] can be used. Below, we have the following theorem regarding the adaptive regret of least square in when the predictor function lies in RKHS.

Theorem 3 Let λ > 0. For online kernel regression with square loss if for all i ∈ [n], y i ∈ [-Y, Y ] then for Algorithm 2 with Kernel-AWV [START_REF] Gammerman | On-line prediction with kernels and the complexity approximation principle[END_REF] with regularization parameter λ as subroutine, we have

s t=r f t (θ t ) - s t=r f t (u) ≤ 8Y 2 (p + 2) log n + λp θ 2 +Y 2 pd ef f (λ, s -r) log e + enκ 2 λ .
where [r, s] ⊆ [n], p ≤ log 2 (s -r + 1) + 1 and

f t (θ) = (y t -φ(x t ) θ) 2 .
With Theorem 5 and Lemma 3, we have the expression for upper bound on both the independent error terms which after combining together bound the overall expected cumulative error. Below, we provide our final bound on the expected cumulative error assuming the capacity condition, i.e., that the effective dimension satisfies

d ef f (λ, n) ≤ (n/λ) β for β ∈ (0, 1). The proof is given in Appendix C. Theorem 4 Let n, m ≥ 1, σ > 0, B > 0, κ > 0, and C n > 0. Let θ 1 , . . . , θ n such that T V (θ 1:n ) ≤ C n and θ t H ≤ B for all t ≥ 1. Assume also that φ(x t ) H ≤ κ for t ≥ 1.
Then , for well chosen η > 0, Alg. 2 with Kernel-AWV using

λ = (n/m) β β+1 satisfies R n ( y 1:n , θ 1:n ) ≤ Õ C 2(β+1) 2β+3 n n 1 2β+3 σ 2 log n δ + B 2 κ 2 + C 2 2β+3 n n 2β+1 2β+3 B 4(β+1) 2β+3 κ 2 2β+3
. with probability at least 1 -δ.

Discussion: To the best of our knowledge, this work is the first extension of non-stationary online regression to non-stationary kernel regression. After carefully looking at the bound on the expected cumulative regret term presented in Theorem 4 and comparing it with that of nonstationary online linear regression (Theorem 2), we find that as β → 0, we have λ → O(1/d) and we would have the similar dependence of C n and n in the expected cumulative error bound for linear and kernel part. However, we have a slightly worse dependent on the variance of the noise σ in the expected cumulative error bound for nonstationary online kernel regression than that of non-stationary online linear regression part. This artefact arises due to difficulty in simultaneously choosing optimal number of restart time m and regularization parameter λ. We believe that the dependence in σ in Theorem 4 can be improved further.

As discussed in [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF], the per round space and time complexities is of order O(n 2 ) for each prediction sequence corresponding to different start times. However, the method can be made computationally more efficient by the use of Nyström approximation [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF].

It is also worth pointing out that the optimal learning rate η only depends on B, κ, δ, and n and can be optimized using standard calibration techniques (e.g., doubling trick). The regularization parameter of λ on the other hand depends on the regularity of the Kernel. It can be calibrated by starting at each time steps t in Alg. 2 several new instances of Kernel-AWV, each run with a different parameter λ in a logarithmic grid.

Experiments

In this section, we evaluate our results on empirical simulations. We compare the theoretical bound with the performance of ARROWS [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] (wherever possible (1 dimension, no input)) and the procedure analyzed here, i.e., IFLH [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] with different subroutines (Online Newton Step [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF], OGD [START_REF] Zinkevich | Online convex programming and generalized infinitesimal gradient ascent[END_REF], or Azoury-Warmuth-Vovk forecaster [START_REF] Katy | Relative loss bounds for on-line density estimation with the exponential family of distributions[END_REF][START_REF] Vovk | Competitive on-line statistics[END_REF]), and online gradient descent with fixed restart [START_REF] Besbes | Non-stationary stochastic optimization[END_REF]. We test the algorithms on two different settings. The first one involves a non-stationary time series with continuous small changes in distribution which we call soft shifts. We use decaying innovation variance in order to observe how the algorithms react to a smooth change in the total variation. The second one involves hard and abrupt changes in distribution at well separated time intervals, we call the hard shifts.

Data Generation

Before presenting the experimental results and plots, we quickly here discuss the data generation process. Details of data generation process in the setting of soft shifting and hard shifting is given below.

Soft Shifts: We let θ 1 , . . . , θ n be a multivariate random walk with exponential decaying variance. We set, θ t = θ t-1 + t with t ∼ N(0, t -α I d ) multivariate normal. The total variation of this time series is

T V = n t=2 ||θ t -θ t-1 || 1 = n t=2 || t || 1 .
Hard Shifts: For generating the data used in hard shifts mechanism, we split the time series θ 1 , ..., θ n ∈ R d into M chunks such that m i is the index of the start of the i th chunk. At the start of each new chunk, all coordinates θ m i (k) for 1 ≤ k ≤ d are sampled from independent Rademacher distributions. The values of θ t are then constant within a chunk. The total variation

is T V = n t=2 ||θ t - θ t-1 || 1 = M i=2 ||θ m i -θ m i-1 || 1 .

1-Dimension (Figs 1 and 2)

We use ARROWS [START_REF] Baby | Online forecasting of total-variation-bounded sequences[END_REF] as our baseline for this part of the experiment, we compare it with our procedure proposed in Section 2, that is IFLH with moving averages as a subroutine. We recall that ARROWS was especially designed for this one dimensional setting in which it achieves the optimal rate. It also requires the variance of the noise σ 2 to be give beforehand which is not the case for our procedure. We average the predictions and the cumulative errors on 10 iterations over the time series. In all of our experiments, we consider the sub-Gaussian noise with standard deviation to be σ = 1. We have y t = θ t + Z t with Z t ∼ N(0, σ 2 ). We generate data by soft shifting and hard shifts mechanism described above.

Soft Shifts: In first part of our experiment, we generate the data by soft shifting mechanism. The parameter α, which controls how much the time-series is non-stationary, is set to be 0.3. This results in a slow decay of total variation of order n t=2 t -α = O(n 1-α ) = O(n 0.7 ) and in an upper-bound of order O(C

2/3 n n 1/3 ) = O(n 1-2 3 α ) = O(n 4/5
). We can see in Figure 2a that IFLH reacts faster to slight changes in the time series yielding a slightly smaller cumulative error than ARROWS.

Hard Shifts: For the second part of the experiment, we generate data using the hard shift mechanism. In Figures 1b and2b we test IFLH and ARROWS on a time series with equal spaced shifts, whereas in Figure 1c and2c, time intervals between shifts grow exponentially with the length of the total number of shifts. It is clear from the plots that IFLH reacts faster than ARROWS to abrupt changes and manages to adapt better to stationary portions of the time series.

Online linear regression

We test IFLH [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] on the online linear regression setting with three different subroutines: Online Gradient Descent (OGD), Online Newton Steps (ONS) as well as AWV (online ridge regression). We chose these subroutines because they are well used by the online learning community for standard stationary online linear regression. Note that ONS and AWV achieve optimal regret while this is not the case for OGD which cannot take advantage of the exp-concavity of the square loss. We compare their performances with the Online Gradient Descent with fixed restart of Besbes et al. [START_REF] Besbes | Non-stationary stochastic optimization[END_REF]. We use as batch size their theoretical result of σ -1 √ n log n/T V . We again consider two data generation mechanism as described above (soft shifting and hard shifting) to generate decision vectors θ t for all t. We take the sub-Gaussian noise to be multivariate normal with Σ = I d . We have

{m i = 100i | 1 ≤ i ≤ 100}, d = 10 (b) {m i = 2 i | 1 ≤ i ≤ 14}, d = 2 (c) {m i = 2 i | 1 ≤ i ≤ 14}, d = 10
Y t = X t θ t + Z t
with Z t ∼ N(0, Σ). We take X t to be multivariate uniformly distributed random variables X t ∼ U (-1, 1). The expected cumulative error of OGD with fixed restart grows at a rate greater than the the theoretical upper bound of O(d 1/3 n 1/3 T V 2/3 ) proved in this paper. IFLH algorithms regrets stay below the theoretical upper bound.

Soft shifts: In Figure 3, we vary the noise decaying parameter α as well as the dimension d of the time series. We can clearly notice the better performance of IFLH algorithms especially with ONS and AWV as subroutine. When α = 2 for instance, the sequence of θ t quickly converges and OGD with fixed restart continues on resetting which leads to the high divergence of its regret.

Hard shifts: In experiment 4a, we use fixed size chunks. The OGD with fixed restart algorithm performs well since the sizes of the chunks are constant and adopting a fixed restart window strategy corresponds to the setting. IFLH algorithms reacts faster to these changes and have a slightly lower regret. In 4b and 4c, we use an exponentially growing size partitions. OGD with fixed restart's regret grows at a rate bigger than the boundary line of O(d 1/3 n 1/3 T V 2/3 ). IFLH algorithms conserve a regret rate of this order.

Proof

Let ỹt be the estimate of the restarted moving average forecaster defined in Eq. ( 4) at time t. Let m ≥ 1 be the total number of batches and 1 = t 1 ≤ • • • ≤ t m+1 = n + 1 and batches be numbered as 1, • • • , m where m is the total number of batches. By Equation ( 3), the total variation of ground truth within batch i is fixed and is bounded by Cn+B m for each i, i.e. if the time interval of batch i is denoted by [t i , t i+1 -1] then by Inequality (3)

t i+1 -2 t=t i |θ t -θ t+1 | ≤ C n m .
Let us fix a batch i ∈ {1, . . . , m}. By (4), the cumulative error within the batch equals

R i := t i+1 -1 t=t i E (ỹ t -θ t ) 2 (4) = E (ȳ t i-1 :(t i -1) -θ t i ) 2 + t i+1 -1 t=t i +1 E ȳt i :(t-1) -θ t 2 = E ( θt i-1 :(t i -1) + Zt i-1 :(t i -1) -θ t i ) 2 + t i+1 -1 t=t i +1 E θt i :(t-
R i = ( θt i-1 :(t i -1) -θ t i ) 2 + E Z2 t i-1 :(t i -1) + t i+1 -1 t=t i +1 θt i :(t-1) -θ t 2 + E Z2 t i :(t-1) ≤ ( θt i-1 :(t i -1) -θ t i ) 2 + σ 2 t i -t i-1 + t i+1 -1 t=t i +1 θt i :(t-1) -θ t 2 + σ 2 t -t i
Assuming θ 0 = 0, and summing across all bins yields that the cumulative error is upper-bounded by,

R n (ỹ 1:n , θ 1:n ) := m i=1 R i ≤ m i=1 ( θt i-1 :(t i -1) -θ t i ) 2 + m i=1 t i+1 -1 t=t i +1 θt i :(t-1) -θ t 2 + m i=1 t i+1 t=t i +1 σ 2 t -t i t 1 =1 ≤ |θ 1 | 2 + m i=2 θt i-1 :(t i -1) -θ t i 2 + m i=1 t i+1 -1 t=t i +1 θt i :(t-1) -θ t 2 + m i=1 t i+1 t=t i +1 σ 2 t -t i ≤ B 2 + m i=2 θt i-1 :(t i -1) -θ t i 2 + m i=1 t i+1 -1 t=t i +1 θt i :(t-1) -θ t 2 + 2mσ 2 (2 + log n)
Then, because for all i ≥ 1 and t i+1 ≥ t ≥ t i ,

θt i :(t-1) -θ t = 1 t -t i t-1 s=t i θ s -θ t Jensen ≤ 1 t -t i t-1 s=t i θ s -θ t ≤ max s∈{t i ,...,t-1} |θ s -θ t | = max s∈{t i ,...,t-1} t-1 r=s θ r -θ r+1 ≤ max s∈{t i ,...,t-1} t-1 r=s θ r -θ r+1 = t-1 s=t i |θ s -θ s+1 | , we have R n (ỹ 1:n , θ 1:n ) ≤ B 2 + m i=2 t i -1 s=t i-1 |θ s -θ s+1 | 2 + m i=1 t i+1 -1 t=t i +1 t-1 s=t i |θ s -θ s+1 | 2 + 2mσ 2 (2 + log n) ≤ B 2 + C 2 n + m i=1 t i+1 -1 t=t i +1 t-1 s=t i |θ s -θ s+1 | 2 + 2mσ 2 (2 + log n) .
Therefore, using Inequality (3),

R n (ỹ 1:n , θ 1:n ) ≤ B 2 + C 2 n + m i=1 t i+1 -1 t=t i +1 C n m 2 + 2mσ 2 (2 + log n) ≤ B 2 + C 2 n + C 2 n m 2 m i=1 t i+1 -t i + 2mσ 2 (2 + log n) ≤ B 2 + C 2 n + nC 2 n m 2 + 2mσ 2 (2 + log n).

Now in the above equation, the choice

m = nC 2 n σ 2 (2+log n) 1/3 yields R n (ỹ 1:n , θ 1:n ) ≤ B 2 + C 2 n + 2n 1/3 C 2/3 n σ 4/3 (2 + log n) 2/3 . (12)
Remark 1 Since, we also have the boundedness assumption here on each θ i such that |θ i | ≤ B hence, it is easy to see that the bound given in the above result in Theorem 5 can be written as

R n (ỹ 1:n , θ 1:n ) ≤ B 2 + 2BC n + 2n 1/3 C 2/3 n σ 4/3 (2 + log n) 2/3 . (13)
B Non-Stationary Online Linear Regression

B.1 Bias-variance decomposition for online linear regression

Lemma 4 (Restatement of Lemma 1) For any sequence of functions gt : R d → R for t ∈ [n] independent of Z t for all t, the cumulative error R n ( y 1:n , g 1:n ) can be decomposed as follows:

R n ( y 1:n , g 1:n ) = n t=1 E ( y t -y t ) 2 -(g t (x t ) -y t ) 2 + n t=1 E (g t (x t ) -g t (x t )) 2 .
Proof Let t ≥ 1. Since y -t -g t (x t ) = Z t , which is zero mean and independent from g t (x t ) -g t (x t ), we have

E g t (x t )-y t 2 = E g t (x t )-g t (x t )+g t (x t )-y t 2 = E g t (x t )-g t (x t ) 2 +E g t (x t )-y t 2 . (14)
Therefore, by definition (8) of the cumulative error

R n ( y 1:n , g 1:n ) (8) = n t=1 E ( g t (x t ) -g t (x t )) 2 (14)
= n t=1 E ( g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 = n t=1 E ( g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 + (g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 = n t=1 E ( g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 + n t=1 E (g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 = n t=1 E ( g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 + n t=1 E gt (x t ) 2 -2g t (x t )y t -g t (x t ) 2 + 2g t (x t )y t (7) = n t=1 E ( g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 + n t=1 E[g t (x t ) 2 ] -2E[g t (x t )(g t (x t ) + Z t )] -E[g t (x t ) 2 ] + E[2g t (x t )(g t (x t ) + Z t )] = n t=1 E ( g t (x t ) -y t ) 2 -(g t (x t ) -y t ) 2 + n t=1 E (g t (x t ) -g t (x t )) 2 ,
where the last line of the proof comes from the fact that gt (x t ) is independent of Z t for all t.

B.2 Approximation error of the hypothetical forecaster

Lemma 5 (Restatement of Lemma 2) Let X, B > 0. Assume that x t ≤ X and θ t ≤ B for all t ∈ [n]. Then, there exists a sequence of restarts

1 = t 1 < • • • < t m = n + 1 such that n t=1 (x t θt -x t θ t ) 2 = m j=1 t j+1 -1 t=t j ( θt j :(t j+1 -1) -θ t ) x t 2 ≤ X 2 n C n m 2 + 4X 2 B 2 m ,
where θt := θt j :(t j+1 -1) for t j ≤ t ≤ t j+1 -1 and θt j :

(t j+1 -1) = 1 t j+1 -t j t j+1 -1 t=t j θ t .

Proof

Let m ∈ [n] be the total number of batches. Let 1 = t 1 ≤ • • • ≤ t m+1 = n + 1 be such that the total variation of the ground truth with each batch i is at most

(C n + B)/m, that is for all i ∈ [m] t i+1 -2 t=t i θ t -θ t+1 1 ≤ C n m . (15)
Therefore,

t i+1 -1 t=t i E (x t θt i :(t i+1 -1) -x t θ t ) 2 ≤ t i+1 -1 t=t i E θt i :(t i+1 -1) -θ t 2 2 x t 2 ≤ X 2 t i+1 -1 t=t i θt i :(t i+1 -1) -θ t 2 2 ≤ 4X 2 B 2 + X 2 t i+1 -2 t=t i θt i :(t i+1 -1) -θ t 2 2 ≤ 4X 2 B 2 + X 2 t i+1 -2 t=t i θt i :(t i+1 -1) -θ t 2 1 .
But, since for all i ≥ 1 and all t ∈ {t i , . . . , t i+1 -2}

θt i :(t i+1 -1) -θ t 1 Jensen ≤ 1 t i+1 -t i t i+1 -1 s=t i θ s -θ t 1 ≤ max t i ≤s≤t i+1 -1 θ s -θ t 1 ≤ t i+1 -2 t=t i θ t -θ t+1 1 (15) ≤ C n m , it yields t i+1 -1 t=t i E (x t θt i :(t i+1 -1) -x t θ t ) 2 ≤ 4X 2 B 2 + X 2 (t i+1 -t i -1) C n m 2 . ( 16
)
Summing over all batches i = 1, . . . , m concludes the proof.

B.3 Dynamic regret bound for IFLH with Online Newton Step

We present here a result from Zhang et al. [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] on the adaptive regret of Algorithm 2 that will be usefull for our regret analysis. Let us first recall their setting on non stationary online convex optimization.

Let Ω ⊂ R d be a convex compact subset of R

B.4 Proof of Theorem 2

Theorem 7 (Restatement of Theorem 2) Let n, m ≥ 1, σ > 0, B > 0, X > 0, and C n > 0. Let θ 1 , . . . , θ n such that T V (θ 1:n ) ≤ C n and θ t ≤ B. Assume that x t ≤ X for all t ≥ 1. Then, Alg. 2 [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF] with Online Newton Step [START_REF] Hazan | Logarithmic regret algorithms for online convex optimization[END_REF] as subroutine and well-tuned learning rate η > 0 satisfies

R n ( y 1:n , θ 1:n ) d 1/3 n 1/3 C 2/3 n (X 2 σ 2 B + X 2 B 2 ) 1/3 ,
with high probability.

Proof As discussed before, here the goal is to control the expected cumulative error with respect to the unobserved outputs ỹt = x t θ t = y t -Z t . Our prediction for θ t at any time instant t is denoted as θ t . Hence, the prediction for ỹt is given by y t = θ t x t and the expected cumulative error R n ( y 1:n , θ 1:n ) can be written as

R n ( y 1:n , θ 1:n ) = n t=1 E ( y t -ỹt ) 2 = n t=1 E ( θ t -θ t ) x t 2 .
Let 1 = t 1 ≤ • • • ≤ t m+1 = n + 1, θt , and θt i :(t i+1 -1) be defined as in Lemma 5. Applying Lemma 1 with gt (x t ) = θ t x t for all t, we have

R n ( y 1:n , θ 1:n ) = n t=1 E (x t θ t -y t ) 2 -(x t θt -y t ) 2 + n t=1 (x t θt -x t θ t ) 2 ≤ n t=1 E (x t θ t -y t ) 2 -(x t θt -y t ) 2 + X 2 n C n m 2 + 4X 2 B 2 m , (17)
where the second inequality is by Lemma 5. Now, we can upper-bound the first term of the right-handside by applying Theorem 6 with f t (θ) = (x t θ -y t ) 2 . Then,

∇f t (θ) = 2(x t θ -y t )x t = 2(x t θ -x t θ t -Z t )x t = 2(x t x t (θ -θ t )) -2Z t x t .
Since for all t ≥ 1, Z t are σ-subGaussian with zero-mean, we have We consider this favorable event until the end of the proof. In particular, this implies that G = 4X 2 B 2 + 2σX log n δ and that all losses f t are α-exp-concave with any parameter α ≤ 16B 2 X 2 + 2σ 2 log n δ -1 .

|Z t | ≤
Applying Theorem 6, for the choice η = α in Alg. 2, we thus get

n t=1 E (x t θ t -y t ) -(x t θt -y t ) 2 = m i=1 t i+1 -1 t=t i E (x t θ t -y t ) -(x t θt i :(t i+1 -1) -y t ) 2 Thm. 6 ≤ m (5d + 1)( log 2 n + 1) + 2 α + 5d( log 2 n + 1)GB .
Finally, substituting G and α, and plugging back into Inequality (17), we get

R n ( y 1:n , θ 1:n ) ≤ m (5d + 3)(16B 2 X 2 + 2σ 2 log n δ ) + 5d(4X 2 B + 2σX log n δ )B log 2 n +Xn C n m 2 + 4X 2 B 2 m , with probability greater than 1 -δ. Choosing m = Õ n 1/3 C 2/3 n d 1/3 (X 2 σ 2 B+X 2 B 2 ) 1/3 we get, R n ( y 1:n , θ 1:n ) ≤ Õ d 1/3 n 1/3 C 2/3 n (X 2 σ 2 B + X 2 B 2 ) 1/3 .
with high probability.

C Non-Stationary Online Kernel Regression

Below, we provide two results from Jézéquel et al. [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF] for online kernel regression with square loss. Kernel-AWV Jézéquel et al. [START_REF] Jézéquel | Efficient improper learning for online logistic regression[END_REF] computes the following estimator.

θ t = arg min θ∈H t-1 s=1 (y s -θ φ(x t )) 2 + λ θ 2 + (φ(x t ) θ) 2 , (18)
where φ : R d → H and H is RKHS corresponding to kernel K.

( y t -y t ) 2 -(θ φ(x t ) -y t ) 2 ≤ λ θ 2 H + Y 2 n k=1 log 1 + λ k (K nn ) λ
where K nn is defined as (K nn ) i,j = φ(x i ), φ(x j ) and λ k (K nn ) denotes the k-th largest eigenvalue of K nn .

Theorem 9 (Proposition 2, [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF]) For all n ≥ 1, λ > 0 and all input sequences

x 1 , • • • , x n ∈ X, n k=1 log 1 + λ k (K nn ) λ ≤ log e + enκ 2 λ d ef f (λ).
where κ = sup x∈X K(x, x) and d ef f (λ) := Tr(K nn (K nn + λI n ) -1 ).

Before proceeding to the next result, we reiterate our definition of time dependent effective dimension

d ef f (λ, s, r) = T r(K s-r,s-r (K s-r,s-r + λI) -1 ) (19)
where by abuse of notation K s-r,s-r = φ(x i ) φ(x j ) for r ≤ i ≤ s and r ≤ j ≤ s. It is also important to note that for each fixed r, d ef f (λ, s, r) is an increasing function of s -r, so that we assume that their exists an upper-bound such that for all 1 ≤ r ≤ s ≤ n,

d ef f (λ, s, r) ≤ d ef f (λ, s -r) ,
which only depends on s -r.

(θ) = (y t -φ(x t ) θ) 2 .
Proof Following the proof of Theorem 1 from [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF], we know that there exists p segments

I j = [t j , τ t j ], j ∈ [p]
with p ≤ log 2 (s -r + 1) + 1, such that t 1 = r, t j+1 = τ t j + 1, j ∈ [p -1] and τ tp ≥ s. Also, the expert (or subroutine) A t j corresponds to Kernel-AWV started at round t j and stopped at round τ t j . We denote θ

t j t j , • • • , θ t j
τt j as the sequence of solutions generated by the subroutine A t j . In other words, θ t j t

denotes the prediction at round t output by an instance of Kernel-AWV started at time t j . Following the proof of Theorem 1 of [START_REF] Zhang | Dynamic regret of strongly adaptive methods[END_REF], we have

p-1 j=1   τt j t=t j f t ( θ t ) -f t (θ t j t )   + s t=tp f t ( θ t ) -f t (θ tp t ) ≤ 1 α p j=1 log t j + 2 α s t=r+1 1 t ≤ p + 2 α log n, (20)
where α is the exp-concavity parameter of the functions f t that will be fixed later. From Theorem 8 and 9, for any j ∈ [p -1], the regret of the subroutine A t j can be upper-bounded as

τt j t=t j f t (θ t j t ) -f t (θ) ≤ λ θ 2 + Y 2 d ef f (λ, t j , τ t j ) log e + enκ 2 λ ≤ λ θ 2 + Y 2 d ef f (λ, τ t j -t j ) log e + enκ 2 λ ≤ λ θ 2 + Y 2 d ef f (λ, s -r) log e + enκ 2 λ .
Similarly for j = p, we have

s t=tp f t (θ tm t ) -f t (θ) ≤ λ θ 2 + Y 2 d ef f (λ, s -r) log e + enκ 2 λ .
Combining everything together, we have

s t=r f t ( θ t ) -f t (θ) ≤ p + 2 α log n + λp θ 2 + Y 2 pd ef f (λ, s -r) log e + enκ 2 λ .
For square loss with bounded output domain i.e. y i ∈ [-Y, Y ] for all i ∈ [n], the square loss is α-exp-concave with α = 1/8B where θt := θt j :(t j+1 -1) for t j ≤ t < t j+1 , C n ≥ n t=2 θ t -θ t-1 H , and θt j :(t j+1 -1) = 1 t j+1 -t j t j+1 -1 t=t j θ t .

Proof Let m be the total number of batches and 1 = t 1 ≤ • • • ≤ t m+1 = n + 1 such that for each batch i ∈ [m] the total variation within the batch is upper-bounded as

t j+1 -2 i=t j θ t -θ t+1 H ≤ C n m .
Following the proof of Lemma 5, we get for all i ∈ [m] where the last inequality is obtained similarly to [START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF]. Summing over the batches i = 1, . . . , m concludes the proof. Let m to be fixed later and let θt and 1 = t 1 < • • • < t m = n + 1, for t ∈ {t j , • • • , t j+1 } be as defined in Lemma 3. Applying Lemma 1 with gt (x t ) = θ t φ(x t ) for all t, followed by Lemma 6, we get R n ( y 1:n , θ From the capacity condition, we know that there exists β ∈ (0, 1) such that for all λ > 0 and n ≥ 1

d ef f (λ, n) ≤ n λ β .
Hence, using (log 2 (n) + 4) log n ≤ 8 log 2 n for n ≥ 1 and log 2 (n) + 2 ≤ 5 log n for n ≥ 2 (the error bound is true for n = 1), we get

T 1 ≤

Figure 1 :

 1 Figure 1: Examples of predictions obtained by the two considered algorithms (ARROWS in red and IFLH in blue) together with the one dimensional time-series to be predicted.

Figure 2 : 3 n

 23 Figure 2: Cumulative errors suffered in average over 10 runs by the two algorithms (ARROWS in red, IFLH in blue) together with the upper bound of order O(n 1/3 C 2/3 n ) for one dimensional time-series.

  (a) α = 1 and d = 2 (b) α = 2 and d = 2 (c) α = 2 and d = 10

Figure 3 :

 3 Figure 3: Performances of online linear regression IFLH algorithms and OGD with fixed restart on the time series generated with the Soft Shifts for various dimension d and parameter α

Figure 4 :

 4 Figure 4: Performances of online linear regression IFLH algorithms and OGD with fixed restart on the time series generated with the Hard Shifts.

Theorem 11 (,

 11 Restatement of Theorem 4) Let n, m ≥ 1, σ > 0, B > 0, κ > 0, and C n > 0. Assume thatd ef f (λ, r, s) ≤ s -r λ β , for all 1 ≤ r ≤ s ≤ n. Let θ 1 , . . . , θ n such that T V (θ 1:n ) ≤ C n and θ t H ≤ B for all t ≥ 1. Assume also that φ(x t ) H ≤ κ for t ≥ 1.Then , for well chosen η > 0, Alg. 2 with Kernel-AWV usingλ = (n/m) β β+1 and m := O C y 1:n , θ 1:n ) ≤ Õ Cwith probability at least 1 -δ.Proof Recall that the cumulative error R n ( y 1:n , θ 1:n ) can be written asR n ( y 1:n , θ 1:n ) = n t=1 E ( y t -y t ) 2 = n t=1 E ( θ t -θ t ) φ(x t ) 2 .

)T 1 + κ 2 n C n m 2 + κ 2 B 2

 1222 Last line comes from the Jensen's inequality. In the above equation, we choose λ = n m β β+1 to get the following, mY 2 log 2 n + B 2 log n m Inequality (21), it yields R n ( y 1:n , θ 1:n ) mY 2 log 2 n + B 2 log n m

  means t s=t x s and where we used that y t = θ t + Z t for any t. Using that Z t are i.i.d. random variables with E[Z t ] = 0 and E[Z 2 t ] ≤ σ 2 , we have by bias-variance decomposition

1)

+ Zt i :(t-1) -θ t 2

where the notation xt:t

  d . A sequence of convex loss functions f t : Ω → R d is sequentially optimized as follows. At each round t = 1, . . . , n, a learner chooses a parameter θ t ∈ Ω, then observes a subgradient ∇f t (θ t ) and updates θ t+1 . Learner's goal is to minimize his adaptive regret defined as the maximum static regret over intervals of length τ ≥ 1 Let n, d ≥ 1, Ω ⊆ R d , and G, B, α > 0. Let f 1 , . . . , f n : Ω → R d be a sequence of α-exp-concave loss functions such that ∇f t (θ) ≤ G for all θ ∈ Ω and 1 ≤ t ≤ n.

			s+τ -1		s+τ -1
	SA-Regret(n, τ ) := max 1≤s≤n-τ	t=s	f t (θ t ) -min θ∈Ω	t=s	f t (θ) .
	Theorem 6 (Theorem 1, [22]) Then, Algorithm 2 (i.e., Alg. 1 of Zhang et al. [22] with K = 2) with η = α and Online Newton Step as
	subroutine satisfies
	SA-Regret(T, τ ) ≤	(5d + 1)( log 2 τ + 1) + 2 α	+ 5d( log 2 τ + 1)GB log n = O d log 2 n ,
	for any τ ∈ [n].

  with probability at least 1 -δ. Hence, with probability at least 1 -δ, for all t ∈ [n] and all θ ≤ B|y t | = |θ t x t + Z t ≤ BX + 2σ log

	2σ log	n δ	,	for all t = 1, . . . , n ,
			n δ	and	∇f t (θ) ≤ 4X 2 B 2 + 2σX log	n δ	.

  Theorem 8 (Proposition 1,[START_REF] Jézéquel | Efficient online learning with kernels for adversarial large scale problems[END_REF]) Let λ, Y ≥ 0, X ⊂ R d and Y ⊂ [-Y, Y ].For any RKHS H, for n ≥ 1, for any arbitrary sequence of observations (x 1 , y 1 ), • • • , (x n , y n ) ∈ X × Y, the regret of Kernel-AWV (Equation (18),[START_REF] Jézéquel | Efficient improper learning for online logistic regression[END_REF]) is upper-bounded for all θ ∈ H as

	n
	R n (θ) :=
	t=1

  Theorem 10 (Restatement of Theorem 3) For online kernel regression with square loss if for all t ∈ [n], y t ∈ [-Y, Y ], then for the Alg. 2 with Kernel-AWV[START_REF] Jézéquel | Efficient improper learning for online logistic regression[END_REF] as subroutine with parameter λ > 0, we have for all 1 ≤ r ≤ s ≤ n and all θ ∈ H

	s t=r	f t (θ t ) -	s t=r	f t (θ) ≤ 8Y 2 (p + 2) log n + λp θ 2 + Y 2 pd ef f (λ, s -r) log e +	enκ 2 λ	,
	where p ≤ log 2 (s -r + 1) + 1 and f t

  Lemma 6 (Restatement of Lemma 3) Let B, κ > 0. Assume that φ(x t ) 2 ≤ κ 2 , and θ t H ≤ B for all t. Then, there exists a sequence of restarts1 = t 1 < • • • < t m = n + 1 such that ) θt -φ(x t ) θ t ) 2 =

	n	(φ(x t m	t j+1 -1	( θt
	t=1	j=1	t=t j

2

. Hence, substituting the value

s t=r f t ( θ t ) -f t (θ) ≤ 8Y 2 (p + 2) log n + λp θ 2 + Y 2 pd ef f (λ, s -r) log e + enκ 2 λ . j :(t j+1 -1) -θ t ) φ(x t ) 2 ≤ κ 2 n C n m 2 + 4κ 2 B 2 m ,

  E (φ(x t ) θt i :(t i+1 -1) -φ(x t ) θ t ) 2 ≤ ≤ 4κ 2 B 2 + κ 2 (t i+1 -t i -1)C n m

	t i+1 -1	t i+1 -1
		E θt i :(t i+1 -1) -θ t	2 H φ(x t ) 2 H
	t=t i	t=t i
		t i+1 -1
		≤ κ 2	E θt i :(t i+1 -1) -θ t	2 H
		t=t i

2

,

  ) θ t -y t ) 2 -(φ(x t ) θt -y t ) 2Now, we upper-bound T 1 the first term of the right-hand-side by applying Theorem 3. We only need to compute the upper-bound Y which will hold with high probability. Since for all t ≥ 1, Z t are σ-subGaussian with zero-mean, we haveE (φ(x t ) θ t -y t ) 2 -(φ(x t ) θt i :(t i+1 -1) -y t ) 2(22)≤ 8mY 2 (log 2 (n) + 4) log n + mλ(log 2 (n) + 2)B 2 + Y 2 (log 2 (n) + 2) log e + enκ 2 λ m i=1 d ef f (λ, t i+1 -t i )

	m	t i+1 -1
	=
	i=1	t=t i
					+κ 2 n	C n m	2	+ 4κ 2 B 2 m.	(21)
		:=T 1
		|Z t | ≤ 2σ log	n δ	,	for all t = 1, . . . , n ,

1:n ) = n t=1 E (φ(x t ) θ t -y t ) 2 -(φ(x t ) θt -y t ) 2 + n t=1 E (φ(x t ) θt -φ(x t ) θ t ) 2 ≤ n t=1 E (φ(x t

with probability at least 1 -δ. We consider this favorable high probability event until the end of the proof. Hence,

|y t | = |θ t φ(x t ) + Z t | ≤ Bκ + 2σ log n δ := Y for all t ∈ [n]

. Therefore, Theorem 3 entails

T 1 := n t=1 E (φ(x t ) θ t -y t ) 2 -(φ(x t ) θt -y t ) 2

  64mY 2 log 2 n + λmB 2 log n + 5Y 2 log n log e + enκ 2 λ 2 log 2 n + λmB 2 log n + 5Y 2 log n log e +

		m	t j+1 -t j
		j=1	λ
	enκ 2 λ	m 1-β n λ

β ≤ 64mY β .

Throughout the paper, the notation denotes a rough inequality which is up to universal multiplicative or additive constants and poly-logarithmic factors in n.

Acknowledgments This work was funded in part by the French government under management of Agence Nationale de la Recherche as part of the "Investissements d'avenir" program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA Institute).

Appendix

A Warmup : One Dimensional Time Series Theorem 5 (Approximation error) Let n, m ≥ 1, σ > 0, and C n > 0. Assume that 1 ≤ t 1 , . . . , t m+1 = n + 1 are defined such that (3) holds for each 1 ≤ i ≤ m. Then, for any sequence θ 1 , . . . , θ n such that T V (θ 1:n ) ≤ C n and |θ 1 | ≤ B, the hypothetical forecasts ỹt defined in Equation (4) satisfy

Therefore, optimizing m :=