André Ferrari
email: andre.ferrari@univ-cotedazur.fr

Cédric Richard
email: cedric.richard@univ-cotedazur.fr

Anthony Bourrier
email: anthony.bourrier@thalesaleniaspace.com

Ikram Bouchikhi
email: ikram.bouchikhi@univ-cotedazur.fr

Online change-point detection with kernels

Keywords: Non-parametric change-point detection, reproducing kernel Hilbert space, kernel least-mean-square algorithm, online algorithm, convergence analysis

published or not. The documents may come

Introduction

From a statistical perspective, a change-point is defined as a time instant at which some properties of a signal change, that is, the observations belong to one state up to that point, and belong to an other state after it. This change can be caused by external events, as well as by sharp transitions in the dynamics of the signal, either way it can hold critical information. Among possible applications of change point detection (CPD) we can mention medical monitoring [START_REF] Yuan | Automatic seizure detection using diffusion distance and BLDA in intracranial EEG[END_REF][START_REF] Gajic | Detection of epileptiform activity in EEG signals based on time-frequency and non-linear analysis[END_REF][START_REF] D'rozario | An automated algorithm to identify and reject artefacts for quantitative EEG analysis during sleep in patients with sleepdisordered breathing[END_REF], finance [START_REF] Bolton | Statistical fraud detection: A review[END_REF] and network security [START_REF] Tartakovsky | Detection of intrusions in information systems by sequential change-point methods[END_REF][START_REF] Yan | Catching instant messaging worms with change-point detection techniques[END_REF]. We refer interested readers to, e.g., [START_REF] Pimentel | A review of novelty detection[END_REF] or [START_REF] Truong | Selective review of offline change point detection methods[END_REF] for comprehensive reviews of CPD algorithms.

CPD algorithms can be classified, based on what is assumed to be known about the data distribution, as parametric or non-parametric. Parametric approaches assume that a model describing the data distributions of the different sists of estimating the density ratio over consecutive intervals of the time series data, was introduced [START_REF] Bouchikhi | Nonparametric online change-point detection with kernel LMS by relative density ratio estimation[END_REF]. In this algorithm, the model parameters are estimated in an online and adaptive way similar to the Kernel Least Mean Squares (KLMS) algorithm [START_REF] Liu | The kernel least-mean-square algorithm[END_REF]. The methodology showed promising and reliable detection results. In [START_REF] Bouchikhi | Kernel based online change point detection[END_REF] the authors proposed to modify the original cost function used in [START_REF] Bouchikhi | Nonparametric online change-point detection with kernel LMS by relative density ratio estimation[END_REF] in order to further improve the performance and achieve unbiasedness of the algorithm, referred to as NOUGAT (Nonparametric Online chanGepoint detection AlgoriThm).

The main contribution of this paper are as follows. After introducing the proposed algorithm (Section 2), we provide a theoretical analysis of its stochastic behavior by deriving models for the mean and the variance of the detection statistics, in the absence and the presence of a change-point (Section 3). These models are useful for several purposes: i) to assess the detection performance of the algorithm; ii) for detector design and optimization. Then we demonstrate the accuracy of these models, and we present performance comparisons with state-of-the-art algorithms on simulated and real data sets (Section 4). Finally we conclude this paper with recommendations for future research (Section 5).

The NOUGAT algorithm

In this section, we formulate the CPD problem. We review the proposed method and the online algorithm denoted as NOUGAT. Then we introduce the detection statistic. Finally we briefly discuss related works.

Problem formulation

We aim at detecting change-points in the distribution of independent random variables {y t } t∈N , y t ∈ R k , by estimating a model g(•) for r(y)-1, where r(y) = p test (y)/p ref (y) is the density ratio between the probability density p test (y) of the data on a test interval:

Y test t = (y t-(Ntest-1) , . . . , y t-1 , y t ) ∈ R k×Ntest (1)
and the probability density p ref (y) of the data on a reference interval:

Y ref t = (y t-(N ref +Ntest-1) , , . . . , y t-Ntest ) ∈ R k×N ref (2)
where N test and N ref are the number of samples in the test and reference intervals, respectively. Note that, contrary to RuLSIF [START_REF] Liu | Change-point detection in time series data by relative density-ratio estimation[END_REF], r(y) -1 is preferred to r(y) because it leads to an unbiased estimator under the no change-point hypothesis as we shall see later.

In the general case of a scalar time series {y t } t∈N , as commonly reported in the literature, we propose to proceed by considering y t = (y t , y t+1 , . . . , y t+k-1 ) ∈ R k

to take into account any dependence that may exist between successive y t .

Density-ratio estimation

The problem addressed in this paper consists of estimating a model g(•) for r(y) -1. It can be solved by fitting g(y) to r(y) -1 with respect to the squared loss:

C(g) = 1 2 E p ref (y) {(r(y) -1 -g(y)) 2 } (4)
Note that, as in [START_REF] Liu | Change-point detection in time series data by relative density-ratio estimation[END_REF], the expectation operator is defined with respect to the reference interval. By expanding (4) and then using r(y)p ref (y) = p test (y), we obtain:

C(g) = 1 2 E p ref (y) {g 2 (y)} -E ptest(y) {g(y)} + E ptest(y) {g(y)} + C (5)
where C denotes a constant value. Approximating the expected values in [START_REF] Tartakovsky | Detection of intrusions in information systems by sequential change-point methods[END_REF] by their empirical averages over the reference and test intervals data Y ref t and Y test t for any fixed t, leads to the following empirical optimization problem:

min g∈H   1 2N ref t-Ntest i=t-(N ref +Ntest-1) g 2 (y i ) - 1 N test t i=t-(Ntest-1) g(y i ) + 1 N ref t-Ntest i=t-(N ref +Ntest-1) g(y i ) + ν Ω( g H )   (6)
where H denotes an arbitrary reproducing kernel Hilbert space of real-valued functions on R. Let κ(• , •) be the reproducing kernel of H. The term ν Ω( g H ) with ν ≥ 0 is a regularization term added to promote smoothness of the solution. By virtue of the Representer Theorem [START_REF] Schölkopf | A generalized representer theorem[END_REF], any function g(•) of H that minimizes (6) can be expressed as a kernel expansion in terms of available data:

