Sequential Learning of Principal Curves: Summarizing Data Streams on the Fly

When confronted with massive data streams, summarizing data with dimension reduction methods such as PCA raises theoretical and algorithmic pitfalls. A principal curve acts as a nonlinear generalization of PCA, and the present paper proposes a novel algorithm to automatically and sequentially learn principal curves from data streams. We show that our procedure is supported by regret bounds with optimal sublinear remainder terms. A greedy local search implementation (called slpc, for sequential learning principal curves) that incorporates both sleeping experts and multi-armed bandit ingredients is presented, along with its regret computation and performance on synthetic and real-life data.


Introduction
Numerous methods have been proposed in the statistics and machine learning literature to sum up information and represent data by condensed and simpler-to-understand quantities.Among those methods, Principal Component Analysis (PCA) aims at identifying the maximal variance axes of data.This serves as a way to represent data in a more compact fashion and hopefully reveal as well as possible their variability.PCA has been introduced by Pearson (1901) and Spearman (1904) and further developed by Hotelling (1933).This is one of the most widely used procedures in multivariate exploratory analysis targeting dimension reduction or features extraction.Nonetheless, PCA is a linear procedure and the need for more sophisticated nonlinear techniques has led to the notion of principal curve.Principal curves may be seen as a nonlinear generalization of the first principal component.The goal is to obtain a curve which passes "in the middle" of data, as illustrated by Figure 1.This notion of skeletonization of data clouds has been at the heart of numerous applications in many different domains, such as physics (Friedsam and Oren, 1989;Brunsdon, 2007), character and speech recognition (Reinhard and Niranjan, 1999;Kégl and Krzyżak, 2002), mapping and geology (Banfield and Raftery, 1992;Stanford and Raftery, 2000;Brunsdon, 2007), to name but a few.

Earlier works on principal curves
The original definition of principal curve dates back to Hastie and Stuetzle (1989).A principal curve is a smooth (C ∞ ) parameterized curve f (s) = (f 1 (s), . . ., f d (s)) in R d which does not intersect itself, has finite length inside any bounded subset of R d and is self-consistent.This last requirement means that f (s) = E[X|s f (X) = s], where X ∈ R d is a random vector and the so-called projection index s f (x) is the largest real number s minimizing the squared Euclidean distance between f (s) and x, defined by Self-consistency means that each point of f is the average (under the distribution of X) of all data points projected on f , as illustrated by Figure 2.However, an unfortunate conse-Figure 2: A principal curve and projections of data onto it.
quence of this definition is that the existence is not guaranteed in general for a particular distribution, let alone for an online sequence for which no probabilistic assumption is made.Kégl (1999) proposed a new concept of principal curves which ensures its existence for a large class of distributions.Principal curves f are defined as the curves minimizing the expected squared distance over a class F L of curves whose length is smaller than L > 0, namely, f ∈ arg inf where If E X 2 2 < ∞, f always exists but may not be unique.In practical situation where only i.i.d copies X 1 , . . ., X n of X are observed, Kégl (1999) considers classes F k,L of all polygonal lines with k segments and length not exceeding L, and chooses an estimator fk,n of f as the one within F k,L which minimizes the empirical counterpart of ∆(f ).It is proved in Kégl et al. (2000) that if X is almost surely bounded and k ∝ n 1/3 , then ∆ fk,n − ∆ (f ) = O n −1/3 .
As the task of finding a polygonal line with k segments and length at most L that minimizes ∆ n (f ) is computationally costly, Kégl et al. (2000) proposes the Polygonal Line algorithm.This iterative algorithm proceeds by fitting a polygonal line with k segments and considerably speeds up the exploration part by resorting to gradient descent.The two steps (projection and optimization) are similar to what is done by the k-means algorithm.However, the Polygonal Line algorithm is not supported by theoretical bounds and leads to variable performance depending on the distribution of the observations.
As the number k of segments plays a crucial role (a too small k leads to a poor summary of data while a too large k yields overfitting, see Figure 3), Biau and Fischer (2012) aim to fill the gap by selecting an optimal k from both theoretical and practical perspectives.Their approach relies strongly on the theory of model selection by penalization introduced  by Barron et al. (1999) and further developed by Birgé and Massart (2007).By considering countable classes {F k, } k, of polygonal lines with k segments, total length ≤ L and whose vertices are on a lattice, the optimal ( k, ˆ ) is obtained as the minimizer of the criterion crit(k, ) = ∆ n fk, + pen(k, ), is a penalty function where δ stands for the diameter of observations, w k, denotes the weight attached to class F k, and with constants c 0 , c 1 , c 2 depending on δ, the maximum length L and the dimension of observations.Biau and Fischer (2012) then prove that where Σ is a numerical constant.The expected loss of the final polygonal line fk , ˆ is close to the minimal loss achievable over F k, up to a remainder term decaying as 1/ √ n.

