Hervé Cardot

Peggy C Énac

Pierre-André Zitt
email: pierre-andre.zitt@u-bourgogne.fr

Peggy Cardot

Cenac

Efficient and fast estimation of the geometric median in Hilbert spaces with an averaged stochastic gradient algorithm

Keywords: CLT, functional data, geometric quantiles, high dimension, L 1 -median, online algorithms, recursive estimation, Robbins-Monro algorithm, spatial median

With the progress of measurement apparatus and the development of automatic sensors it is not unusual anymore to get large samples of observations taking values in high dimension spaces such as functional spaces. In such large samples of high dimensional data, outlying curves may not be uncommon and even a few individuals may corrupt simple statistical indicators such as the mean trajectory. We focus here on the estimation of the geometric median which is a direct generalization of the real median in metric spaces and has nice robustness properties. The geometric median being defined as the minimizer of a simple convex functional that is differentiable everywhere when the distribution has no atom, it is possible to estimate it with online gradient algorithms. Such algorithms are very fast and can deal with large samples. Furthermore they also can be simply updated when the data arrive sequentially. We state the almost sure consistency and the L 2 rates of convergence of the stochastic gradient estimator as well as the asymptotic normality of its averaged version. We get that the asymptotic distribution of the averaged version of the algorithm is the same as the classic estimators which are based on the minimization of the empirical loss function. The performances of our averaged sequential estimator, both in terms of computation speed and accuracy of the estimations, are evaluated with a small simulation study. Our approach is also illustrated on a sample of more than 5000 individual television audiences measured every second over a period of 24 hours.

Introduction

With the progress of measurement apparatus, the development of automatic sensors and the increasing storage performances of computers it is not unusual anymore to get large samples of functional observations. For example Cardot et al. (2010a) analyze a sample of more than 18000 electricity consumption curves measured every half hour over a period of two weeks. Our study is motivated by the estimation of the central point of a sample of n = 5423 vectors of R d , with d = 86400, which correspond to individual television audiences measured every second over a period of 24 hours.

In such large samples of high dimensional data, outlying curves may not be uncommon and a even few individuals may corrupt simple statistical indicators such as the mean trajectory or the principal components [START_REF] Gervini | Robust functional estimation using the median and spherical principal components[END_REF]). Detecting these atypical curves automatically is not straightforward in such a high dimensional and large sample context and considering directly robust techniques is an interesting alternative. There are many robust location indicators in the multivariate setting [START_REF] Small | A survey of multidimensional medians[END_REF]) but most of them require high computational efforts to be estimated, even for small sample sizes, when the dimension is relatively large. For example, [START_REF] Fraiman | Trimmed means for functional data[END_REF] have extended the notion of trimmed means to a functional context in order to get robust estimators of the mean profile.

In order to deal with the dimensionality issue and to reduce the computation time, [START_REF] Cuevas | Robust estimation and classification for functional data via projection-based depth notions[END_REF] have proposed random projection techniques in the context of maximal depth estimators and studied their properties via simulation studies. Note that sub-sampling approaches based on survey sampling with unequal probability sampling designs have also been proposed in the literature in order to reduce the computational time [START_REF] Chaouch | Design-based estimation for geometric quantiles with application to outliers detection[END_REF]).

We focus here on the geometric median, also called L 1 -median or spatial median, which is a direct generalization of the real median proposed by [START_REF] Haldane | Note on the median of a multivariate distribution[END_REF] and whose properties have been studied in details by [START_REF] Kemperman | The median of a finite measure on a Banach space[END_REF]. It can be defined even if the random variable does not have a finite first order moment and it has nice robustness properties since its breakdown point is equal to 0.5. As noted in [START_REF] Small | A survey of multidimensional medians[END_REF], one drawback of the geometric median is that it is not affine equivariant. Nevertheless, it is invariant to translation and scale changes and thus is well adapted to functional data which are observed with the same units at each instant of time. In a functional context, consistent estimators of the L 1 -median have been proposed by [START_REF] Kemperman | The median of a finite measure on a Banach space[END_REF], [START_REF] Cadre | Convergent estimators for the L 1 -median of a Banach valued random variable[END_REF] and [START_REF] Gervini | Robust functional estimation using the median and spherical principal components[END_REF]. Iterative estimation algorithms have been developed by [START_REF] Gower | Algorithm as 78: The mediancentre[END_REF], [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF] in the multivariate setting and by [START_REF] Gervini | Robust functional estimation using the median and spherical principal components[END_REF] for functional data. This latter algorithm requires to invert at each step matrices whose dimension is equal to the dimension d of the data and thus requires important computational efforts. The algorithm proposed by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF] is much faster and only requires O(nd) operations at each iteration, where n is the sample size. Nevertheless, these estimation procedures are not adapted when the data arrive sequentially, they need to store all the data and they cannot be simply updated.

In this paper, we explore another direction. The geometric median being defined as the minimizer of a simple functional that is differentiable everywhere when the distribution has no atom, it is possible to estimate it with online gradient algorithms. Such algorithms are very fast and can be simply updated when the data arrive sequentially. There is a vast literature on stochastic gradient algorithms which mainly focus on finite dimensional situations (see [START_REF] Kushner | Stochastic Approximation Methods for Constrained and Unconstrained Systems[END_REF], [START_REF] Ruppert | A Newton-Raphson version of the multivariate Robbins-Monro procedure[END_REF], [START_REF] Benveniste | Adaptive Algorithms and Stochastic Approximations[END_REF][START_REF] Ljung | Stochastic Approximation and Optimization of Random Systems[END_REF], [START_REF] Duflo | of Applications of Mathematics[END_REF], [START_REF] Kushner | Stochastic approximation and recursive algorithms and applications[END_REF], [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF] in the multivariate case, and [START_REF] Arnaudon | Stochastic algorithms for computing means of probability measures[END_REF] on manifolds). The literature is much less abundant when one has to consider online observations taking values in a functional space (usually an infinite dimensional Banach or Hilbert space) and most works focus on linear algorithms [START_REF] Walk | An invariance principle for the Robbins-Monro process in a Hilbert space[END_REF], [START_REF] Dippon | The averaged Robbins-Monro method for linear problems in a Banach space[END_REF], [START_REF] Smale | Online learning algorithms[END_REF]).

It is also known in the multivariate setting that averaging procedures can lead to efficient estimation procedure under additional assumptions on the noise and when the target is defined as the minimizer of a strictly convex function [START_REF] Polyak | Acceleration of stochastic approximation[END_REF], [START_REF] Pelletier | Asymptotic almost sure efficiency of averaged stochastic algorithms[END_REF]). There is little work on averaging when considering random variables taking values in Hilbert spaces and, as far as we know, they only deal with linear algorithms (Dippon and [START_REF] Dippon | The averaged Robbins-Monro method for linear problems in a Banach space[END_REF]). Nevertheless, it has been noted in an empirical study whose aim was to estimate the geometric median with functional data (Cardot et al. (2010a)) that averaging could improve in an important way the accuracy of the estimators.