g(•, θ) = t i=t-(N ref +Ntest-1) θ i κ(• , y i ) (7)
where the θ i are parameters to be learned. This model cannot be trained efficiently in an online framework, as it needs to update both {y i } and θ as time t progresses. A standard strategy in the literature is to substitute {y i } in [START_REF] Pimentel | A review of novelty detection[END_REF] by a fixed dictionary of size L, {y ωi } L i=1 , whose elements are chosen according to some sparsification rule [START_REF] Rojo-Álvarez | Digital Signal Processing with Kernel Methods[END_REF] to represent the input data space accurately, resulting in a fixed order model,

g(•, θ) = L i=1 θ i κ ωi (•) (8)
where κ ωi (•) = κ(•, y ωi ), for all i ∈ {1, . . . , L}, are the elements of the dictionary, and

κ ω (•) = [κ ω1 (•), . . . , κ ω L (•)] .
Substituting ( 8) into ( 6), assuming a ridge parameter space regularization [START_REF] Schölkopf | A generalized representer theorem[END_REF], and minimizing (6) w.r.t. θ, we find that the optimal parameter vector θt is the solution of the following strictly convex quadratic optimization problem:

θt = arg min θ∈R L J t (θ) with J t (θ) = 1 2 θ H ref t θ + θ e • t + ν 2 θ 2 (9)
where

e • t = h ref t -h test t ( 10
)
and

h test t = 1 N test t i=t-(Ntest-1) κ ω (y i ) ( 11
)
h ref t = 1 N ref t-Ntest i=t-(Ntest+N ref -1) κ ω (y i ) (12)
H ref t = 1 N ref t-Ntest i=t-(Ntest+N ref -1) κ ω (y i )κ ω (y i ) (13)

Online density-ratio estimation

Let θ t be an estimate of the parameter vector of the density ratio model at time instant t. When t → t + 1, according to [START_REF] Basseville | Detection of Abrupt Changes -Theory and Application[END_REF], θ t+1 should be computed, as proposed in RuLSIF [START_REF] Liu | Change-point detection in time series data by relative density-ratio estimation[END_REF], by updating first (11)-( 13) and then minimizing the updated criterion J t+1 (θ). In order to reduce the computational cost, we propose as an alternative strategy to compute θ t+1 by updating θ t based on a gradient descent step of J t+1 (θ):