Motivation
The big data paradigm-where collecting, storing and analyzing massive amounts of large and complex data becomes the new standard-commands to revisit some of the classical statistical and machine learning techniques.The tremendous improvements of data acquisition infrastructures generates new continuous streams of data, rather than batch datasets.
This has drawn a large interest to sequential learning.Extending the notion of principal curves to the sequential settings opens immediate practical application possibilities.As an example, path planning for passengers' location can help taxi companies to better optimize their fleet.Online algorithms that could yield instantaneous path summarization would be adapted to the sequential nature of geolocalized data.Existing theoretical works and practical implementations of principal curves are designed for the batch setting (Kégl, 1999;Kégl et al., 2000;Kégl and Krzyżak, 2002;Sandilya and Kulkarni, 2002;Biau and Fischer, 2012) and their adaptation to the sequential setting is not a smooth process.As an example, consider the algorithm in Biau and Fischer (2012).It is assumed that vertices of principal curves are located on a lattice, and its computational complexity is of order O(nN p ) where n is the number of observations, N the number of points on the lattice and p the maximum number of vertices.When p is large, running this algorithm at each epoch yields a monumental computational cost.In general, if data is not identically distributed or even adversary, algorithms that originally worked well in the batch setting may not be ideal when cast onto the online setting (see Cesa-Bianchi and Lugosi, 2006, Chapter 4).
To the best of our knowledge, very little effort has been put so far into extending principal curves algorithms to the sequential context (to the notable exception of Laparra and Malo, 2016, in a fairly different setting and with no theoretical results).The present paper aims at filling this gap: our goal is to propose an online perspective to principal curves by automatically and sequentially learning the best principal curve summarizing a data stream.Sequential learning takes advantage of the latest collected (set of) observations and therefore suffers a much smaller computational cost.
Sequential learning operates as follows: a blackbox reveals at each time t some deterministic value x t , t = 1, 2, . . ., and a forecaster attempts to predict sequentially the next value based on past observations (and possibly other available information).The performance of the forecaster is no longer evaluated by its generalization error (as in the batch setting) but rather by a regret bound which quantifies the cumulative loss of a forecaster in the first T rounds with respect to some reference minimal loss.In sequential learning, the velocity of algorithms may be favored over statistical precision.An immediate use of aforecited techniques (Kégl et al., 2000;Sandilya and Kulkarni, 2002;Biau and Fischer, 2012) at each time round t (treating data collected until t as a batch dataset) would result in a monumental algorithmic cost.Rather, we propose a novel algorithm which adapts to the sequential nature of data, i.e., which takes advantage of previous computations.
The contributions of the present paper are twofold.We first propose a sequential principal curves algorithm, for which we derive regret bounds.We then move towards an implementation, illustrated on a toy dataset and a real-life dataset (seismic data).The sketch of our algorithm procedure is as follows.At each time round t, the number of segments of k t is chosen automatically and the number of segments k t+1 in the next round is obtained by only using information about k t and a small amount of past observations.The core of our procedure relies on computing a quantity which is linked to the mode of the so-called Gibbs quasi-posterior and is inspired by quasi-Bayesian learning.The use of quasi-Bayesian estimators is especially advocated by the PAC-Bayesian theory which originates in the machine learning community in the late 1990s, in the seminal works of Shawe-Taylor andWilliamson (1997) andMcAllester (1999a,b).The PAC-Bayesian theory has been successfully adapted to sequential learning problems, see for example Li et al. (2018) for online clustering.
The paper is organized as follows.Section 2 presents our notation and our online principal curve algorithm, for which we provide regret bounds with sublinear remainder terms in Section 3. A practical implementation is proposed in Section 4 and we illustrate its performance on synthetic and real-life data sets in Section 5. Proofs to all original results claimed in the paper are collected in Section 6.Our goal is to learn a time-dependent polygonal line which passes through the "middle" of data and gives a summary of all available observations x 1 , . . ., x t−1 (denoted by (x s ) 1:(t−1) hereafter) before time t.Our output at time t is a polygonal line ft ∈ F p depending on past information (x s ) 1:(t−1) and past predictions ( fs ) 1:(t−1) .When x t is revealed, the instantaneous loss at time t is computed as