The paper is organized as follows. We first fix notations, give some properties of the geometric median and present our stochastic gradient algorithm as well as its averaged version. We also note that our study extends directly to the estimation of geometric quantiles defined by [START_REF] Chaudhuri | On a geometric notion of quantiles for multivariate data[END_REF]. In a third section we state the almost sure consistency and the L 2 rates of convergence of the stochastic gradient estimators as well as the asymptotic normality of its averaged version. We get that the asymptotic distribution of the averaged version of the algorithm is the same as the classic estimators. A fourth section is devoted to a small simulation study which aims at comparing the performances of our estimator with the static algorithm developed by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF]. The comparison is performed according to two points of view, for the same sample size and for the same computation time. We also analyze a real example with a large sample of individual television audiences measured every second over a period of 24 hours. The proofs are gathered in Section 6.

The algorithms and some properties of the geometric median

Definitions and assumptions

Let H be a separable Hilbert space such as R d or L 2 (I), for some closed interval I ⊂ R. We denote by ., . its inner product and by • the associated norm.

The geometric median m of a random variable X taking values in H is defined by (see [START_REF] Kemperman | The median of a finite measure on a Banach space[END_REF]):

m := arg min u∈H E [ X -u -X ] . (1)
Note that this general definition (1) does not assume the existence of the first order moment of X . We suppose from now on that the following assumptions are fulfilled.

A1.

The random variable X is not concentrated on a straight line: for all v ∈ H, there is w ∈ H such that v, w = 0 and Var( w, X ) > 0.

A2. The law of X is a mixing of two "nice" distributions :

µ X = λµ c + (1 -λ)µ d ,
where

-µ c is not strongly concentrated around single points: if B(0, A) is the ball {α ∈ H, α ≤ A}, and Y is a random variable with law µ c , ∀A, ∃C A ∈ [0, ∞), ∀α ∈ B(0, A), E Y -α -1 ≤ C A . -µ d is a discrete measure, µ d = ∑ i p i δ α i .
We denote by D the support of µ d and assume that m / ∈ D.

As shown in [START_REF] Kemperman | The median of a finite measure on a Banach space[END_REF], assumption (A1) ensures that the median m is uniquely defined. The second assumption could probably be relaxed, but it is general enough for most natural examples. As noted in [START_REF] Chaudhuri | Multivariate location estimation using extension of R-estimates through U-statistics type approach[END_REF], the conditions on µ c are satisfied when H = R d , with d ≥ 2, whenever µ c has a bounded density on every compact subset of R d . More precisely, this property is closely related to small ball probabilities since

E Y -α -1 = ∞ 0 P Y -α ≤ t -1 dt. If P [ Y -α ≤ ǫ] ≤ Cǫ d
, for some small ǫ and some positive constant C, it is easy to check that

E Y -α -β < ∞, whenever 0 ≤ β < d.
When H = L 2 (I), the dimension is not finite and small ball probabilities have been derived for some particular classes of Gaussian processes (see [START_REF] Nazarov | Exact l 2 -small ball asymptotics of gaussian processes and the spectrum of boundary-value problems[END_REF] for a recent reference). In this case, by symmetry of the distribution, the median m is equal to the mean, and many processes satisfy, for positive constants C 1 , C 2 , C 3 , C 4 which depend on the process under study,

P [ Y -m ≤ ǫ] ≤ C 1 ǫ C 4 exp(-C 2 ǫ -C 3 ), (2)
so that E Ym -β < ∞, for all positive β. Similar properties of shifted small balls, for α close to m, can be found in [START_REF] Li | Gaussian processes: Inequalities, small ball probabilities and applications[END_REF].

Some convexity and robustness properties of the median

In this section we derive quantitative convexity bounds which will be useful in the proofs.

As a consequence, we are also able to bound for the gross sensitivity error, which is a classical robustness indicator (see [START_REF] Huber | Robust Statistics[END_REF]).

Recalling the definition of the median (eq. ( 1)), let us denote by G : H → R the function we would like to minimize:

G(α) := E [ X -α -X ] . ( 3
)
This function is convex since it is a convex combination of convex functions. However, convexity is not sufficient to get the convergence of the algorithm. Under assumptions (A1) and (A2) this function can be decomposed in two parts:

G(α) = λG c (α) + (1 -λ)G d (α),
where the discrete part G d (α) = ∑ i p i ( x iαx i ) has been isolated. The first part is Fréchet differentiable everywhere [START_REF] Kemperman | The median of a finite measure on a Banach space[END_REF]), so G is differentiable except on D, the support of the discrete part µ d . We denote by

Φ = λΦ c + (1 -λ)Φ d its Fréchet derivative, Φ(α) := ∇ α G = -E X -α X -α .
Remark 1. It will be useful to define Φ on the set D. If x ∈ D, we define G x by "forgetting" x,

G x (y) = ∑ i,x i =x p i ( x i -y -x i ).
This function is Fréchet differentiable in x, and we let

Φ d (x) = ∑ i,x i =x p i x -x i x -x i .
In the vocabulary of convex analysis, we have just chosen a particular subgradient of G on the set D of non-differentiability. It is easily seen that:

∀x, y, G(y) -G(x) ≥ Φ(x), y -x , (4)
which asserts that Φ is a subgradient. A short proof of inequality (4) is given in Section 6.1.

The median m is then the unique solution of the nonlinear equation,

Φ(α) = 0. ( 5
)
To exhibit some useful strong convexity and robustness properties of the median we need to introduce the Hessian of functional G, for α ∈ H \ D. It is denoted by Γ α , maps H to H and it is easy to check (see [START_REF] Koltchinskii | M-estimation, convexity and quantiles[END_REF] for the multivariate case and Gervini (2008) for the functional one) that

Γ α = E 1 X -α I H - (X -α) ⊗ (X -α) X -α 2 ,
where I H is the identity operator in H and u ⊗ v(h) = u, h v, for u, v and h belonging to H. The operator Γ α is not compact but it is bounded when E Xα -1 < ∞.

If we define h = h/ h , and P h the projection onto the orthogonal complement of h,

h, Γ α h = h 2 E 1 α -X 1 - h, α -X 2 α -X 2 = h 2 E 1 α -X Ph(α -X) 2 α -X 2 . ( 6
)
We can now state a strong convexity property of functional G which can be seen as an extension to an infinite dimensional setting of Proposition 4.1 in [START_REF] Koltchinskii | M-estimation, convexity and quantiles[END_REF].