θ t+1 = θ t -µ∇J t+1 (θ t ) (14) = θ t -µ (H ref t+1 + νI)θ t + e • t+1 (15)
where µ > 0 is a small step size, and ∇J t+1 (θ t ) denotes the gradient of J t+1 (•) evaluated at θ t . The resulting algorithm shares similarities with the KLMS algorithm [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. The convergence behavior of the KLMS was analyzed in the case of a fixed dictionary in [START_REF] Chen | Convergence analysis of kernel LMS algorithm with pre-tuned dictionary[END_REF], and in a more general case in [START_REF] Parreira | Stochastic behavior analysis of the Gaussian kernel least-mean-square algorithm[END_REF]. Additional constraints such as sparsity have been also considered [START_REF] Gao | Kernel LMS algorithm with forward-backward splitting for dictionary learning[END_REF].

In practice, updating model g(•, θ t ) at each time instant t is a two-stage process that consists of updating both the dictionary {y ωi } L i=1 and the order L of the kernel expansion (8), followed by the update of parameter vector θ t .

Dictionary update

Numerous strategies of dictionary learning have been introduced in the online kernel filtering literature. They consist of building the dictionary {y ωi } L i=1 sequentially, by inserting selected samples y i that improve the representation of input data according to some criterion. For instance, the Approximate Linear Dependency (ALD) [START_REF] Engel | The kernel recursive least-squares algorithm[END_REF] criterion checks whether, in feature space H, the new candidate element κ(• , y t+1 ) can be well approximated by a linear combination of the elements κ(• , y ωi ) which are already in the dictionary. If not, it is added to the dictionary. The coherence rule [START_REF] Richard | Online prediction of time series data with kernels[END_REF] was introduced to avoid the computational complexity inherent to ALD. It is now considered as a state-of-the-art strategy and widely used as such. Defined by:

η = max i =j |κ(y ωi , y ωj )|,
coherence η reflects the largest correlation between the dictionary elements. The coherence rule for kernel-based dictionary selection consists of inserting y t+1 in dictionary {y ωi } L i=1 provided that its coherence remains below a threshold η 0 preset by the user: max

y ω i ∈{y ω i } L i=1 |κ(y t+1 , y ωi )| ≤ η 0 (16)
In [START_REF] Richard | Online prediction of time series data with kernels[END_REF] the authors show that the dimension of dictionaries determined with rule ( 16) is finite due to the compactness of the input space.

NOUGAT Algorithm

Depending on whether the new sample y t+1 has been inserted into the dictionary, or not, parameter vector θ t is updated similarly to [START_REF] Richard | Online prediction of time series data with kernels[END_REF]. At each time instant t, given θ t , we propose as a test statistic to consider the average of the (shifted by 1) density ratio estimators over the test interval, namely:

g t = 1 N test t i=t-(Ntest-1) g(y i , θ t ) = θ t h test t (17)
CPD is then performed by comparing g t + 1 to a given threshold ξ. The corresponding NOUGAT algorithm is described in Alg. 1.

Related works

Iteration (15) turns out to be related to the classical Geometric Moving Average algorithm (GMA) proposed in [START_REF] Roberts | Control chart tests based on geometric moving averages[END_REF]. GMA monitors a geometrically weighted estimate of the mean of y t and detects a change when the estimated mean deviates from its nominal value. Without loss of generality, the mean in the observation space can be replaced by the mean E{κ ω (y)} in the feature space defined by mapping κ ω (•), leading to:

ϑ t+1 = (1 -α)ϑ t + ακ ω (y t+1 ) (18)
However, as pointed out in [START_REF] Keriven | NEWMA: A new method for scalable model-free online change-point detection[END_REF], a drawback of GMA is that it requires to know the nominal value of E{κ ω (y)} in order to be able to calculate the associated test statistic: ϑ t -E{κ ω (y)} 2 .

Algorithm 1: NOUGAT Algorithm

1:
Step size µ, initial dictionary ω, regularization ν, thresholds η 0 and ξ 2: for t = 1, 2, . . . do if max y ω i ∈{y ω i } L i=1 |κ(y t+1 , y ωi )| > η 0 then 5:

# the dictionary remains unchanged and θ t is updated using [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF] 6: # y t+1 is added to the dictionary and θ t is updated # compute the test statistic and test 13: end if 17: end for To solve this problems in the GMA framework, a natural approach consists of comparing the estimates of E{κ ω (y t )} on two sliding windows, namely, the reference interval (2) and the test interval (1), as proposed in the Moving Average (MA) algorithm described in [START_REF] Keriven | NEWMA: A new method for scalable model-free online change-point detection[END_REF] which tracks:

θ t+1 = θ t -µ (H ref t+1 + νI)θ t + e •
9: L ← L + 1, ω L+1 = t + 1 10: θ t+1 = θ t 0 -µ (H ref t+1 + νI) θ t 0 + e •
g t+1 = θ t+1 h test
e • t 2 = h test t -h ref t 2 (19)
The approach implemented by NOUGAT differs in so far as, instead of calculating a deviation between two quantities estimated over the test and reference intervals, it estimates a unique statistic r(y) over the two intervals which is inherently equal to 1 under the null hypothesis. Note that an alternative approach, called NEWMA, proposed recently in [START_REF] Keriven | NEWMA: A new method for scalable model-free online change-point detection[END_REF], consists of testing the deviation between two GMA with different forgetting factors. The GMA with the smallest forgetting factor is used to provide an estimation of the in-control quantity. Contrarily to NOUGAT these algorithms do not explicitly take into account the covariance of the data in the feature space. From this point of view NOUGAT is similar to the Kernel Fisher Discriminant Ratio [START_REF] Harchaoui | Kernel change-point analysis[END_REF] but with a much lower computation footprint since it does not require the inversion of a covariance matrix for each t. Concerning the memory resources, NOUGAT requires to buffer N ref + N test data points. The B-statistics CPD algorithm [START_REF] Zaremba | B-tests: Low variance kernel twosample tests[END_REF] also relies on a single sliding test window, but, conversely, on multiple reference windows which considerably increases the memory requirement.

In Section 4.2 we shall compare the detection performance and computational load of the three algorithms mentioned above with the same memory footprint, namely, MA, NOUGAT and RuLSIF.

Theoretical analysis

In this section we analyze the stochastic behavior of the proposed algorithm, and derive conditions for its stability in the mean and mean square sense, in the absence and presence of a change-point. To make the analysis tractable, we shall conduct it in the case of a pre-tuned dictionary, i.e., a fixed dictionary of size L is assumed to be available beforehand. This means that L is fixed and the {y ωi } L i=1 are assumed to be deterministic. The classical Modified Independence Assumption (MIA) [START_REF] Minkoff | Comment on the "Unnecessary assumption of statistical independence between reference signal and filter weights in feedforward adaptive systems[END_REF], which assumes that H ref t+1 and θ t are statistically independent, will also be considered. Although not true in general, this assumption is commonly used to analyze adaptive constructions since it allows to simplify the derivations without constraining the conclusions. There are several results in the adaptation literature that show that performance results that are obtained under this assumption match well the actual performance of the algorithms when the step-size is sufficiently small.

Using the update rule [START_REF] Schölkopf | Estimating the support of a high-dimensional distribution[END_REF], we obtain the following recursion for θ t :

θ t+1 = I -µ(H ref t+1 + νI) θ t -µe • t+1 (20)
Define:

h test t = E ptest(y) {h test t } (21)
h ref t = E p ref (y) {h ref t } (22)
H ref t = E p ref (y) {H ref t } (23)
Taking the expected values on both sides of (20) and using the MIA we obtain the mean weight model:

m θ,t+1 = I -µ(H ref t+1 + νI) m θ,t + µ(h test t+1 -h ref t+1 ) (24)
We denote by C θ,t the correlation matrix of the weight vector θ t :

C θ,t = E{θ t θ t }
Estimating the variance of the test statistics requires a model for matrix C θ,t .

Post-multiplying (20) by its transpose, taking the expectation, and using the MIA, we obtain the following recursive expression:

C θ,t+1 = (1-µν) 2 C θ,t -µ(1-µν)(H ref t+1 C θ,t + C θ,t H ref t+1 ) + µ 2 (T +Q+Z +Z )-µ(1 -µν)(N +N ) (25)
where:

T = E{H ref t+1 θ t θ t H ref t+1 } (26)
Q = E{e • t+1 e • t+1 } (27)
Z = E{e • t+1 θ t H ref t+1 } (28)
N = E{e • t+1 θ t } (29)
In the general, all these matrices can depend on t. To simplify the notations, this dependence is dropped.

Stochastic behavior analysis under the null hypothesis

Mean analysis

Under the null hypothesis we have:

h ref t = h test t = E p ref (y) {κ ω (y)} = E ptest(y) {κ ω (y)} = h H ref t = E p ref (y) {κ ω (y) κ ω (y)} = H
and the mean weight model [START_REF] Liu | Change-point detection in time series data by relative density-ratio estimation[END_REF] simplifies to:

m θ,t+1 = I -µ(H + νI) m θ,t (30)
The mean stability of the algorithm is then ensured by using a step size µ that satisfies:

µ < 2 ζ max {H + νI} (31)
where ζ max {•} stands for the maximal eigenvalue of its matrix argument. Under this assumption m θ,t → 0. When y is Gaussian distributed, analytical expressions of h and H for a Gaussian reproducing kernel can be derived; see Appendix A.

Taking the expectation of ( 17) and assuming that θ t and h test t are independent, we get the mean of the test statistics g t :

E{g t } = h m θ,t (32)
The necessary independence assumption together with the MIA will be validated by computer simulations. Assuming (31) holds, under the null hypothesis, the asymptotic unbiasedness of the estimator implies lim t→∞ E{g t } = 0. When initializing [START_REF] Shawe-Taylor | Kernel Methods for Pattern Analysis[END_REF] with θ 0 = 0, namely, m θ,0 = 0, equation ( 30) implies m θ,t = 0 for all t. As a consequence E{g t } = 0, which means that the estimation of the density ratio r(y t ) = 1 is unbiased under the null hypothesis for all t.

Mean squared analysis

The general model of C θ,t in (25) depends on the matrices T , Q, Z and N , defined in ( 26)- [START_REF] Rojo-Álvarez | Digital Signal Processing with Kernel Methods[END_REF]. These matrices can be computed under the null hypothesis as follows.

• Denoting c θ,t = vec(C θ,t ) where vec(•) refers to the standard vectorization operator that stacks the columns of a matrix on top of each other, using the MIA and vec(ABC) = (C ⊗ A) vec(B) with ⊗ the Kronecker product, we find:

T = 1 N ref vec -1 (Γc θ,t ) + (N ref -1) HC θ,t H ( 33
)
where Γ is the (L 2 × L 2 ) matrix defined by:

Γ = E p ref (y) {κ ω (y)κ ω (y) ⊗ κ ω (y)κ ω (y)} (34)
The expression of Γ is given in Appendix B.

• Under the null hypothesis, Q is given by:

Q = N ref + N test N ref N test (H -hh ) (35)
• In the same way as T , we find that:

Z = 1 N ref vec -1 (∆m θ,t ) -hm θ,t H ( 36
)
where ∆ is the (L 2 × L) matrix defined by:

∆ = E p ref (y) {κ ω (y)κ ω (y) ⊗ κ ω (y)} (37)
The expression of ∆ is given in Appendix B.

• Application of the MIA implies N = 0.

The variance of the test statistics g t in ( 17) can be calculated using the independence assumption required previously for the computation of its mean. In particular:

var{g t } = E{g(y t ) 2 } -E{g(y t )} 2 = 1 N test tr(HC θ,t ) -(h m θ,t ) 2 (38)
The mean term m θ,t in (36), [START_REF] Schilling | Multivariate two-sample tests based on nearest neighbors[END_REF] equals zero when θ 0 = 0 and can be neglected for large values of t according to the mean analysis in Section 3.1.1. As a consequence, setting Z = 0, vectorizing (25) and using standard results on Kronecker product leads to the following proposition.

Proposition 1. Under the null hypothesis, neglecting m θ,t and assuming the MIA holds, c θ,t = vec(C θ,t ) verifies:

c θ,t+1 = Sc θ,t + µ 2 vec(Q) ( 39
)
with:

S = (1 -µν) 2 I + µ 2 N ref (Γ + (N ref -1)H ⊗ H) -µ(1 -µν)(H ⊕ H)
where

H ⊕ H = H ⊗ I + I ⊗ H.
The variance of the test statistics (17) is given by:

var{g t } = 1 N test tr(HC θ,t ) (40)
The stability of matrix S then ensures the mean-square stability of the algorithm. If the algorithm is mean-square stable, then, c θ,t converges to:

c θ,∞ = µ 2 (I -S) -1 vec(Q) (41)
The asymptotic variance of the test statistics directly derives from this result. Assuming a small step size µ we have:

S = I -2µνI -µH ⊕ H + o(µ)
Replacing in [START_REF] Chen | Sequential change-point detection based on nearest neighbors[END_REF] and using classical properties of Kronecker products, the asymptotic variance simplifies to:

var{g ∞ } = µ N test tr H vec -1 (2νI + H ⊕ H) -1 vec(Q) + o(µ) (42) = µ N test vec(H)(2νI + H ⊕ H) -1 vec(Q) + o(µ) (43)
Note that the rightmost term in [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF] 21)-( 23) depend on time t. We assume that input data y i are i.i.d with y i ∼ p 0 (•) before the change, and i.i.d with y i ∼ p 1 (•) after the change. where n 1 = t -t 0 + 1, and n 0 = N test -n 1 . In that case:

• If t < t 0 : h test t = h ref t = h 0 and H ref t = H 0 . • If t 0 ≤ t ≤ t 0 + N test -1:
h test t = 1 N test t i=t-(Ntest-1) E{κ ω (y i )} = 1 N test n 0 h 0 + n 1 h 1 time t 0 t reference test y i ∼ p 0 (•) y i ∼ p 1 (•)
y i ∼ p 0 (•) y i ∼ p 1 (•)
h ref t = h 0 and H ref t = H 0 . • If t 0 + N test ≤ t ≤ t 0 + N test + N ref -1:
the reference interval contains samples from both distributions; see figure 2. In the same way we find:

h ref t = 1 N ref (n 0 h 0 + n 1 h 1 )
H ref t = 1 N ref (n 0 H 0 + n 1 H 1 )
where:

n 1 = t -(t 0 + N test ) + 1, and n 0 = N ref -n 1 . • If t ≥ t 0 + N ref + N test , h test t = h ref t = h 1 and H ref t = H 1 .
h 0 , h 1 , H 0 and H 1 can be computed using the expressions in Appendix A when p 0 is N (µ 0 , R 0 ) and p 1 is N (µ 1 , R 1 ).

Mean analysis

A recursive model of m θ,t can be obtained by replacing h ref t , h test t , and H ref t by their expressions over time in [START_REF] Liu | Change-point detection in time series data by relative density-ratio estimation[END_REF].