Notation
( In what follows, we investigate regret bounds for the cumulative loss based on (2).Given a measurable space Θ (embedded with its Borel σ-algebra), we let P(Θ) denote the set of probability distributions on Θ, and for some reference measure π, we let P π (Θ) be the set of probability distributions absolutely continuous with respect to π.
For any k ∈ 1, p , let π k denote a probability distribution on F k,L .We define the prior where w 1 , . . ., w p ≥ 0 and k∈ 1,p w k = 1.
We adopt a quasi-Bayesian-flavored procedure: consider the Gibbs quasi-posterior (note that this is not a proper posterior in all generality, hence the term "quasi") as advocated by Audibert (2009) and Li et al. (2018) who then consider realisations from this quasi-posterior.In the present paper, we will rather focus on a quantity linked to the mode of this quasi-posterior.Indeed, the mode of the quasi-posterior ρt+1 is arg min , where (i) is a cumulative loss term, (ii) is a term controlling the variance of the prediction f to past predictions fs , s ≤ t, and (iii) can be regarded as a penalty function on the complexity of f if π is well chosen.This mode hence has a similar flavor to follow the best expert or follow the perturbed leader in the setting of prediction with experts (see Hutter andPoland, 2005 andCesa-Bianchi andLugosi, 2006, Chapters 3 and4) if we consider each f ∈ F p as an expert which always delivers constant advice.These remarks yield Algorithm 1.

Regret bounds for sequential learning of principal curves
We now present our main theoretical results.
The expectation of the cumulative loss of polygonal lines f1 , . . ., fT is upper-bounded by the smallest penalised cumulative loss over all k ∈ {1, . . ., p} up to a multiplicative term (1 + c 0 (e − 1)η) which can be made arbitrarily close to 1 by choosing a small enough η.However, this will lead to both a large h(f )/η in S T,h,η and a large 1 η (1 + ln f ∈Fp e −h(f ) ).In addition, another important issue is the choice of the penalty function h.For each f ∈ F p , h(f ) should be large enough to ensure a small f ∈Fp e −h(f ) while not too large to avoid overpenalization and a larger value for S T,h,η .We therefore set for each f with k segments (where |M | denotes the cardinality of a set M ) since it leads to The penalty function h(f , where c 1 , c 2 , c 3 are constants depending on R, d, δ, p (this is proven in Lemma 3, in Section 6).We therefore obtain the following corollary.
Corollary 1.Under the assumptions of Theorem 1, let and we conclude by setting .
Sadly, Corollary 1 is not of much practical use since the optimal value for η depends on inf f ∈Fp T t=1 ∆(f , x t ) which is obviously unknown, even more so at time t = 0. We therefore provide an adaptive refinement of Algorithm 1 in the following Algorithm 2.
Algorithm 2 Sequentially and adaptively learning principal curves 7: End for Theorem 2. For any sequence where c 1 , c 2 , c 3 are constants depending on R, d, δ, ln p.Let π(z) = e −z 1 {z>0} and where t ≥ 1 and c 0 = d(2R + δ) 2 .Then the procedure described in Algorithm 2 satisfies The message of this regret bound is that the expected cumulative loss of polygonal lines f1 , . . ., fT is upper-bounded by the minimal cumulative loss over all k ∈ {1, . . ., p}, up to an additive term which is sublinear in T .The actual magnitude of this remainder term is √ kT .When L is fixed, the number k of segments is a measure of complexity of the retained polygonal line.This bound therefore yields the same magnitude than (1) which is the most refined bound in the literature so far (Biau and Fischer, 2012, where the optimal values for k and L are obtained in a model selection fashion).

Implementation
The argument of the infimum in Algorithm 2 is taken over F p = ∪ p k=1 F k,L which has a cardinality of order |Q δ | p , making any greedy search largely time-consuming.We instead turn to the following strategy: given a polygonal line ft ∈ F kt,L with k t segments, we consider, with a certain proportion, the availability of ft+1 within a neighbourhood U( ft ) (see the formal definition below) of ft .This consideration is well suited for the principal curves setting since if observation x t is close to ft , one can expect that the polygonal line which well fits observations x s , s = 1, . . ., t lies in a neighbourhood of ft .In addition, if each polygonal line f is regarded as an action, we no longer assume that all actions are available at all times, and allow the set of available actions to vary at each time.This is a model known as "sleeping experts (or actions)" in prior work (Auer et al., 2003;Kleinberg et al., 2008).In this setting, defining the regret with respect to the best action in the whole set of actions in hindsight remains difficult since that action might sometimes be unavailable.Hence it is natural to define the regret with respect to the best ranking of all actions in the hindsight according to their losses or rewards, and at each round one chooses among the available actions by selecting the one which ranks the highest.Kleinberg et al. (2008) introduced this notion of regret and studied both the full-information (best action) and partial-information (multi-armed bandit) settings with stochastic and adversarial rewards and adversarial action availability.They pointed out that the EXP4 algorithm (Auer et al., 2003) attains the optimal regret in adversarial rewards case but has a runtime exponential in the number of all actions.Kanade et al. (2009) considered full and partial information with stochastic action availability and proposed an algorithm that runs in polynomial time.
In what follows, we materialize our implementation by resorting to "sleeping experts" i.e., a special set of available actions that adapts to the setting of principal curves.
Let σ denote an ordering of |F p | actions, and A t a subset of the available actions at round t.We let σ(A t ) denote the highest ranked action in A t .In addition, for any action f ∈ F p we define the reward r f ,t of f at round t, t ≥ 0 by It is clear that r f ,t ∈ (0, c 0 ).The convention from losses to gains is done in order to facilitate the subsequent performance analysis.The reward of an ordering σ is the cumulative reward of the selected action at each time and the reward of the best ordering is max σ T t=0 r σ(At),t (respectively, E max σ T t=1 r σ(At),t when A t is stochastic).
Our procedure starts with a partition step which aims at identifying the "relevant" neighbourhood of an observation x ∈ R d with respect to a given polygonal line, and then proceeds with the definition of the neighbourhood of an action f .We then provide the full implementation and prove a regret bound.
Partition For any polygonal line f with k segments, we denote by V = (v 1 , . . ., v k+1 ) its vertices and by s i , i = 1, . . ., k the line segments connecting v i and v i+1 .In the sequel, we use f ( V) to represent the polygonal line formed by connecting consecutive vertices in V if no confusion arises.Let V i , i = 1, . . ., k + 1 and S i , i = 1, . . ., k be the Voronoi partitions of R d with respect to f , i.e., regions consisting of all points closer to vertex v i or segment s i .Figure 5 shows an example of Voronoi partition with respect to f with 3 segments.
Neighbourhood For any x ∈ R d , we define the neighbourhood N(x) with respect to f as the union of all Voronoi partitions whose closure intersects with two vertices connecting the projection f (s f (x)) of x to f .For example, for the point x in Figure 5, its neighbourhood N(x) is the union of S 2 , V 3 , S 3 and V 4 .In addition, let N t (x) = {x s ∈ N (x) , s = 1, . . ., t.} be the set of observations x 1:t belonging to N (x) and Nt (x) be its average.Let D(M ) = sup x,y∈M ||x − y|| 2 denote the diameter of set M ⊂ R d .We finally define the local grid We can finally proceed to the definition of the neighbourhood U( ft ) of ft .Assume ft , v jt:kt+1 ), where vertices of (ii) belong to Q δ,t (x t ) while those of (i) and (iii) do not.The neighbourhood U( ft ) consists of f sharing vertices (i), (iii) with ft , but can be equipped with different vertices (ii) in Q δ,t (x t ), i.e., same number of segments, j t − i t + 1 increase segments by 1 unit.
In Algorithm 3, we initiate the principal curve f1 as the first component line segment whose vertices are the two farthest projections of data x 1:t 0 (t 0 can be set to 2 or 3 in practice) on the first component line.The reward of f at round t in this setting is therefore r f ,t = c 0 − ∆(f , x t 0 +t ).Algorithm 3 has an exploration phase (when I t = 1) and an exploitation phase (I t = 0).In the exploration phase, it is allowed to observe rewards of all actions and to choose an optimal perturbed action from the set F p of all actions.In the exploitation phase, only rewards of a part of actions can be accessed and rewards of others are estimated by a constant, and we update our action from the neighbourhood U ft−1 Algorithm 3 A locally greedy algorithm to sequentially learn principal curves 1: Input parameters: p > 0, R > 0, L > 0, > 0, α > 0, 1 > β > 0 and any penalty function h 2: Initialization: Given (x t ) 1:t 0 , obtain f1 as the first principal component 3: For t = 2, . . ., T

5:
Let i.e., sorting all f ∈ F p in descending order according to their perturbed cumulative reward till t − 1.

6:
If I t = 1, set A t = F p and ft = σt (A t ) and observe r ft,t 7: rf,t = r f ,t for f ∈ F p .

8:
If I t = 0, set A t = U( ft−1 ), ft = σt (A t ) and observe r ft,t 9: where H t denotes all the randomness before time t and cond (t) = f ∈ F p : P ft = f |H t > β .In particular, when t = 1, we set rf,1 = r f ,1 for all f ∈ F p , U f0 = ∅ and rσ 1 (U( f0)) ,1 ≡ 0. |F p | when p is large.In addition, this local search will be enough to account for the case when x t locates in U ft−1 .The parameter β needs to be carefully calibrated since it should not be too large to ensure that the condition cond(t) is non-empty, otherwise all rewards are estimated by the same constant and thus lead to the same descending ordering of tuples for both t−1 s=1 rf,s , f ∈ F p and t s=1 rf,s , f ∈ F p .Therefore, we may face the risk of having ft+1 in the neighbourhood of ft even if we are in the exploration phase at time t + 1. Conversely, very small β could result in large bias for the estimation of r f ,t .Note that the exploitation phase is close yet different to the label efficient prediction (Cesa-Bianchi et al., 2005, Remark 1.1) since we allow an action at time t to be different from the previous one.Neu and Bartók (2013) have proposed the Geometric Resampling method to estimate the conditional probability P ft = f |H t since this quantity often does not have an explicit form.However, due to the simple exponential distribution of z f chosen in our case, an explicit form of P ft = f |H t is straightforward.
Then the procedure described in Algorithm 3 satisfies the regret bound 4 ).
The proof of Theorem 3 is presented in Section 6.The regret is upper bounded by a term of order is the price to pay for the local search (with a proportion 1 − ) of polygonal line ft in the neighbourhood of the previous ft−1 .If = 1, we would have that ĉ0 = c 0 and the last two terms in the first inequality of Theorem 3 would vanish, hence the upper bound reduces to Theorem 2. In addition, our algorithm achieves an order that is smaller (from the perspective of both the number |F p | of all actions and the total rounds T ) than Kanade et al. (2009) since at each time, the availability of actions for our algorithm can be either the whole action set or a neighbourhood of the previous action while Kanade et al. (2009) consider at each time only partial and independent stochastic available set of actions generated from a predefined distribution.