Proposition 2.1. Recall that B(0, A) is the ball of radius A in H. Under assumptions A1 and A2, there is a strictly positive constant c A , such that:

∀α ∈ B(0, A) \ D, ∀h ∈ H, c A h 2 ≤ h, Γ α h ≤ C A h 2 .
In other words, G is strictly convex in H and it is strongly convex on any bounded set, as shown in the following corollary.

Corollary 2.2. Assume hypotheses of Proposition 2.1 are fulfilled. For any strictly positive A, there is a strictly positive constant c A such that:

∀α 1 , α 2 ∈ B(0, A) 2 , Φ(α 2 ) -Φ(α 1 ), α 2 -α 1 ≥ c A α 2 -α 1 2 .
As a particular case of Proposition 2.1, we get that there exist two strictly positive con-

stants 0 < c m ≤ C m ≤ E X -m -1 < ∞, such that c m h 2 ≤ h, Γ m h ≤ C m h 2 . ( 7
)
As noted in [START_REF] Kemperman | The median of a finite measure on a Banach space[END_REF], the geometric median has a 50 % breakdown point.

Furthermore, an immediate consequence of ( 7) is that operator Γ m has a bounded inverse.

Thus, the gross error sensitivity, which is also a classical indicator of robustness, is bounded for the median in a separable Hilbert space. Indeed, thanks to the expression derived in [START_REF] Gervini | Robust functional estimation using the median and spherical principal components[END_REF], it is bounded as follows,

sup z∈H Γ -1 m z -m z -m ≤ 1 c m .

The algorithms

Given X 1 , X 2 , . . . , X n , n independent copies of X, a natural estimator of m is the solution m n of the empirical version of (5),

n ∑ i=1 X i -m n X i -m n = 0.
The solution m n is defined implicitly and is found by iterative algorithms.

We propose now an alternative and simple estimation algorithm which can be seen as a stochastic gradient algorithm [START_REF] Ruppert | A Newton-Raphson version of the multivariate Robbins-Monro procedure[END_REF]; [START_REF] Duflo | of Applications of Mathematics[END_REF]) and is defined as follows

Z n+1 = Z n + γ n X n+1 -Z n X n+1 -Z n = Z n -γ n U n+1 , (8)
with a starting point that can be random and bounded, e.g. Z 0 = X 0 1 { X 0 ≤M} for some positive constant M fixed in advance, or deterministic. If X n+1 = Z n , we set U n+1 = 0 and Z n+1 = Z n so the algorithm does not move. The sequence of descent steps γ n controls the convergence of the algorithm. The direction U n+1 is an "estimate" of the gradient Φ of G at Z n since the conditional expectation given the sequence of σ-algebra

F n = σ(Z 1 , . . . , Z n ) = σ(X 1 , . . . , X n ) satisfies E [U n+1 |F n ] = Φ(Z n ). (9)
Note that our particular choice of subgradient Φ ensures that this equality always holds.

Defining now by ξ the sequence of "errors" in these estimates,

ξ n+1 = Φ(Z n ) -U n+1 , (10)
algorithm ( 8) can also be seen as a non linear Robbins-Monro algorithm,

Z n+1 = Z n + γ n (-Φ(Z n ) + ξ n+1 ) . (11)
Thanks to ( 9) and ( 10), the sequence (ξ n ) is a sequence of martingale differences. Let us note that the bracket of the associated martingale satisfies,

E ξ n+1 2 F n = E U n+1 2 F n + Φ(Z n ) 2 -2 Φ(Z n ), E [U n+1 |F n ] = E U n+1 2 F n -Φ(Z n ) 2 ≤ 1 -Φ(Z n ) 2 ≤ 1. ( 12
)
Our second algorithm consists in averaging all the estimated past values,

Z n+1 = Z n + 1 n + 1 Z n+1 -Z n with Z 0 = 0, so that Z n = 1 n ∑ n i=1 Z i .
Remark 2. An extension of the notion of quantiles in Euclidean and Hilbert spaces has been proposed by [START_REF] Chaudhuri | On a geometric notion of quantiles for multivariate data[END_REF]. In such spaces, quantiles are associated to a direction and a magnitude specified by a vector v ∈ H, such that v < 1. The geometric quantile of X, say m v , corresponding to direction v and magnitude v is defined, uniquely under previous assumptions, by

m v = arg min u∈H E [ X -u + X -u, v ] .
If v = 0 one recovers the geometric median. When v is close to one, m v is a (directed) extreme quantile. In any case, m v is characterized by:

Φ v (m v ) = Φ(m v ) -v = 0,
so that it can be naturally estimated with the following stochastic algorithm

m v n+1 = m v n + γ n X n+1 -m v n X n+1 -m v n + v ,
as well as with its averaged version.

3 Convergence results

Almost sure convergence of the stochastic gradient algorithm

We first state the almost sure consistency of our sequence of estimators Z n under classical and general assumptions on the descent steps γ n .

Theorem 3.1. If ( A1) and ( A2) hold, and if (γ n ) n∈N satisfies the usual conditions:

∑ n≥1 γ n = ∞, ∑ n≥1 γ 2 n < ∞, then lim n→∞ Z n -m = 0, a.s.

Rates of convergence and asymptotic normality

We present now the rates of convergence of the stochastic gradient algorithm as well as the asymptotic distribution of its averaged version. The proofs are given in Section 6. More specific sequences (γ n ) are considered and we suppose from now on that γ n = c γ n -α , where c γ is a positive constant and α ∈ ( 1 2 , 1). We need one additional assumption to get these rates of convergence:

A3. There is a positive constant A such that ∃C A ∈ [0, ∞), ∀h ∈ B(0, A), E X -(m + h) -2 ≤ C A . ( 13
)
This assumption is not restrictive when the dimension is strictly larger than two as discussed in Section 2.1.

The following proposition states that, on events of arbitrarily high probability, the functional estimator Z n attains the classical rates of convergence in quadratic mean (see [START_REF] Duflo | of Applications of Mathematics[END_REF], theorem 2.2.12) for the multivariate case) up to a logarithmic factor.

Proposition 3.2. Assume (A1), ( A2) and (A3). Then, there exist an increasing sequence of events (Ω N ) N∈N , and constants C N , such that Ω = N∈N Ω N , and

∀N, E 1 Ω N Z n -m 2 ≤ C N γ n ln n ∑ k=1 γ k ≤ C N ln(n) n α . Remark 3. An immediate consequence of Proposition 3.2 is that Z n -m 2 = O P ln(n) n α .
Assumption (A3) is needed to bound the difference between G and its quadratic approximation, in a neighborhood of m as stated in the following Lemma.