• If t < t 0 : m θ,t+1 = I -µ(H 0 + νI) m θ,t (44)
• If t 0 ≤ t ≤ t 0 + N test -1: m θ,t+1 = I -µ(H 0 + νI) m θ,t + µ n 1 N test (h 1 -h 0 ) (45)
• If t 0 + N test ≤ t ≤ t 0 + N test + N ref -1: m θ,t+1 = I -µ ( n 0 N ref H 0 + n 1 N ref H 1 + νI) m θ,t + µ n 0 N ref (h 1 -h 0 ) (46) • If t ≥ t 0 + N ref + N test : m θ,t+1 = I -µ(H 1 + νI) m θ,t (47)
and the mean stability of the algorithm is ensured by using a step size µ that satisfies:

µ < 2 ζ max {H 1 + νI}
The mean of the test statistics ( 17) is, in the presence of a change-point, given by:

E{g t } =          h 0 m θ,t t < t 0 1 N test (n 1 h 1 + n 0 h 0 ) m θ,t t 0 ≤ t < t 0 + N test h 1 m θ,t t ≥ t 0 + N test (48)

Mean squared analysis

The first step consists of the computation of the matrices T , Q, Z and N in the presence of a change point.

• Following the same steps as in [START_REF] Gao | Kernel LMS algorithm with forward-backward splitting for dictionary learning[END_REF], we find:

T = 1 N ref vec -1 (Γc θ,t ) + (N ref -1) H ref t+1 C θ,t H ref t+1 ( 49
)
where Γ is defined in [START_REF] Engel | The kernel recursive least-squares algorithm[END_REF] and depends on t.

• Q can be decomposed as:

Q = Q 1 + Q 2 -(Q 3 + Q 3 ) (50)
where:

Q 1 = E{h test t+1 h test t+1 } (51)
Substituting ( 11) into (51) and expanding the expression we get:

Q 1 = 1 N test H test t+1 + (1 - 1 N test ) h test t+1 (h test t+1 )
where:

H test t = E ptest(y) {H test t }
In the same way, we find:

Q 2 = E{h ref t+1 h ref t+1 } = 1 N ref H ref t+1 + (1 - 1 N ref ) h ref t+1 (h ref t+1 )
and since the samples y i in the reference and test intervals are independent,

Q 3 = E{h test h ref } = h test t+1 (h ref t+1 )
• The matrix Z can be expanded as:

Z = E{h ref t+1 θ t H ref t+1 } -E{h test t+1 θ t H ref t+1 } (52)
The first expectation term in (52) can be computed using the MIA and the vectorization operator:

Z = 1 N ref vec -1 (∆m θ,t ) + (N ref -1) h ref t+1 m θ,t H ref t+1 -h test t+1 m θ,t H ref t+1 ( 53
)
where ∆ is defined in [START_REF] Bernstein | Matrix Mathematics: Theory, Facts, and Formulas[END_REF] and depends on t.

• Using the MIA:

N = (h test t+1 -h ref t+1 ) m θ,t
The last step consists in replacing the expressions of h ref t , h test t , and H ref t as a function of t in T , Q, Z and N . We will denote by Γ 0 (resp. Γ 1 ) and ∆ 0 (resp. ∆ 1 ) matrices Γ and ∆ in (34,37) computed for y i ∼ p 0 (•) (resp. y i ∼ p 1 (•)), see Appendix B. This leads to the following proposition. Proposition 2. Under the assumption of a change-point t 0 and assuming that the MIA holds, C θ,t is given by [START_REF] Bouchikhi | Nonparametric online change-point detection with kernel LMS by relative density ratio estimation[END_REF] where:

• If t < t 0 : T = 1 N ref vec -1 (Γ 0 c θ,t ) + (N ref -1) H 0 C θ,t H 0 Q = N ref + N test N ref N test (H 0 -h 0 h 0 ) Z = 1 N ref vec -1 (∆ 0 m θ,t ) -h 0 m θ,t H 0 N = 0 • If t ≥ t 0 + N test + N ref T = 1 N ref vec -1 (Γ 1 c θ,t ) + (N ref -1) H 1 C θ,t H 1 Q = N ref + N test N ref N test (H 1 -h 1 h 1 ) Z = 1 N ref vec -1 (∆ 1 m θ,t ) -h 1 m θ,t H 1 N = 0
The variance of the test statistics (17) is given by:

N test var{g t } =        tr(H 0 C θ,t ) -(h 0 m θ,t ) 2 t < t 0 tr(H test t C θ,t ) -(h test, t m θ,t ) 2 t 0 ≤ t < t 0 + N test tr(H 1 C θ,t ) -(h 1 m θ,t ) 2 t ≥ t 0 + N test (54)
All these expressions can be further simplified by neglecting m θ,t , specifically when t < t 0 if e.g. θ 0 = 0 and t t 0 assuming mean stability.

Simulation results

The julia code to reproduce all these experiments will be made available at github.com/andferrari.

Model validation

In this subsection, we present Monte Carlo simulations to illustrate the accuracy of the models derived in Section 3. Analytical expressions of the mean and the variance of the detection statistics under the null hypothesis are first considered. The observations y i were zero-mean two-dimensional i.i.d Gaussian vectors, with correlation coefficient equal to 0.25, and standard deviation equal to 0.5. Under these assumptions and for a Gaussian reproducing kernel, expressions of h and H are given in Appendix A. The algorithm parameters were set as follows: the bandwidth of the Gaussian kernel was σ = 0.25, the regularization parameter ν = 10 -3 , the step-size µ = 5.10 -4 . The windows lengths were set to N ref = N test = 250, and the L = 16 dictionary elements were obtained by sampling the same distribution as y i . The results were averaged over 500 Monte Carlo runs.