Numerical experiments
We illustrate the performance of Algorithm 3 on synthetic and real-life data.Our implementation (hereafter denoted by slpc -Sequential Learning of Principal Curves) is conducted with the R language and thus our most natural competitor is the R package princurve (which is the algorithm from Hastie and Stuetzle, 1989).We let p = 20, The spacing δ of the lattice is ajusted with respect to data scale.

Synthetic data in high dimension
We also apply our algorithm on a data set x t ∈ R 6 , t = 1, 2, . . ., 200 in higher dimension.It is generated uniformly along a parametric curve whose coordinates are where t takes 100 equidistant values in [0, 2π].To the best of our knowledge, Hastie and Stuetzle (1989), Kégl (1999) and Biau and Fischer (2012) only tested their algorithm on 2-dimensional data.This example aims at illustrating that our algorithm also works on higher dimensional data.Table 2 shows the regret for the ground truth, princurve and slpc.In addition, Figure 6 shows the behaviour of slpc (green) on each dimension.ground truth princurve slpc 3.290 (0) 14.204 (0) 6.797 (0.409) Table 2: Regret (cumulative loss) on synthetic data in higher dimension (average over 10 trials, with standard deviation in brackets).princurve is deterministic.Seismic data Seismic data spanning long periods of time are essential for a thorough understanding of earthquakes.The "Centennial Earthquake Catalog" (Engdahl and Villaseñor, 2002) aims at providing a realistic picture of the seismicity distribution on Earth.It consists in a global catalog of locations and magnitudes of instrumentally recorded earthquakes from 1900 to 2008.We focus on a particularly representative seismic active zone (a lithospheric border close to Australia) whose longitude is between E130 • to E180 • and latitude between S70 • to N30 • , with T = 218 seismic recordings.As shown in Figure 7, slpc recovers nicely the tectonic plate boundary.Lastly, since no ground truth is available, we use the R 2 coefficient to assess the performance (residuals are replaced by the squared distance between data points and their projections onto the principal curve).The average over 10 trials is 0.990.
Back to the synthetic data setting Figure 8 presents the predicted principal curve ft+1 for both princurve (red) and slpc (green).The output of princurve yields a curve which does not pass in "the middle of data" but rather bends towards the curvature of the data cloud: slpc does not suffer from this behavior.To better illustrate the way slpc works between two epochs, Figure 9 focuses on the impact of collecting a new data point on the principal curve.We see that only a local vertex is impacted, whereas the rest of the principal curve remains unaltered.This cutdown in algorithmic complexity is one the key assets of slpc.magnitude of earthquakes, etc.) used in the present paper may be downloaded from this website.