Lemma 3.3. Assume (A3) is in force. Then, ∀h ∈ B(0, A), Φ(m + h) = Γ m (h) + O h 2 .
Finally, Theorem 3.4 stated below probably gives the most important result of this work.

It is shown that the averaged estimator Z n and the classic static estimator m n have the same asymptotic distribution. Consequently, for large sample sizes, it is possible to get, very quickly, estimators which are as efficient, at first order, as the slower static one m n . Note that the asymptotic distribution of m n has been derived in the multivariate case by Haberman (1989), Theorem 6.1. For variables taking values in a Hilbert space, such asymptotic distribution has only been proved for a particular case, when the support of X is a finite dimensional space (Theorem 6 in Gervini ( 2008)).

Theorem 3.4. Assume (A1), ( A2) and (A3). Then,

√ n Z n -m L ---→ n→∞ N 0, Γ -1 m ΣΓ -1 m , with, Σ = E (X -m) X -m ⊗ (X -m) X -m .
Note that with (7), operator Γ -1 m is well defined, it is bounded and positive.

Illustrations on simulated and real data 4.1 A simulation study

A simple simulation study is performed to check the good behavior of the averaged stochastic estimator and to make a comparison with the static estimator developed by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF]. Two points of view are considered. The first classic one consists in evaluating the performances of these two different approaches for different sample sizes. The second one, which is the point of view that should be adopted when computation time matters, consists in comparing the accuracy of both approaches when the allocated computation time is fixed in advance. We use (R Development Core Team ( 2010)) and the function spatial.median from the library ICSNP to estimate the median with the algorithm developed by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF].

For simplicity, we consider random variables taking values in R 3 and make simulations of Gaussian random vectors with median m = (0, 0, 0) and covariance matrix:

Γ =    3 2 1 2 4 -0.5 1 -0.5 2    .
In order to compare the accuracy of the different algorithms, we compute the following estimation error,

R( m) = m -m , ( 14
)
where m is an estimator of m.

Our averaged estimator depends on the tuning parameters α and c γ which control the

descent steps γ k = c γ k -α .
It is well known that for the particular case α = 1, the choice of parameter c γ is crucial for the convergence and depends on the second derivative of G in m which is unknown in practice. As usually done for such procedures, we fix α = 3/4 and focus on the choice of c γ . We also run in parallel the algorithm for 10 initial points chosen randomly in the sample and then select the best estimate m which corresponds to the minimum value of

α → 1 n n ∑ i=1 ( X i -α -X i ) ,
which is the empirical version of (3).

Fixed sample sizes

We perform 1000 simulations for different sample sizes, n = 250, n = 500 and n = 2000. Table 1 presents basic statistics for the estimation errors (first quartile Q 1 , median and third quartile Q 3 ), according to criterion ( 14), for the algorithm by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF] and our averaged procedure considering different values for c γ ∈ {0.2, 0.6, 1, 2, 5, 10, 15, 25, 50, 75}.

At first, we can note that even for moderate sample sizes the averaged procedure performs well in comparison with the Vardi and Zhang estimator which only gives slightly better estimations. We can also remark that the averaged stochastic estimator is not much sensitive to the tuning parameter c γ which can take values in the interval [2, 75] without modifying the performances of the estimator. As a matter of fact, we noted on simulations that interesting values for c γ are around or above E [ Xm ] , which is about 2.7 for this particular simulation study.

Fixed computation time

Even if both algorithms require computation times which are O(nd) (for n observations in dimension d), the averaged stochastic gradient approach is much faster (on the same computer, with procedures coded in the same language). For example, in previous simulations, if the sample size is n = 1000, the averaged estimator is about 30 times faster.

When the dimension gets larger the difference is even more impressive, as we will see in the next section.

Let us suppose the allocated time for computation is limited and fixed in advance, say 1 second, and compare the sample sizes that can be handled by the different algorithms. The static estimator by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF] can deal with n = 150 observations, whereas our recursive algorithm, coded in the language, can take into account n = 4500 observations so that it will gives much better estimates of the median, as seen in Table 1. Finally, if the algorithm is coded in C and called from , then it is at least 20 times faster than its analogue, so that it can deal with at least n = 90000 observations, during the same second.

Estimation of the median television audience profile

The analysis of audience profiles for different channels, or different days of the year, is an essential tool to understand the consumers' habits as regards television. The French company Médiamétrie provides official television audience rates in France. Médiamétrie works with a panel of about 9000 individuals and the television sets of these individuals are equipped with sensors that measure the audience of the different channels at a second scale.

A sample of around 7000 people is drawn every day in this panel and the television consumption of the people belonging to this sample is recorded every second. The data are then sent sequentially to Médiamétrie during the night. Survey sampling techniques with unequal probability sampling designs are used by Médiamétrie to select the sample and thus the i.i.d assumption is clearly not satisfied. Nevertheless, our aim is just to give an illustration of the ability of our averaged stochastic algorithm to deal with a large sample of very high dimensional data. Moreover, Médiamétrie has noted in these samples the presence of some atypical behaviors so that robust techniques may be helpful.

We focus our study on the estimation of the television audience profile during the 6th september 2010. After removing from the sample people that did not watch television at all on that day, we finally get a sample of size n = 5423. For each element i of the sample, we have a vector X i ∈ {0, 1} 86400 , where 86400 is the number of seconds within a day, and zero values correspond to seconds during the day where i is not watching television.

A classical audience indicator is given by the mean profile, drawn in Figure 1, which is simply the proportion of people watching television at every second over the considered period of time. We compare this classical indicator with the geometric median, whose estimation is drawn in black in Figure 1. We can first note that both estimators have the same shape along time, showing three peaks of audience during the day with higher audience rates between 8 and 10 PM. Estimated values are smaller for the geometric median which is less sensitive to small perturbations and outliers. This also indicates that the distribution of the individual audience curves is not symmetric around the mean profile.

From a computational point of view, even if the database is huge, it takes less than one minute for our algorithm to converge whereas we were not able to perform the estimation with the static estimator developed by [START_REF] Vardi | The multivariate L 1 -median and associated data depth[END_REF] because of memory contraints. The value of the tuning parameter was chosen to be c γ = 400, it leads to a value of about 92 for the empirical loss criterion.

Concluding remarks

The experimental results confirm that averaged recursive estimators of the geometric median relying on stochastic gradient approaches are of particular interest when one has to deal with large samples of data and potential outliers. Furthermore, when the allocated computation time is limited and fixed in advance and the data arrive online these techniques can deal, in a recursive way, with larger sample sizes and finally provide estimations that are much more accurate than static estimation procedures. We have also noted in the simulation experiment that they are not very sensitive to the value of the tuning parameter c γ .