Figures 3 and4 compare respectively the theoretical models of the mean given by ( 30), [START_REF] Parreira | Stochastic behavior analysis of the Gaussian kernel least-mean-square algorithm[END_REF] and the variance given in Proposition 1 of the detection statistics, to Monte Carlo simulations. The asymptotic value of the variance computed from ( 41) is also reported. The initial weight vector was set to θ 0 = (0.3, 0.3) . The simulation results clearly show a good accuracy between the models and the actual performance provided by Monte Carlo simulations. These results also confirm the asymptotic unbiasedness of the estimator: E{g t } converges to 0 as expected, and validate the assumptions used in the derivations. We also provide the histogram of the detection statistics in Figure 5. Contrarily to [START_REF] Bouchikhi | Nonparametric online change-point detection with kernel LMS by relative density ratio estimation[END_REF], the histogram is very close to its Gaussian approximation as reported in this figure. Note that, as proved above, the mean converges towards zero for larger values of t. The accuracy of g t Gaussian approximation is a central result to set the threshold and guarantee a given false alarm rate using, e.g., the asymptotic expression of (38) computed using [START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF]. Figures 6 compares the asymptotic variance of the test statistic computed using [START_REF] Chen | Sequential change-point detection based on nearest neighbors[END_REF], [START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF] and its first order approximation (43) when the step size µ is small (µ = 5.10 -4 in this section). Note that this expression depends on h and H which can be computed by Monte Carlo simulations in the non Gaussian case.

For the second part of the simulations, we inserted a change-point at time instant t 0 = 25 • 10 3 by changing the input vectors correlation coefficient to 0.1 and standard deviation to 0.7. The results for the mean behavior are given in Figure 7, and for the variance in Figure 8. Both figures clearly show that the theoretical curves provided by ( 44)-( 48) and Proposition 2 match well the actual performance provided by Monte Carlo simulations, especially in the vicinity of the change point.

Performances comparison

This section aims to compare the performances of 1) dRuLSIF, a debiased version of RuLSIF obtained solving (9) at each time instant t, 2) NOUGAT, the proposed online version of dRuLSIF, and 3) MA, as defined in [START_REF] Xie | Change-point detection for high-dimensional time series with missing data[END_REF]. Note that all these algorithms share the same memory footprint. var{g } Eqs. [START_REF] Chen | Sequential change-point detection based on nearest neighbors[END_REF][START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF] var{g } first order approximation Eq. ( 43) In order to assess the performance of all these algorithms compared to a non-kernel-based algorithm, we shall now report the detection performance of a nearest-neighbors based CPD algorithm. The algorithm we selected is based on the two-sample test proposed by [START_REF] Schilling | Multivariate two-sample tests based on nearest neighbors[END_REF][START_REF] Henze | A multivariate two-sample test based on the number of nearest neighbor type coincidences[END_REF], and recently considered for CPD in [START_REF] Chen | Sequential change-point detection based on nearest neighbors[END_REF]. At each time instant t, the k-nearest-neighbors algorithm is applied to the samples in the interval (t -N ref -N test + 1, t). The test statistic is related to the number of edges N e of the graph that connect observations in the reference window with observations in the test window. Indeed, when these observations are driven by two different distributions, this graph tends to be clustered with respect to the reference and the test window, and N e is then ideally close to zero. Correspondingly, the test statistic, denoted as k-NN, is N e corrected by its mean value under equal distributions hypothesis, see e.g. [START_REF] Chen | Sequential change-point detection based on nearest neighbors[END_REF].

The observations y t were sampled from a mixture of n k-dimensional Gaussian distributions N k (m q , q -1 C q ), with q = 1, . . . , n. The weights φ q of the mixture model were sampled from a n-dimensional Dirichlet distribution of parameter α. The means m q were sampled from N k (0, I) and the covariance matrices C q from a Wishart distribution with scale matrix I and k + 2 degrees of freedom, that is, W k (I, k + 2).

The change point was set to t 0 = 400, the number of samples to n t = 700, the dimension of measurements was fixed to k = 6, the number of mixture components was set to n = 3 and α = 5. All the parameters (m q , φ q , C q ), with q = 1, . . . , n of the GMM were resampled at time t = t 0 .

For all simulations, we considered a Gaussian kernel. Its bandwidth σ was set using the median trick, that is, the median of the pairwise distances between samples governed by the same distribution as y t under the null hypothesis. A dictionary of L = 80 elements was designed by sampling the same distribution. For all Monte Carlo simulations, these parameters were kept fixed. For all algorithms, the window lengths were set to N ref = N test = 64. The regularization parameter for dRuLSIF and NOUGAT was set to ν = 10 -2 and the step size for NOUGAT was set to µ = 47 • 10 -3 . The number of nearest-neighbors of the k-NN was set to 10. This value was obtained experimentally in order to achieve the best performance.

Detection performance

Figure 9 provides the mean ± standard deviation for the four test statistics, namely, NOUGAT, dRuLSIF, MA and k-NN computed from 10 6 runs. Note that, contrarily to Figure 3, NOUGAT was initialized with θ -1 = 0 to guarantee unbiasedness under the null hypothesis as shown in Section II.A.

When comparing the ratio between the peak at t 0 + N ref (green line) and the noise level for t < t 0 (before the red line), Figure 9 reveals a slight drop in performance of NOUGAT compared to dRuLSIF. A larger loss of performance can be observed with MA and k-NN compared to NOUGAT and dRulSIF. As MA test statistic is the norm of the solution of ( 9) with H ref t = I and ν = 0, it does not take into account correlations in the feature space. In addition, MA does not take advantage of the functional approximation framework as it tests the norm of the parameters vector while NOUGAT approximates the likelihood ratio [START_REF] Yamada | Change-point detection with feature selection in high-dimensional time-series data[END_REF]. Moreover, we can observe a small detection delay between NOUGAT and the other algorithms. This can be explained by the approximate resolution of ( 9) by a gradient descent step in [START_REF] Moskvina | An algorithm based on singular spectrum analysis for change-point detection[END_REF]. This loss of performance of the online NOUGAT algorithm must be put into perspective, given its much lower computational cost compared to e.g. the offline dRuLSIF algorithm, see Section 4.2.2.

To get more insight in the performance of dRuLSIF, NOUGAT and MA, we shall now analyze the Mean Time to False Alarm (MTFA) and the Mean Time to Detection (MTD). Both are usual online performance measures [START_REF] Gustafsson | Adaptive Filtering and Change Detection[END_REF]. Let t a be the time instant of detection and t 0 the change point. They are defined as:

MTD = E{t a -t 0 | t a ≥ t 0 } (55) MTFA = E{t a | t a < t 0 } (56)
Figures 10 and11 provide the MTFA and MTD as a function of the Probability of False Alarm (PFA). The PFA was computed, for each algorithm, as the probability to detect an event at a time instant t a with t a < t 0 . The Probability of Detection (PD) was estimated as the probability to detect at least a change at a time instant t a with t 0 ≤ t a ≤ n t , i.e. the probability that the test statistics is larger than the threshold at least once. Figure 12 provides the Receiver Operating Characteristic (ROC) for the three algorithms.

Figure 10 shows that, for PFA > 0.2, the MTFA for the four algorithms is smaller than 40 samples. This means that the detection thresholds are too small and make the algorithms non-operational due to numerous false alarms.

Focusing on the case PFA < 0.2, we observe in Figure 12 that when PFA < 0.01, as expected from Figure 9, MA and k-NN reach the worst performance: for a given PFA, their PD are the smallest. We note that for 0.01 < PFA < 0.2 the PDs of NOUGAT, dRuLSIF and k-NN are almost equal to 1.

Figure 10 shows that the MTFAs are almost the same with a larger delay of 5 samples for k-NN. Figure 11 shows that NOUGAT MTD is approximately 10 samples larger than dRuLSIF and 5 samples larger than k-NN. The delay of NOUGAT, which can be observed on Figure 9 is, as explained before, due to the online update of NOUGAT. It is worthy to note that the performance of NOUGAT depends on µ and a smaller value would result in a smaller MTD.

Computational cost