Back to seismic data
Daily commute data The identification of segments of personal daily commuting trajectories can help taxi or bus companies to optimise their fleets and increase frequencies on segments with high commuting activity.Sequential principal curves appear to be an ideal tool to address this learning problem: we test our algorithm on trajectory data from the University of Illinois at Chicago2 .The data is obtained from the GPS reading systems carried by two of the lab members during their daily commute for 6 months in the Cook county and the Dupage county of Illinois.Figure 11 presents the learning curves yielded by princurve and slpc on geolocalization data for the first person, on May 30 in the data set.A particularly remarkable asset of slpc is that abrupt curvature in the data sequence is perfectly captured, whereas princurve does not enjoy the same flexibility.Again, we use the R 2 coefficient to assess the performance (where residuals are replaced by the squared distance between data points and their projections onto the principal curve).The average over 10 trials is 0.998.

Proofs
This section contains the proof of Theorem 2 (note that Theorem 1 is a straightforward consequence, with η t = η, t = 0, . . ., T ) and the proof of Theorem 3 (which involves intermediary lemmas).Let us first define for each t = 0, . . ., T the following forecaster sequence Note that f t is an "illegal" forecaster since it peeks into the future.In addition, denote by the polygonal line in F p which minimizes the cumulative loss in the first T rounds plus a penalty term.f is deterministic while f t is a random quantity (since it depends on z f , f ∈ F p drawn from π).If several f attain the infimum, we choose f T as the one having the smallest complexity.We now enunciate the first (out of three) intermediary technical result.
where 1/η −1 = 0 by convention.The second and third inequality is due to respectively the definition of f T and f T .Hence where the second inequality is due to E[Z f T ] = 0 and 1 ηt − 1 η t−1 > 0 for t = 0, 1, . . ., T since η t is decreasing in t in Theorem 2. In addition, for y ≥ 0, one has Hence, for any y ≥ 0 P sup where u = f ∈Fp e −h(f ) .Therefore, we have We thus obtain Next, we control the regret of Algorithm 2.
Lemma 2. Assume that z f is sampled from the symmetric exponential distribution in R, (1 + η t−1 c 0 (e − 1)) E ∆ f t , x t . (7) Proof.Let us denote by the instantaneous loss suffered by the polygonal line ft when x t is obtained.We have where the inequality is due to the fact that ∆(f , x) ≤ d(2R + δ) 2 holds uniformly for any f ∈ F p and x ∈ B(0, √ dR).Finally, summing on t on both sides and using the elementary inequality e x ≤ 1 + (e − 1)x if x ∈ (0, 1) concludes the proof.
Lemma 3.For k ∈ 1, p , we control the cardinality of set f ∈ F p , K(f ) = k as where V d denotes the volume of the unit ball in R d .
Proof.First, let N k,δ denote the set of polygonal lines with k segments and whose vertices are in where the second inequality is a consequence to the elementary inequality p k ≤ pe k k combined with Lemma 2 in Kégl (1999).
We now have all the ingredients to prove Theorem 1 and Theorem 2.
First, combining ( 6) and ( 7) yields that Assume that η t = η, t = 0, . . ., T and h(f and the second inequality is obtained with Lemma 1.By setting Since E ∆( f t , x t ) ≤ c 0 for any t = 1, . . ., T , we have where cond(t) c denotes the complement of set cond(t).The first inequality above is due to the assumption that for all f ∈ A t ∩ cond(t), we have P σt (A t ) = f H t ≥ β.For t = 1, the above inequality is trivial since rσ 1 (U( f0)) ,1 ≡ 0 by its definition.Hence, for t ≥ 1, one has E r ft,t H t = E r σt (Fp),t H t , I t = 1 + (1 − )E r σt (At),t H t , I t = 0 Summing on both sides of inequality (8) over t terminates the proof of Lemma 4. Proof.By the definition of rf,t in Algorithm 3, for any f ∈ F p and t ≥ 1, we have where in the second inequality we use that r f ,t ≤ c 0 for all f and t, and that P ft = f H t ≥ β when f ∈ U ft−1 ∩ cond(t).The rest of the proof is similar to those of Lemma 1 and Lemma 2. In fact, if we define by ∆ (f , x t ) = ĉ0 − rf,t , then one can easily observe the following relation when I t = 1 (similar relation in the case that I t = 0) ft = σt (F p ) = arg max = arg min Then applying Lemma 1 and Lemma 2 on this newly defined sequence ∆ ft , x t , t = 1, . . .T leads to the result of Lemma 5.
The proof of the upcoming Lemma 6 requires the following submartingale inequality: let Y 0 , . . .Y T be a sequence of random variable adapted to random events H 0 , . . ., H T such that for 1 ≤ t ≤ T , the following three conditions hold Then for any λ > 0, .
The proof can be found in Chung and Lu (2006, Theorem 7.3).
Lemma 6. Assume that 0 < β <  Finally, noticing that max f ∈Fp T t=1 (r f ,t − rf,t ) ≤ c 0 T almost surely, we terminate the proof of Lemma 6.
Proof of Theorem 3. Assume that p > 6, T ≥ 2|F p | 2 and let With those values, the assumptions of Lemma 4, Lemma 5 and Lemma 6 are satisfied.
(a) A too small k.(b) Right k.(c) A too large k.

Figure 3 :
Figure 3: Principal curves with different number k of segments.

A
parameterized curve in R d is a continuous function f : I −→ R d where I = [a, b] is a closed interval of the real line.The length of f is given by a sequence of data, where B(c, R) stands for the 2 -ball centered in c ∈ R d with radius R > 0. Let Q δ be a grid over B(0, √ dR), i.e., Q δ = B(0, √ dR) ∩ Γ δ where Γ δ is a lattice in R d with spacing δ > 0. Let L > 0 and define for each k ∈ 1, p the collection F k,L of polygonal lines f with k segments whose vertices are in Q δ and such that L(f ) ≤ L. Denote by F p = ∪ p k=1 F k,L all polygonal lines with a number of segments ≤ p, whose vertices are in Q δ and whose length is at most L. Finally, let K(f ) denote the number of segments of f ∈ F p .This strategy is illustrated by Figure 4.

Figure 5 :
Figure 5: An example of a Voronoi partition.

Figure 6 :
Figure 6: slpc (green line) on synthetic data in higher dimension from different perspectives.Black dots represent recordings x 1:99 , the red dot is the new recording x 200 .
Figure 10 is taken from the USGS website 1 and gives the global earthquakes locations on the period 1900-1999.The seismic data (latitude, longitude,

Figure 7 :
Figure 7: Seismic data.Black dots represent seismic recordings x 1:t , red dot is the new recording x t+1 .
Figure 8: Synthetic data -Black dots represent data x 1:t , red point is the new observation x t+1 .princurve (solid red) and slpc (solid green).
Figure 9: Synthetic data -Zooming in: how a new data point impacts only locally the principal curve.