One could imagine many directions for future research that certainly deserve further attention. Taking advantage of the rapidity of our estimation procedure, one could use resampling techniques, similar to the bootstrap, in order to approximate the asymptotic distribution of the estimator given in Theorem 3.4 and then build pointwise confidence intervals. Proving rigorously the validity of such techniques is far beyond the scope of this paper.

Our procedure can also be extended readily for online clustering, adapting the well known MacQueen algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] be optimized is not convex anymore, it can be proved that stochastic gradient approaches converge almost surely to the set of stationary points [START_REF] Cardot | Fast clustering of large datasets with sequential k-medians : a stochastic gradient approach[END_REF]) and thus are interesting candidates for online clustering.

Another direction of interest is online estimation of the conditional geometric median when real covariates are available. For instance, the age or the size of the city where individual live are known by Médiamétrie and it can be possible to take such information into account in order to get varying time regression models that can also be estimated in a very fast way thanks to sequential approaches.

Finally, as noted in Section 4.2 the independence condition is rarely satisfied for real studies. A direction that deserves further investigation is to determine under which dependence conditions our results, such as Theorems 3.1 and 3.4, remain true.

Proofs

Convexity -Proofs of Proposition 2.1 and Corollary 2.2

We first show that Φ is a subgradient of G. For points x / ∈ D, this is clear since G is Fréchet differentiable.

Pick a point x 0 in D and recall that

Φ d (x 0 ) = ∑ i,x i =x 0 p i x 0 -x i x 0 -x i
. We have,

Φ d (x 0 ), y -x 0 = ∑ i,x i =x 0 p i x 0 -x i , y -x 0 x 0 -x i = ∑ i,x i =x 0 p i x 0 -x i , y -x i x 0 -x i -∑ i,x i =x 0 p i x 0 -x i ≤ ∑ i,x i =x 0 p i y -x i -∑ i,x i =x 0 p i x 0 -x i ≤ G d (y) -G d (x 0 ), so that Φ d is a subgradient of G d .
The upper bound in Proposition 2.1 follows immediately from ( 6) and the assumption (A2).

For the lower bound, thanks to (6), we only need to prove:

∀α ∈ B(0, A), ∀u, u = 1, u, Γ α u = E P u (X -α) 2 X -α 3 ≥ c A , ( 15
)
where P u is the projection on the orthogonal of u. This quantity is small when Xα is in span(u).