This experiment aims at comparing the computational cost of the four algorithms. It is worthy to note that, as long as the averages on the reference and test windows required by dRuLSIF, NOUGAT and MA are computed recursively, their computational cost does not depends on N ref and N test . The data dimension k only intervenes when computing κ ω (•) and penalizes equally the three algorithms. On the contrary, the computational cost of k-NN strongly depends on the size of the windows, typically O((N ref + N test ) log(N ref + N test ) using a KD tree [START_REF] Bentley | Multidimensional binary search trees used for associative searching[END_REF].

Figures 13 depicts the run time for the three kernel-based algorithms as a function of the dictionary size L, when processing 1200 samples of dimension k = 6. For each value of L, the run time was calculated as the median value of 1000 runs on an Intel Core i7 3,5 GHz. The data set used in this experiment, called "Credit Card Fraud Detection", contains the 28 principal components of transactions made by European cardholders in September 2013. The data set is highly unbalanced as it contains 492 frauds out of 284,807 recorded transactions; see [START_REF] Kaggle | Credit card fraud detection[END_REF] for more details. We chose to keep only 2,000 genuine transactions, and we inserted the 492 frauds in order to create two change-points at t 0 = 1000 and t 0 = 1492 in data stream {y t } t∈N . The four most significant principal components were used as inputs (k = 4). The Gaussian kernel with kernel bandwidth σ 2 = 14, and reference and test windows of length N ref = N test = 114, were considered for all algorithms. A regularization term with ν = 10 -2 was used for NOUGAT and dRuLSIF, and the step size of NOUGAT was set to µ = 0.28.

The online dictionary update procedure described in Section 2.3 was used for all algorithms. The coherence threshold was set to η 0 = 0.7, leading to a dictionary size of L = 100. Parameter vector θ -1 was initially set to zero for NOUGAT.

The results provided in Figure 14 show that all the algorithms were able to detect the change-points. As expected, the detected change-points defined by the maximum value of each peak of the test statistics, were all in the vicinity of t 0 + N test . Nevertheless, if MA was able to detect the two change-points marked by red lines in addition to some false positive detections, it suffered from a bias that deviated its static from zero after the first change-point. dRuLSIF hardly detected the first change-point, but successfully detected the second one. NOUGAT detected both change-points with less fluctuations of its detection statistics. Finally, NOUGAT and dRuLSIF test statistics fluctuated around 0 under the null hypothesis. These results highlight the ability of the proposed algorithm to detect consecutive change-points.

Sentiment change detection in Twitter data streams

The data set used in this paragraph, called "Twitter US Airline Sentiment", is available at [START_REF] Kaggle | Twitter us airline sentiment[END_REF]. This data set contains tweets related to US Airline in February 2015, manually tagged as positive, negative and neutral. Raw tweets were first cleaned from non-ASCII characters. Stop words from Natural Language Toolkit (NLTK) corpus were also removed. Finally, tweets were represented, using a frequency-based method, in a linear space of dimension k = 50. The series {y t } t∈N was obtained by concatenating the 9178 negative-tagged tweets, the 2363 positive-tagged tweets and the 3099 neutral-tagged tweets. Parameters were set to: µ = 10 -1 , ν = 5.10 -33 , and N ref = N test = 100. A Gaussian kernel with σ 2 = 1.3 was used, along with an online dictionary learning procedure with a maximal coherence of η 0 = 10 -3 . This resulted in a dictionary of size L = 12. Parameter vector θ -1 was set to zero for NOUGAT.

Figure 15 provides the detection statistics of MA, NOUGAT and dRuLSIF. MA produced 2 false alarms and the variance of its statistics was larger than the other two methods. NOUGAT and dRuLSIF led to similar results. Note that, as expected, for the three methods, the peak at the first (negative/positive) transition was slightly higher than the peak at the second (positive/neutral) transition.

Change detection in satellite telemetry

The data set used in this experiment was provided by Thales Alenia Space. It consists of an electrical current signal produced by a panel of a geostationary satellite. The sampling period of data points is approximately 32 seconds, and the data span a time period of six months. A change point is known to occur at time instant t 0 = 177, 630. Marked by a red line, it represents a drop in the quantity of electrical current produced by the panel due to the loss of solar cells.

Figure 16 (top) partly shows the electrical current signal. The consecutive current drops observed at the beginning of the signal represent each a period of eclipse. These drops were removed using a median filter of length 600, which corresponds to the maximum duration of an eclipse. The filtered signal is shown in Figure 16 (bottom). Vectors y t of dimension k = 10 used as inputs for the detection algorithms were extracted using a sliding window as explained in [START_REF] D'rozario | An automated algorithm to identify and reject artefacts for quantitative EEG analysis during sleep in patients with sleepdisordered breathing[END_REF]. Window lengths N ref = N test = 3000 were used. This value corresponds approximately to a 1-day period, which is sufficient to capture the main stationary characteristics of the signal. These characteristics depend on changes in the distance from the panels to the Sun, and the angle of incidence of the sunlight. An online dictionary learning procedure was used with a maximal coherence value of 0.5. This resulted in a dictionary of size L = 33. The three algorithms produced false alarms. MA had a bias, dRuLSIF and NOUGAT showed similar results but with a much lower computational load for NOUGAT. Note that computational load is a key concern for this application.

Conclusion

We introduced an online kernel-based change-point detection method built upon direct estimation of the density ratio on consecutive time intervals. We analyzed its behavior in the mean and mean square sense. Finally, we evaluated its detection performance and we compared it to state-of-the-art kernel-based methodologies, MA and RuLSIF, and to another approach based on k-NN. We showed that our algorithm has a considerably lower computational complexity than dRuLSIF while ensuring comparable performance. Experiments on realworld data proved the usefulness and efficiency of our algorithm in a number of applications. These applications involved different types of data, namely, text data, raw data, and features extracted from data, showing the interest in using non-parametric techniques to perform change-point detection

We leave for future work the derivation of methods for kernel selection, and the opportunity of using a symmetric detection statistic where covariance information on the test interval would also be considered.

t+1 14 :

 14 if |g t+1 + 1| > ξthen15: flag t + 1 as a change point 16:

Figure 1 :

 1 Figure 1: An illustration of CPD procedure when the test interval contains samples driven by p 0 (•) and p 1 (•).

Figure 2 :

 2 Figure 2: An illustration of CPD procedure when the reference interval contains samples driven by p 0 (•) and p 1 (•).

Figure 3 :

 3 Figure 3: Mean of NOUGAT detection statistics obtained using model (32) and Monte Carlo simulations under the null hypothesis.

Figure 4 :

 4 Figure 4: Variance of NOUGAT detection statistics obtained using model (38) and Monte Carlo simulations under the null hypothesis.

Figure 5 :

 5 Figure 5: Histogram of NOUGAT detection statistics under the null hypothesis.

Figure 6 :

 6 Figure 6: Asymptotic variance of the test statistic compared to its first order approximation, as a function of the step size µ.

Figure 7 :

 7 Figure 7: Mean of NOUGAT detection statistics obtained using the model (48) and Monte Carlo simulations. The change is identified by the green line.

Figure 9 :

 9 Figure 9: Mean of the test statistic (± standard deviation) for dRuLSIF, NOUGAT and MA. The change point t 0 is located at the red line and t 0 + N ref at the green line.

Figure 10 :

 10 Figure 10: MTFA as a function of PFA for NOUGAT, MA and dRuLSIF.

Figure 11 :

 11 Figure 11: MTD as a function of PFA for NOUGAT, MA and dRuLSIF.

Figure 12 :Figure 13 :

 1213 Figure 12: ROC curve for NOUGAT, MA and dRuLSIF.

Figure 14 :

 14 Figure 14: Credit card fraud detection. The red lines correspond to t 0 + Ntest.

4. 3 .

 3 Experiments with real data 4.3.1. Credit card fraud detection

Figure 15 :

 15 Figure 15: Sentiment change detection in Twitter data stream. The red lines correspond to t 0 + N Test .

Figure 16 :

 16 Figure 16: Top: Telemetry. Bottom: Signal pre-processed by a median filter.

Figure 17 :

 17 Figure 17: CPD in satellite telemetry data. The red line corresponds to t 0 + N Test .

  Under the assumption of the presence of a single change-point t 0 , the analysis is conducted by comparing each time instant t to t 0 , as h ref t , h test

	t	and H ref t defined
	by (

can be efficiently computed as the solution of the Lyapunov equation (νI + H)X + X(νI + H) = Q , see

[START_REF] Bernstein | Matrix Mathematics: Theory, Facts, and Formulas[END_REF] Proposition 7.2.4

].

3.2. Stochastic behavior analysis in the presence of a change-point

  the test interval contains samples from both distributions; see figure 1. According to (21):

  Variance of NOUGAT detection statistics obtained using model (54) and Monte Carlo simulations. The change is identified by the green line.

		0.0004						Monte Carlo Model Asymptotic variance
	var{g t }	0.0002 0.0003		9.0×10 1.0×10 1.1×10 1.2×10 1.3×10 1.4×10	2.1×102.4×102.7×103.0×103.3×10
		0.0001
		0.0000	0	1×10	2×10	t	3×10	4×10	5×10
	Figure 8: 200 0.1 0.3 0.4 0.5 0.6 NOUGAT 0.2 0.0 200 0 4 5 dRuLSIF 3 2 1 200 0.5 0.6 0.7 MA 0.4 0.3 0.2 200 0 50 100 150 k-NN	300 300 300 300	400 400 400 400		500 500 500 500	600 600 600 600	700 700 700 700

Acknowledgement

This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.

Ikram Bouchikhi was partly funded by Thales Alenia Space TAS.

• If t 0 ≤ t ≤ t 0 + N test -1:

Appendix A. Computation of H and h

Considering the Gaussian reproducing kernel:

the entries of ( 23) are given by:

[H] ,q = e - and those of h by:

with , q ∈ {1, . . . , L}. These expectations can be computed for Gaussian distributed entries y i ∼ N (µ, R) using the moment generating function of a quadratic form of a Gaussian vector [START_REF] Omura | Some useful probability distributions[END_REF]:

Appendix B. Computation of Γ and ∆

The (r, s)-th entry of the (q, n)-th block of matrix Γ, that is, E{κ ωq (y i )κ ωn (y i )κ ω (y i )κ ω (y i ) } is given by:

Similarly, the r-th entry of the (q, n)-th block of ∆, that is, E{κ ωq (y i )κ ωn (y i )κ ω (y i )}, is given by: ∆ (q-1)L+r,n = e -