Recall that (by (A1)), X is not supported on a line. Consider the set of subspaces K ⊂ H satisfying: ∀x ∈ K, Var( x, X ) = 0. Suppose that this set is non-empty, and let H ′ be a maximal element in it (this exists by Zorn's lemma). The orthogonal of H ′ has at least dimension 2 (otherwise, we get a contradiction to A1). Let v 1 , v 2 be two orthogonal vectors

in H ′⊥ . Let v t = cos(t)v 1 + sin(t)v 2 . The map t → Var( v t , X )
is continuous on a compact set. Its minimum cannot be zero (since this would contradict the maximality of H ′ ). Therefore there exists a c such that, for all unit v in the plane spanned by

(v 1 , v 2 ), Var( X, v ) ≥ c.
The orthogonal of u (an hyperplane) and the (2-dimensional) plane spanned by v 1 and v 2 necessarily intersect: there exists a unit vector v ∈ span(v 1 , v 2 ) such that u, v = 0. Therefore, for all y ∈ H, P u y 2 ≥ y, v 2 . In particular,

P u (X -α) 2 ≥ X -α, v 2 .
Suppose first that X is a.s. bounded by K. Then

E v, X -α 2 X -α 3 ≥ 1 (A + K) 3 E v, X -α 2 .
It is easily seen that the last term is bounded below by Var( v, X ) ≥ c and (15) holds with

c A = 1 (K + A) 3 c.
To get rid of the boundedness assumption on X, we can just choose K large enough so that

Var( v, X1 X ≤K ) is strictly positive for v = v 1 , v 2 .
The corollary is a consequence of Proposition 2.1 and of the fact that Φ is a subgradient. Indeed, the inequality holds for Φ c by interpolation: for an elementary proof, define

α t = (1 -t)α 1 + tα 2 , and write Φ(α 2 ) -Φ(α 1 ) = 1 0 f ′ (t)dt where f (t) = Φ(α t ). One can then apply (15), t by t, with α = α t and u = α 2 -α 1 α 2 -α 1 . Moreover, Φ d (α 2 ) -Φ d (α 1 ), α 2 -α 1 = Φ d (α 2 ), α 2 -α 1 + Φ d (α 1 ), α 1 -α 2 ≥ G(α 2 ) -G(α 1 ) + G(α 1 ) -G(α 2 ) = 0,
where the second line follows from (4). Since Φ = λΦ c + (1λ)Φ d , Corollary 2.2 is proved.

Proof of Theorem 3.1.

The proof of Theorem 3.1 follows a classical strategy and consists of two steps. Lemma 6.1. Under the hypotheses of Theorem 3.1, there is a random variable V such that, E |V| 2 < ∞, and

lim n→∞ Z n -m 2 = V, a.s. Proof of Lemma 6.1. Let us consider V n := Z n -m 2 . Recall that Z n+1 = Z n -γ n Φ(Z n ) + γ n ξ n+1 (cf. (11)). Therefore V n+1 = Z n -m -γ n Φ(Z n ) 2 + γ 2 n ξ n+1 2 + 2γ n ξ n+1 , Z n -γ n Φ(Z n ) .
If we condition with respect to F n , the last term disappears since (ξ n ) is a martingale difference sequence and it comes:

E [V n+1 |F n ] = Z n -m -γ n Φ(Z n ) 2 + γ 2 n E ξ n+1 2 F n = Z n -m 2 -2γ n Z n -m, Φ(Z n ) + γ 2 n Φ(Z n ) 2 + E ξ n+1 2 F n = V n -2γ n Z n -m, Φ(Z n ) + γ 2 n , (16)
where we used the definition of V n and (12) for the last term. Since G is convex, using Corollary 2.2, we get:

Z n -m, Φ(Z n ) = Z n -m, Φ(Z n ) -Φ(m) ≥ 0. Therefore, for all n, E [V n+1 |F n ] ≤ V n + γ 2 n .
From the Robbins Siegmund theorem (see for instance (Duflo, 1997, page 18)), we deduce that (V n ) converges almost surely to V. Moreover, we note that

Z n -m is bounded in L 2 , ∀n, V n = E Z n -m 2 ≤ E Z 0 -m 2 + ∞ ∑ k=1 γ 2 k < ∞, (17)
whenever

E Z 0 -m 2 < ∞, which is satisfied for example if Z 0 = X 0 1 { X 0 ≤M} , with M < ∞.
We can now give the proof the theorem.

Proof of Theorem 3.1. Lemma 6.1 shows that the sequence V n converges almost surely. Let us check now that its limit is zero. Let us take expectations in equation ( 16):

E [V n+1 ] = E [V n ] + γ 2 n -2γ n E [ Φ(Z n ), Z n -m ] = E [V 0 ] + n ∑ k=1 γ 2 k -2 n ∑ k=1 γ k E [ Φ(Z k ), Z k -m ] . The sequence ∑ n k=1 γ k E [ Φ(Z k ), Z k -m ]
has positive terms, and is bounded above by

E [V 0 ] + ∑ ∞ k=1 γ 2 k , therefore it converges. This implies in particular that ∞ ∑ n=1 γ n Φ(Z n ), Z n -m < +∞ a.s. ( 18
)
This convergence cannot happen unless Z n converges to m. Indeed, for each ǫ ∈]0, 1[, let us introduce the set

Ω ǫ = ω ∈ Ω : ∃n ǫ (ω) ≥ 1, ∀n ≥ n ǫ (ω), ǫ 2 < V n (ω) < ǫ -2 .
For ω ∈ Ω ǫ , we have with Corollary 2.2,

∑ n≥1 γ n Φ(Z n (ω)), Z n (ω) -m ≥ ∑ n≥n ǫ (ω) γ n inf ǫ< α-m <ǫ -1 Φ(α), α -m = ∞,
which contradicts (18) unless P(Ω ǫ ) = 0. Since V n converges a.s. to a finite limit, and

{lim V n ∈ [c, c -1 ]} ⊂ Ω c/2
, the only possible limit is zero: lim Z nm = 0, a.s.

Proof of Lemma 3.3 and Proposition 3.2

Proof of Lemma 3.3. Consider, for h ∈ B(0, A), the function

f h (t) = Φ(m + th), defined for t ∈ [0, 1]. We have f h (0) = Φ(m) = 0 and f h (1) = Φ(m + h). It is also clear that the first order derivative f ′ h (t) of function f h satisfies f ′ h (t) = Γ m+th .
Consequently, a Taylor expansion with integral remainder of f h about t = 0 gives us

Φ(m + h) = Φ(m) + 1 0 Γ m+th (h) dt.
By Lemma 5.7 in [START_REF] Chaudhuri | Multivariate location estimation using extension of R-estimates through U-statistics type approach[END_REF], there is a constant M A such that for all t ∈ [0, 1],

Γ m+th -Γ m L ≤ M A h
where . L is the usual norm for bounded linear operators. Since Φ(m) = 0, one gets

Φ(m + h) -Γ m (h) ≤ sup t∈[0,1] Γ m+th -Γ m L h ≤ M A h 2 ,
and this concludes the proof.

Proof of Proposition 3.2. The proof is composed of 5 steps.

Step 1 -a spectral decomposition. Recall that Γ m is:

Γ m = E 1 X -m I H - (X -m) ⊗ (X -m) X -m 2 = E X -m -1 I H -E 1 X -m (X -m) ⊗ (X -m) X -m 2 = E X -m -1 I H -∆ m . (19)
Since Γ m is bounded and symmetric, it is self-adjoint. Moreover, the operator ∆ m defined by ( 19) is trace class: it is self-adjoint, non negative, and if (e j ) is an orthonormal basis,

∑ j e j , ∆ m e j = ∑ j E X -m, e j 2 X -m 3 ≤ E 1 X -m < ∞.
Therefore ∆ m is compact, and there is an increasing sequence of eigenvalues (λ j ), with possible repetitions, and an orthonormal basis (v j ) of eigenvectors in H such that:

∀j ∈ N, Γ m v j = λ j v j , σ(Γ m ) = {λ j , j ∈ N} ∪ E X -m -1 , λ j j→∞ --→ E X -m -1 .
Moreover, thanks to (7), the smallest eigenvalue λ min of Γ m is strictly positive. For simplicity of notation, we rewrite this decomposition as follows,

Γ m x = ∑ λ∈Λ λ e λ , x e λ , x ∈ H,
where Λ is the multiset {λ j , j ∈ N}, that can account for eigenspaces of dimension larger than 1.

In the following, we will need the operators:

α k = I H -γ k Γ m , β n = α n α n-1 • • • α 1 .
Since Γ m is bounded, these operators are well defined. Introducing the sequence of real functions, for n ∈ N,

f n (x) = n ∏ k=1 (1 -γ k x),
we see that f n (•) and f -1 n (•) are well defined on σ(Γ m ), provided γ n E Xm -1 < 1, which we can assume without loss of generality. Elementary analysis shows that there exist constants c 1 , C 2 , C 3 such that:

∀x ∈ σ(Γ m ), c 1 exp (-s n x) ≤ f n (x) ≤ C 2 exp (-s n x) , s n - c γ 1 -α n 1-α ≤ C 3 , (20)
where we recall that s n = ∑ n k=1 γ k , and γ k = c γ k -α . Then each operator β n can be also expressed as follows:

β n x = ∑ λ∈Λ f n (λ) e λ , x e λ , x ∈ H,
their inverses are bounded operators, and satisfy:

β -1 n x = ∑ λ∈Λ f -1 n (λ) e λ , x e λ .
Step 2 -Decomposition of the algorithm. Let us rewrite the algorithm in the following way

Z n+1 = Z n + γ n ξ n+1 -γ n Φ(Z n ) = Z n + γ n ξ n+1 -γ n (Γ m (Z n -m) + δ n ) where δ n = Φ(Z n ) -Γ m (Z n -m)
is the difference between the gradient of G and the gradient of its quadratic approximation. Therefore:

∀k, Z k+1 -m = α k (Z k -m) + γ k ξ k+1 -γ k δ k (21) Rewriting α n-1 α n-2 • • • α k+1 as β n-1 β -1 k , we get by induction, Z n -m = β n-1 (Z 1 -m) + β n-1 M n -β n-1 R n-1 , (22)
where

R n = n-1 ∑ k=1 γ k β -1 k δ k M n = n-1 ∑ k=1 γ k β -1 k ξ k+1 .
The first two terms of ( 22) are what we would get if G was exactly quadratic: a deterministic gradient part going to m, and a noise part; R n is the error term. We will look at each of these terms in turn.

Step 3 -The deterministic term. We want to bound β n-1 (Z 1m). The asymptotic behaviour of f n in eq. ( 20) implies that

β n-1 ≤ C 2 exp (-s n λ min ) ,
where λ min > 0 is the smallest eigenvalue of Γ m . Therefore

E β n-1 (Z 1 -m) 2 ≤ C exp -2n 1-α E Z 1 -m 2 . ( 23
)
Step 4 -The martingale. The fact that the β k are operators (instead of real numbers) makes matters more complicated. To deal with this problem, we use the spectral decomposition of the sequence of self-adjoint operators (β k ).

More precisely, we decompose M n = ∑ λ∈Λ e λ , M n e λ = ∑ λ M λ n e λ . For each λ ∈ Λ, M λ n is a martingale, and

E[(M λ n ) 2 ] = ∑ k≤n-1 γ 2 k f -2 k (λ)E ξ k+1 , e λ 2 F k ,
since E [ ξ k ′ , e λ ξ k+1 , e λ |F k ] = 0 when k ′ < k + 1. Summing now over λ ∈ Λ, we get:

E β n-1 M n 2 = ∑ λ f 2 n-1 (λ)E (M λ n ) 2 , ≤ ∑ λ ∑ k≤n-1 γ 2 k f n-1 (λ) f k (λ) 2 E ξ k+1 , e λ 2 .
However, for any k, n, and any λ ∈ Λ,

f n-1 (λ) f k (λ) = n-1 ∏ j=k+1 1 -λγ j ≤ f n-1 (λ min ) f k (λ min ) .
This uniformity in λ allows us to reconstruct E ξ k+1 2 , which is bounded by 1, thanks to (12). We obtain:

E β n-1 M n 2 ≤ ∑ k≤n-1 γ 2 k f n-1 (λ min ) f k (λ min ) 2 ∑ λ E ξ k+1 , e λ 2 , ≤ ∑ k≤n-1 γ 2 k f n-1 (λ min ) f k (λ min ) 2 E ξ k+1 2 , ≤ ∑ k≤n-1 γ 2 k f n-1 (λ min ) f k (λ min ) 2 .
Now we use the bounds (20) on β n :

E β n-1 M n 2 ≤ C 2 2 c 2 1 ∑ k≤n-1 γ 2 k exp - n ∑ j=k+1 γ j ≤ C ∑ k≤n-1 γ 2 k exp - 1 1 -α n 1-α -k 1-α . ( 24
)
The exponential terms are very small when k is much smaller than n, therefore we isolate the last terms. To do that, we choose l(n) such that,

l(n) 1-α = n 1-α -c α ln(n) , (25)
with c α to be chosen later. The first part of the sum (24) (for k ≤ l(n)) gives us:

∑ k≤l(n) γ 2 k exp - 1 1 -α n 1-α -k 1-α ≤ ∑ k≤l(n) γ 2 k exp - c α 1 -α ln(n) ≤ c 2 γ ∑ k≤n k -2α-cα 1-α . ( 26
)
This can be made smaller than any prescribed inverse power of n, if we choose c α large enough. In the second part of the sum (24), for k > l(n), we bound the exponential by 1 and γ k by γ l(n) :

∑ k>l(n) γ 2 k exp - 1 1 -α n 1-α -k 1-α ≤ (n -l(n))γ 2 l(n) .
The number of terms nl(n) is equivalent to c α 1-α ln(n)n α , and γ l(n) ∼ c γ n -α . Therefore, the whole second term is equivalent to c ln(n)n -α , where c depends on c α and c γ . For c α large enough, this dominates the first term (26). Finally we get:

E β n-1 M n 2 ≤ C ln(n) n α . ( 27
)
Step 5 -the error term and the conclusion. The error term is

β n-1 ∑ n k=1 γ k β -1 k δ k , where δ k = Φ(Z k ) -Γ m (Z k -m). With Lemma 3.3, we get that ∃r, C r ∀k, Z k -m ≤ r =⇒ δ k ≤ C r Z k -m 2 . ( 28
)
Since Z n converges a.s. to m, we deduce two things about δ k : it is almost surely bounded, and (28) becomes a.s. eventually true. To use these facts we introduce the following sequence of events:

Ω N =      ω, ∀n ≥ N, ∀k ≥ n -l(n), Z k (ω) -m ≤ 1/K and δ k (ω) ≤ C r Z k (ω) -m 2 ∀k, δ k (ω) ≤ N.      ,
for a value of K to be chosen later, and l(n) defined by ( 25). This sequence is increasing and Ω N = Ω; from now on we work on Ω N .

Once more, since β n-1 β -1 k is very small when k is much smaller than n, only the last terms in the sum defining R n matter. This is why we re-use the definition of l(n) and cut the sum in two parts. For ω ∈ Ω N , and n ≥ N,

β n-1 R n 2 ≤ n ∑ k=1 γ k β n-1 β -1 k δ k 2 ≤ 2N 2 l(n) ∑ k=1 γ k β n-1 β -1 k 2 + 2C 2 r n ∑ k=l(n)+1 γ k Z k -m 2 2 ≤ 2N 2 l(n) ∑ k=1 γ k β n-1 β -1 k 2 + 2 C 2 r K 2 (n -l(n))γ l(n) n ∑ k=l(n)+1 γ k Z k -m 2 .
where we used the crude bound δ k ≤ N in the first part, and for the second part,

β n-1 β -1 k ≤ 1 and the definition of Ω N .
As before, it is easy to see that the first term is bounded by any prescribed inverse power of n, say n -42 . For the second term, we already know that (nl(n))γ l(n) is bounded. Therefore, on Ω N and for n ≥ N,

β n-1 R n 2 ≤ CN 2 n 42 + C K 2 n ∑ k=l(n)+1 γ k Z k -m 2 . ( 29
)
Combining now ( 22), ( 23), ( 27) and (29), we get, for n ≥ N and some new constant

C E 1 Ω N Z n -m 2 ≤ C ln(n) n α + C K 2 n ∑ k=l(n)+1 γ k E 1 Ω N Z k -m 2 ≤ C ln(n) n α + C ′ K 2 sup l(n)<k≤n E 1 Ω N Z k -m 2 . Let us choose K such that K 2 ≥ 2C ′ . Then ∀n ≥ N, E 1 Ω N Z n -m 2 ≤ C ln(n) n α + 1 2 max l(n)<k≤n E 1 Ω N Z k -m 2 . Defining u n = E 1 Ω N Z n -m 2 , this reads: ∀n ≥ N, u n ≤ C ln(n) n α + 1 2 max l(n)<k≤n u k . ( 30
)
Let us prove by induction that, for some N ′ large enough, and for

C ′′ = 4C, ∀n ≥ N ′ , u n ≤ C ′′ ln(n) n α . Suppose that u k ≤ C ′′ ln(k) k α
for all k ≤ n, and let us prove that u n+1 ≤ C ′′ ln(n+1) (n+1) α . Using (30), we know that:

u n+1 ≤ C ln(n + 1) (n + 1) α + 1 2 max l(n+1)<k≤n+1 u k ≤ C ln(n + 1) (n + 1) α + 1 2 max u n+1 , C ′′ ln(l(n + 1)) l(n + 1) α
If the max on the right hand side is u n+1 , we get: 1 2 u n+1 ≤ C ln(n + 1) (n + 1) α which is the desired result since 2C ≤ 4C = C ′′ . If this is not the case, then the max is

C ′′ ln(l(n+1)) l(n+1) α
. However, l(n + 1) ∼ (n + 1) so for n larger than some N ′ , ln(l(n+1 This concludes the induction step and the proof of Proposition 3.2.

Proof of Theorem 3.4

We use the same decomposition as in [START_REF] Pelletier | Asymptotic almost sure efficiency of averaged stochastic algorithms[END_REF].

Γ m T k = ξ k+1 -δ k + 1 γ k (T k -T k+1 ) .
Summing these equalities, it comes,

nΓ m T n = n ∑ k=1 1 γ k (T k -T k+1 ) - n ∑ k=1 δ k + M n+1 .
Applying Abel's transform, and dividing by √ n yields:

√ nΓ m T n = 1 √ n T 1 γ 1 - T n+1 γ n + n ∑ k=2 T k 1 γ k - 1 γ k-1 - n ∑ k=1 δ k + 1 √ n M n+1 .
To prove that last term is a martingale for which the CLT holds,

M n √ n L ---→ n→∞ N (0, Σ) ,
we need to check that the assumptions of Theorem 5.1 in [START_REF] Jakubowski | Tightness criteria for random measures with application to the principle of conditioning in Hilbert spaces[END_REF] are fulfilled.

We first have that the martingale difference sequence is a.s. bounded, ∀n ξ n ≤ 2. Let us define

Σ n = E [ξ n+1 ⊗ ξ n+1 |F n ] ,
which can also be decomposed as follows

Σ n = E (X -Z n ) X -Z n ⊗ (X -Z n ) X -Z n F n -Φ(Z n ) ⊗ Φ(Z n ).
Since Φ(m) = 0, we have by a direct computation,

Φ(Z n ) ≤ E 2 X -m Z n -m .
Using now, for (a, b) ∈ H × H, the inequality a ⊗ b L ≤ a b , where a ⊗ b L is the usual the norm for linear operators, we directly get, with Theorem 3.1,

Φ(Z n ) ⊗ Φ(Z n ) L → 0, a.s.
With similar arguments, it is easy to show that

Σ -E (X -Z n ) X -Z n ⊗ (X -Z n ) X -Z n F n L ≤ 2E (X -Z n ) X -Z n - (X -m) X -m F n ≤ 4E 1 X -m Z n -m ,
so that Σ n -Σ L → 0 a.s., when n tends to infinity. Then condition 5.2 in [START_REF] Jakubowski | Tightness criteria for random measures with application to the principle of conditioning in Hilbert spaces[END_REF] is satisfied and is a consequence of a direct application of Chow's Lemma, see for instance (Duflo, 1997, page 22). Now, it remains to prove that

1 √ n T n+1 γ n - n ∑ k=2 T k 1 γ k - 1 γ k-1 + n ∑ k=1 δ k P ---→ n→∞ 0. ( 31
)
Let us denote by A n , A ′ n and A ′′ n the three terms. Recall that E 1

Ω N T n 2 ≤ C N ln(n)
n α , thanks to Proposition 3.2. For the first term A n =

T n+1 √ nγ n
, we have:

E 1 Ω N A n 2 ≤ C ′ N n 2α-1 ln(n) 2 n 2α = C ′ N ln(n) 2 n .
Therefore A n P ---→ n→∞ 0.

Let us turn to the second term A ′ n . Since γ -1 kγ -1 k-1 ≤ 2αc -1 γ k α-1 , we have, for two positive constants C 0 , C 1 ,

E A ′ n 1 Ω N ≤ 2αc -1 γ √ n ∑ k≤n E [1 Ω N T k ] k α-1 ≤ C 0 √ n ∑ k≤n ln(k)k α/2-1
≤ C 1 ln(n)n α/2-1/2 , which goes to zero since α < 1. Therefore A ′ Finally, for the last term A ′′ n , since there exists a positive constant C 2 such that δ k ≤ C 2 Z km 2 , we have:

E 1 Ω N A ′′ n ≤ 1 √ n ∑ k≤n E 1 Ω N T k 2 ≤ C N √ n ∑ k≤n ln(k)k -α .
Since the right hand side term converges to zero (as can be seen e.g. by Kronecker's lemma, using the fact that α > 1/2), C n P ---→ n→∞ 0, therefore (31) holds, and Theorem 3.4 is finally proved.

Figure 1 :

 1 Figure 1: Estimations of the mean and the geometric median audiences, at a second scale, during the 6th september 2010.

  It consists in linearizing the target function Φ around the true value m. Recall the following decomposition of the error (21),∀k, Z k+1m = (I Hγ k Γ m )(Z km) + γ k ξ k+1γ k δ k ,where ξ k is a martingale difference sequence and δ k are error terms,δ k = Φ(Z k ) -Γ m (Z km). Defining now, T n := Z nm, T n := Z nm and M n+1 := n ∑ k=1 ξ k+1 ,and rearranging the previous expression, we obtain:

Table 1 :

 1 Comparison of the estimation errors for different sample sizes

			n=250		n=500		n=2000
	Estimator	[Q1 median Q3] [Q1 median Q3] [Q1 median Q3]
	c γ = 0.2	0.45	0.60	0.80 0.38	0.53	0.69 0.25	0.35	0.47
	c γ = 0.6	0.21	0.29	0.40 0.15	0.21	0.29 0.06	0.09	0.12
	c γ = 1	0.15	0.22	0.31 0.11	0.16	0.21 0.05	0.08	0.10
	c γ = 2	0.15	0.21	0.30 0.09	0.15	0.20 0.05	0.07	0.10
	c γ = 5	0.13	0.19	0.25 0.09	0.13	0.18 0.04	0.06	0.09
	c γ = 10	0.13	0.18	0.25 0.09	0.13	0.18 0.04	0.06	0.09
	c γ = 15	0.12	0.18	0.25 0.09	0.13	0.18 0.04	0.06	0.08
	c γ = 25	0.13	0.19	0.26 0.09	0.13	0.18 0.04	0.06	0.09
	c γ = 50	0.13	0.19	0.26 0.09	0.13	0.18 0.04	0.06	0.09
	c γ = 75	0.14	0.20	0.27 0.09	0.14	0.19 0.05	0.07	0.09
	Vardi & Zhang 0.12	0.18	0.25 0.09	0.12	0.17 0.04	0.06	0.08

Acknowledgements. We would like to thank the referees for their helpful and valuable suggestions. We also thank the company Médiamétrie for allowing us to illustrate our methodologies with their data.