Keywords: AMS 2000 subject classifications: 62J05, 62J07 Linear regression, Generalization error, Shrinkage, PAC-Bayesian theorems, Risk bounds, Robust statistics, Resistant estimators, Gibbs posterior distributions, Randomized estimators, Statistical learning theory

ROBUST LINEAR LEAST SQUARES REGRESSION

By Jean-Yves Audibert * , † , and Olivier Catoni ‡, §

We consider the problem of robustly predicting as well as the best linear combination of d given functions in least squares regression, and variants of this problem including constraints on the parameters of the linear combination. For the ridge estimator and the ordinary least squares estimator, and their variants, we provide new risk bounds of order d/n without logarithmic factor unlike some standard results, where n is the size of the training data. We also provide a new estimator with better deviations in presence of heavy-tailed noise. It is based on truncating differences of losses in a min-max framework and satisfies a d/n risk bound both in expectation and in deviations. The key common surprising factor of these results is the absence of exponential moment condition on the output distribution while achieving exponential deviations. All risk bounds are obtained through a PAC-Bayesian analysis on truncated differences of losses. Experimental results strongly back up our truncated min-max estimator. Our statistical task. Let Z 1 = (X 1 , Y 1 ), . . . , Z n = (X n , Y n ) be n ≥ 2 pairs of input-output and assume that each pair has been independently drawn from the same unknown distribution P . Let X denote the input space and let the output space be the set of real numbers R, so that P is a probability distribution on the product space Z X × R. The target of learning algorithms is to predict the output Y associated with an input X for pairs Z = (X, Y ) drawn from the distribution P . The quality of a (prediction) function f : X → R is measured by the least squares risk:

R(f ) E Z∼P [Y -f (X)] 2 .
Through the paper, we assume that the output and all the prediction functions we consider are square integrable. Let Θ be a closed convex set of R d , and ϕ 1 , . . . , ϕ d be d prediction functions. Consider the regression model

F = f θ = d j=1 θ j ϕ j ; (θ 1 , . . . , θ d ) ∈ Θ .
The best function f * in F is defined by

f * = d j=1 θ * j ϕ j ∈ argmin f ∈F R(f ).
Such a function always exists but is not necessarily unique. Besides it is unknown since the probability generating the data is unknown. We will study the problem of predicting (at least) as well as function f * . In other words, we want to deduce from the observations Z 1 , . . . , Z n a function f having with high probability a risk bounded by the minimal risk R(f * ) on F plus a small remainder term, which is typically of order d/n up to a possible logarithmic factor. Except in particular settings (e.g., Θ is a simplex and d ≥ √ n), it is known that the convergence rate d/n cannot be improved in a minimax sense (see [START_REF] Tsybakov | Optimal rates of aggregation[END_REF], and [START_REF] Yang | Aggregating regression procedures for a better performance[END_REF] for related results). More formally, the target of the paper is to develop estimators f for which the excess risk is controlled in deviations, i.e., such that for an appropriate constant κ > 0, for any ε > 0, with probability at least 1ε,

(1.1) R( f ) -R(f * ) ≤ κ d + log(ε -1 ) n .
Note that by integrating the deviations (using the identity E(W ) = +∞ 0 P(W > t)dt which holds true for any nonnegative random variable W ), Inequality (1.1) implies

(1.2) ER( f ) -R(f * ) ≤ κ (d + 1) n .
In this work, we do not assume that the function

f (reg) : x → E[Y |X = x],
which minimizes the risk R among all possible measurable functions, belongs to the model F. So we might have f * = f (reg) and in this case, bounds of the form

(1.3) ER( f ) -R(f (reg) ) ≤ C R(f * ) -R(f (reg) ) + κ d n ,
with a constant C larger than 1 do not even ensure that ER( f ) tends to R(f * ) when n goes to infinity. This kind of bounds with C > 1 have been developed to analyze nonparametric estimators using linear approximation spaces, in which case the dimension d is a function of n chosen so that the bias term R(f * ) -R(f (reg) ) has the order d/n of the estimation term (see [START_REF] Györfi | A Distribution-Free Theory of Nonparametric Regression[END_REF] and references within). Here we intend to assess the generalization ability of the estimator even when the model is misspecified (namely when R(f * ) > R(f (reg) )). Moreover we do not assume either that Yf (reg) (X) and X are independent or that Y has a subexponential tail distribution: for the moment, we just assume that Yf * (X) admits a finite second order moment in order that the risk of f * is finite. Several risk bounds with C = 1 can be found in the literature. A survey on these bounds is given in [1, Section 1]. Let us mention here the closest bound to what we are looking for. From the work of Birgé and Massart [START_REF] Birgé | Minimum contrast estimators on sieves: exponential bounds and rates of convergence[END_REF], we may derive the following risk bound for the empirical risk minimizer on a L ∞ ball (see Appendix B of [START_REF] Audibert | Robust linear regression through PAC-Bayesian truncation[END_REF]).

Theorem 1.1. Assume that F has a diameter H for L ∞ -norm, i.e., for any

f 1 , f 2 in F, sup x∈X |f 1 (x) -f 2 (x)| ≤ H
and there exists a function f 0 ∈ F satisfying the exponential moment condition:

(1.4) for any x ∈ X , E exp A -1 Y -f 0 (X) X = x ≤ M,
for some positive constants A and M . Let

B = inf φ 1 ,...,φ d sup θ∈R d -{0} d j=1 θ j φ j 2 ∞ θ 2 ∞
where the infimum is taken with respect to all possible orthonormal bases of

F for the dot product (f 1 , f 2 ) → E f 1 (X)f 2 (X)
(when the set F admits no basis with exactly d functions, we set B = +∞). Then the empirical risk minimizer satisfies for any ε > 0, with probability at least 1ε:

R( f (erm) ) -R(f * ) ≤ κ(A 2 + H 2 ) d log[2 + ( B/n) ∧ (n/d)] + log(ε -1 ) n ,
where κ is a positive constant depending only on M .

The theorem gives exponential deviation inequalities of order at worse d log(n/d)/n, and asymptotically, when n goes to infinity, of order d/n. This work will provide similar results under weaker assumptions on the output distribution.

Notation. When Θ = R d , the function f * and the space F will be written f * lin and F lin to emphasize that F is the whole linear space spanned by ϕ 1 , . . . , ϕ d :

F lin = span{ϕ 1 , . . . , ϕ d } and f * lin ∈ argmin f ∈F lin R(f ).
The Euclidean norm will simply be written as • , and •, • will be its associated inner product. We will consider the vector valued function ϕ :

X → R d defined by ϕ(X) = ϕ k (X) d k=1
, so that for any θ ∈ Θ, we have f θ (X) = θ, ϕ(X) .

The Gram matrix is the d × d-matrix Q = E ϕ(X)ϕ(X) T . The empirical risk of a function f is r(f ) = 1 n n i=1 f (X i ) -Y i 2
and for λ ≥ 0, the ridge regression estimator on F is defined by f

(ridge) = f θ(ridge) with θ(ridge) ∈ arg min θ∈Θ r(f θ ) + λ θ 2 ,
where λ is some nonnegative real parameter. In the case when λ = 0, the ridge regression f (ridge) is nothing but the empirical risk minimizer f (erm) .

Besides, the empirical risk minimizer when Θ = R d is also called the ordinary least squares estimator, and will be denoted by f (ols) .

In the same way, we introduce the optimal ridge function optimizing the expected ridge risk: f = f θ with

(1.5) θ ∈ arg min θ∈Θ R(f θ ) + λ θ 2 .
Finally, let Q λ = Q + λI be the ridge regularization of Q, where I is the identity matrix.

Why should we be interested in this task. There are four main reasons. First we intend to provide a non-asymptotic analysis of the parametric linear least squares method. Secondly, the task is central in nonparametric estimation for linear approximation spaces (piecewise polynomials based on a regular partition, wavelet expansions, trigonometric polynomials. . . )

Thirdly, it naturally arises in two-stage model selection. Precisely, when facing the data, the statistician has often to choose several models which are likely to be relevant for the task. These models can be of similar structure (like embedded balls of functional spaces) or on the contrary of very different nature (e.g., based on kernels, splines, wavelets or on a parametric approach). For each of these models, we assume that we have a learning scheme which produces a 'good' prediction function in the sense that it predicts as well as the best function of the model up to some small additive term. Then the question is to decide on how we use or combine/aggregate these schemes. One possible answer is to split the data into two groups, use the first group to train the prediction function associated with each model, and finally use the second group to build a prediction function which is as good as (i) the best of the previously learnt prediction functions, (ii) the best convex combination of these functions or (iii) the best linear combination of these functions. This point of view has been introduced by Nemirovski in [START_REF] Nemirovski | Topics in non-parametric statistics. Ecole d'Eté de Probabilités de Saint-Flour XXVIII-1998[END_REF] and optimal rates of aggregation are given in [START_REF] Tsybakov | Optimal rates of aggregation[END_REF] and references within. This paper focuses more on the linear aggregation task (even if (ii) enters in our setting), assuming implicitly here that the models are given in advance and are beyond our control and that the goal is to combine them appropriately.

Finally, in practice, the noise distribution often departs from the normal distribution. In particular, it can exhibit much heavier tails, and consequently induce highly non Gaussian residuals. It is then natural to ask whether classical estimators such as the ridge regression and the ordinary least squares estimator are sensitive to this type of noise, and whether we can design more robust estimators.

Outline and contributions. Section 2 provides a new analysis of the ridge estimator and the ordinary least squares estimator, and their variants. Theorem 2.1 provides an asymptotic result for the ridge estimator while Theorem 2.2 gives a non asymptotic risk bound for the empirical risk minimizer, which is complementary to the theorems put in the survey section. In particular, the result has the benefit to hold for the ordinary least squares estimator and for heavy-tailed outputs. We show quantitatively that the ridge penalty leads to an implicit reduction of the input space dimension. Section 3 shows a non asymptotic d/n exponential deviation risk bound under weak moment conditions on the output Y and on the d-dimensional input representation ϕ(X).

The main contribution of this paper is to show through a PAC-Bayesian analysis on truncated differences of losses that the output distribution does not need to have bounded conditional exponential moments in order for the excess risk of appropriate estimators to concentrate exponentially. Our results tend to say that truncation leads to more robust algorithms. Local robustness to contamination is usually invoked to advocate the removal of outliers, claiming that estimators should be made insensitive to small amounts of spurious data. Our work leads to a different theoretical explanation. The observed points having unusually large outputs when compared with the (empirical) variance should be down-weighted in the estimation of the mean, since they contain less information than noise. In short, huge outputs should be truncated because of their low signal-to-noise ratio.

Ridge regression and empirical risk minimization. We recall the definition

F = f θ = d j=1 θ j ϕ j ; (θ 1 , . . . , θ d ) ∈ Θ ,
where Θ is a closed convex set, not necessarily bounded (so that Θ = R d is allowed). In this section, we provide exponential deviation inequalities for the empirical risk minimizer and the ridge regression estimator on F under weak conditions on the tail of the output distribution.

The most general theorem which can be obtained from the route followed in this section is Theorem 1.5 of the supplementary material [START_REF] Audibert | Supplement to "Robust linear least squares regression[END_REF]. It is expressed in terms of a series of empirical bounds. The first deduction we can make from this technical result is of asymptotic nature. It is stated under weak hypotheses, taking advantage of the weak law of large numbers.

Theorem 2.1. For λ ≥ 0, let f be its associated optimal ridge function (see (1.5)). Let us assume that

E ϕ(X) 4 < +∞, (2.1) and E ϕ(X) 2 f (X) -Y 2 < +∞. (2.2) Let ν 1 > • • • > ν d be the eigenvalues of the Gram matrix Q = E ϕ(X)ϕ(X) T , and let Q λ = Q + λI be the ridge regularization of Q. Let us define the ef- fective ridge dimension D = d i=1 ν i ν i + λ 1(ν i > 0) = Tr (Q + λI) -1 Q = E Q -1/2 λ ϕ(X) 2 .
When λ = 0, D is equal to the rank of Q and is otherwise smaller. For any ε > 0, there is n ε , such that for any n ≥ n ε , with probability at least

1 -ε, R( f (ridge) ) + λ θ(ridge) 2 ≤ min θ∈Θ R(f θ ) + λ θ 2 + 30 E Q -1/2 λ ϕ(X) 2 f (X) -Y 2 E Q -1/2 λ ϕ(X) 2 D n + 1000 sup v∈R d E v, ϕ(X) 2 f (X) -Y 2 E v, ϕ(X) 2 + λ v 2 log 3ε -1 n ≤ min θ∈Θ R(f θ ) + λ θ 2 + ess sup E [Y -f (X)] 2 X 30D + 1000 log 3ε -1 n .
Proof. See Section 1 of the supplementary material [START_REF] Audibert | Supplement to "Robust linear least squares regression[END_REF].

This theorem shows that the ordinary least squares estimator (obtained when Θ = R d and λ = 0), as well as the empirical risk minimizer on any closed convex set, asymptotically reaches a d/n speed of convergence under very weak hypotheses. It shows also the regularization effect of the ridge regression. There emerges an effective dimension D, where the ridge penalty has a threshold effect on the eigenvalues of the Gram matrix.

Let us remark that the second inequality stated in the theorem provides a simplified bound which makes sense only when

ess sup E Y -f (X) 2 | X < +∞, implying that f -f (reg) ∞ < +∞.
We chose to state the first inequality as well, since it does not require such a tight relationship between f and f (reg) . On the other hand, the weakness of this result is its asymptotic nature: n ε may be arbitrarily large under such weak hypotheses, and this happens even in the simplest case of the estimation of the mean of a real valued random variable by its empirical mean (which is the case when d = 1 and ϕ(X) ≡ 1).

Let us now give some non asymptotic rate under stronger hypotheses and for the empirical risk minimizer (i.e., λ = 0).

Theorem 2.2. Assume that E [Y -f * (X)] 4 < +∞ and B = sup f ∈span{ϕ 1 ,...,ϕ d }-{0} f 2 ∞ /E[f (X) 2 ] < +∞.
Consider the (unique) empirical risk minimizer

f (erm) = f θ(erm) : x → θ(erm) , ϕ(x) on F for which θ(erm) ∈ span{ϕ(X 1 ), . . . , ϕ(X n )} 1 .
For any values of ε and n such that 2/n ≤ ε ≤ 1 and

n > 1280B 2 3Bd + log(2/ε) + 16B 2 d 2 n ,
with probability at least 1ε,

(2.3) R( f (erm) ) -R(f * ) ≤ 1920 B E [Y -f * (X)] 4 3Bd + log(2ε -1 ) n + 4Bd n 2 . 1 When F = F lin , we have θ(erm) = X + Y, with X = (ϕj(Xi)) 1≤i≤n,1≤j≤d , Y = [Yj] n j=1
and X + is the Moore-Penrose pseudoinverse of X.

Proof. See Section 1 of the supplementary material [START_REF] Audibert | Supplement to "Robust linear least squares regression[END_REF].

It is quite surprising that the traditional assumption of uniform boundedness of the conditional exponential moments of the output can be replaced by a simple moment condition for reasonable confidence levels (i.e., ε ≥ 2/n). For highest confidence levels, things are more tricky since we need to control with high probability a term of order [r(f * ) -R(f * )]d/n (see Theorem 1.6). The cost to pay to get the exponential deviations under only a fourth-order moment condition on the output is the appearance of the geometrical quantity B as a multiplicative factor.

To better understand the quantity B, let us consider two cases. First, consider that the input is uniformly distributed on X = [0, 1], and that the functions ϕ 1 , . . . , ϕ d belong to the Fourier basis. Then the quantity B behaves like a numerical constant. On the contrary, if we take ϕ 1 , . . . , ϕ d as the first d elements of a wavelet expansion, the more localized wavelets induce high values of B, and B scales like √ d, meaning that Theorem 2.2 fails to give a d/n-excess risk bound in this case. This limitation does not appear in Theorem 2.1.

To conclude, Theorem 2.2 is limited in at least four ways: it involves the quantity B, it applies only to uniformly bounded ϕ(X), the output needs to have a fourth moment, and the confidence level should be as great as ε ≥ 2/n. These limitations will be addressed in the next section by considering a more involved algorithm.

3. A min-max estimator for robust estimation. This section provides an alternative to the empirical risk minimizer with non asymptotic exponential risk deviations of order d/n for any confidence level. Moreover, we will assume only a second order moment condition on the output and cover the case of unbounded inputs, the requirement on ϕ(X) being only a finite fourth order moment. On the other hand, we assume here that the set Θ of the vectors of coefficients is bounded. The computability of the proposed estimator and numerical experiments are discussed at the end of the section.

3.1. The min-max estimator and its theoretical guarantee. Let α > 0, λ ≥ 0, and consider the truncation function:

ψ(x) =      -log 1 -x + x 2 /2 0 ≤ x ≤ 1, log(2)
x ≥ 1, -ψ(-x)

x ≤ 0, For any θ, θ ′ ∈ Θ, introduce

D(θ, θ ′ ) = nαλ θ 2 -θ ′ 2 + n i=1 ψ α Y i -f θ (X i ) 2 -α Y i -f θ ′ (X i ) 2 .
We recall that f = f θ with θ ∈ arg min θ∈Θ R(f θ ) + λ θ 2 , and that the effective ridge dimension is defined as

D = E Q -1/2 λ ϕ(X) 2 = Tr (Q + λI) -1 Q = d i=1 ν i ν i + λ 1(ν i > 0) ≤ d,
where

ν 1 ≥ • • • ≥ ν d are the eigenvalues of the Gram matrix Q = E ϕ(X)ϕ(X) T .
Let us assume in this section that

(3.1) E [Y -f (X)] 4 < +∞,
and that for any j ∈ {1, . . . , d},

(3.2) E ϕ j (X) 4 < +∞. Define S = f ∈ F lin : E f (X) 2 = 1 , (3.3) σ = E [Y -f (X)] 2 = R( f ), (3.4) χ = max f ∈S E f (X) 4 , (3.5) κ = E [ϕ(X) T Q -1 λ ϕ(X)] 2 E ϕ(X) T Q -1 λ ϕ(X) , (3.6) κ ′ = E [Y -f (X)] 4 E [Y -f (X)] 2 = E [Y -f (X)] 4 σ 2 , (3.7) T = max θ∈Θ,θ ′ ∈Θ λ θ -θ ′ 2 + E [f θ (X) -f θ ′ (X)]
(3.9) α = 1 2χ 2 √ κ ′ σ + √ χT 2 1 - cκχD n ,
for any estimator f θ satisfying θ ∈ Θ a.s., for any ε > 0 and any λ ≥ 0, with probability at least 1ε, we have

R(f θ) + λ θ 2 ≤ min θ∈Θ R(f θ ) + λ θ 2 + 1 nα max θ 1 ∈Θ D( θ, θ 1 ) -inf θ∈Θ max θ 1 ∈Θ D(θ, θ 1 ) + cκκ ′ Dσ 2 n + 8χ log(ε -1 ) n + c ′ κ 2 D 2 n 2 2 √ κ ′ σ + √ χT 2 1 - cκχD n .
Proof. See Section 2 of the supplementary material [START_REF] Audibert | Supplement to "Robust linear least squares regression[END_REF].

By choosing an estimator such that max

θ 1 ∈Θ D( θ, θ 1 ) < inf θ∈Θ max θ 1 ∈Θ D(θ, θ 1 ) + σ 2 D n ,
Theorem 3.1 provides a non asymptotic bound for the excess (ridge) risk with a D/n convergence rate and an exponential tail even when neither the output Y nor the input vector ϕ(X) have exponential moments. This stronger non asymptotic bound compared to the bounds of the previous section comes at the price of replacing the empirical risk minimizer by a more involved estimator. Section 3.3 provides a way of computing it approximately. Theorem 3.1 requires a fourth order moment condition on the output. In fact, one can replace (3.1) by the following second order moment condition on the output: for any j ∈ {1, . . . , d},

E ϕ j (X) 2 [Y -f (X)] 2 < +∞,
and still obtain a D/n excess risk bound. This comes at the price of a more lengthy formula, where terms with κ ′ become terms involving the quantities 2 . This can be seen by not using Cauchy-Schwarz's inequality in (2.5)and (2.6).

max f ∈S E f (X) 2 [Y -f (X)] 2 and E ϕ(X) T Q -1 ϕ(X)[Y -f (X)]

3.2.

The value of the uncentered kurtosis coefficients χ and κ. We see that the speed of convergence of the excess risk in Theorem 3.1 (page 10) depends on three kurtosis like coefficients, χ, κ and κ ′ . The third, κ ′ , is concerned with the noise, conceived as the difference between the observed output Y and its best explanation f (X) according to the ridge criterion. The aim of this section is to study the order of magnitude of the two other coefficients χ and κ, which are related to the design distribution.

χ = sup E u, ϕ(X) 4 1/2 ; u ∈ R d , E u, ϕ(X) 2 ≤ 1 , and κ = D -1 E Q -1/2 λ ϕ(X) 4 1/2 .
We will review a few typical situations.

3.2.1. Gaussian design. Let us assume first that ϕ(X) is a multivariate centered Gaussian random variable. In this case, its covariance matrix coincides with its Gram matrix Q 0 and can be written as

Q 0 = U -1 Diag(ν i , i = 1, . . . , n) U,
where U is an orthogonal matrix. Using U , we can introduce

W = U Q -1/2 λ ϕ(X).
It is also a Gaussian vector, with covariance

Diag ν i /(λ + ν i ), i = 1, . . . , d . Moreover, since U is orthogonal, W = Q -1/2 λ
ϕ(X) , and since (W i , W j ) are uncorrelated when i = j, they are independent, leading to

E Q -1/2 λ ϕ(X) 4 = E d i=1 W 2 i 2 = d i=1 E W 4 i + 2 1≤i<j≤d E W 2 i E W 2 j = D 2 + 2D 2 ,
where

D 2 = d i=1 ν 2 i λ + ν i 2 .
Thus, in this case

κ = 1 + 2D 2 D -2 ≤ 1 + 2ν 1 (λ + ν 1 )D ≤ √ 3.
Moreover, as for any value of u, u, ϕ(X) is a Gaussian random variable, χ = √ 3. This situation arises in compressed sensing using random projections on Gaussian vectors. Specifically, assume that we want to recover a signal f ∈ R M that we know to be well approximated by a linear combination of d basis vectors f 1 , . . . , f d . We measure n ≪ M projections of the signal f on i.i.d. M -dimensional standard normal random vectors X 1 , . . . , X n : Y i = f, X i , i = 1, . . . , n. Then, recovering the coefficient θ 1 , . . . , θ d such that f = d j=1 θ j f j is associated to the least squares regression problem Y ≈ d j=1 θ j ϕ j (X), with ϕ j (x) = f j , x , and X having a M -dimensional standard normal distribution.

3.2.2. Independent design. Let us study now the case when almost surely ϕ 1 (X) ≡ 1 and ϕ 2 (X), . . . , ϕ d (X) are independent. To compute χ, we can assume without loss of generality that ϕ 2 (X), . . . , ϕ d (X) are centered and of unit variance, since this renormalization is precisely the linear transformation that turns the Gram matrix into the identity matrix. Let us introduce

χ * = max j=1,...,d E ϕ j (X) 4 1/2 E ϕ j (X) 2 ,
with the convention 0 0 = 0. A computation similar to the one made in the Gaussian case shows that

κ ≤ 1 + (χ 2 * -1)D 2 D -2 ≤ 1 + (χ 2 * -1)ν 1 (λ + ν 1 )D ≤ χ * .
Moreover, for any

u ∈ R d such that u = 1, E u, ϕ(X) 4 = d i=1 u 4 i E(ϕ i (X) 4 ) + 6 1≤i<j≤d u 2 i u 2 j E ϕ i (X) 2 E ϕ j (X) 2 + 4 d i=2 u 1 u 3 i E ϕ i (X) 3 ≤ χ 2 * d i=1 u 4 i + 6 i<j u 2 i u 2 j + 4χ 3/2 * d i=2 |u 1 u i | 3 ≤ sup u∈R d + , u =1 χ 2 * -3 d i=1 u 4 i + 3 d i=1 u 2 i 2 + 4χ 3/2 * u 1 d i=2 u 3 i ≤ 3 3/2 4 χ 3/2 * + χ 2 * , χ 2 * ≥ 3, 3 + χ 2 * -3 d , 1 ≤ χ 2 * < 3. Thus in this case χ ≤      χ * 1 + 3 3/2 4 √ χ * 1/2 , χ * ≥ √ 3, 3 + 3 3/2 4 χ 3/2 * + χ 2 * -3 d 1/2 , 1 ≤ χ * < √ 3.
If moreover the random variables ϕ 2 (X), . . . , ϕ d (X) are not skewed, in the sense that E ϕ j (X) 3 = 0, j = 2, . . . , d, then

   χ = χ * , χ * ≥ √ 3, χ ≤ 3 + χ 2 * -3 d 1/2 , 1 ≤ χ * < √ 3.
3.2.3. Bounded design. Let us assume now that the distribution of ϕ(X) is almost surely bounded and nearly orthogonal. These hypotheses are suited to the study of regression in usual function bases, like the Fourier basis, wavelet bases, histograms, or splines.

More precisely, let us assume that P( ϕ(X) ≤ B) = 1 and that for some positive constant A and any u ∈ R d ,

u ≤ AE u, ϕ(X) 2 1/2 .
This appears as some stability property of the partial basis ϕ j with respect to the L 2 -norm, since it can also be written as

d j=1 u 2 j ≤ A 2 E d j=1 u j ϕ j (X) 2 , u ∈ R d .
In terms of eigenvalues, A -2 can be taken to be the lowest eigenvalue ν d of the Gram matrix Q. The value of A can also be deduced from a condition saying that ϕ j are nearly orthogonal in the sense that

E ϕ j (X) 2 ≥ 1, and E ϕ j (X)ϕ k (X) ≤ 1 -A -2 d -1 .
In this situation, the chain of inequalities

E u, ϕ(X) 4 ≤ u 2 B 2 E u, ϕ(X) 2 ≤ A 2 B 2 E u, ϕ(X) 2 2
shows that χ ≤ AB. On the other hand,

E Q -1/2 λ ϕ(X) 4 = E sup u, ϕ(X) 4 ; u ∈ R d , Q 1/2 λ u ≤ 1 ≤ E sup u 2 B 2 u, ϕ(X) 2 ; Q 1/2 λ u ≤ 1 ≤ E sup 1 + λA 2 -1 A 2 B 2 Q 1/2 λ u 2 u, ϕ(X) 2 ; Q 1/2 λ u ≤ 1 ≤ A 2 B 2 1 + λA 2 E Q -1/2 λ ϕ(X) 2 = A 2 B 2 D 1 + λA 2 ,
showing that κ ≤ AB (1 + λA 2 )D .

For example, if X is the uniform random variable on the unit interval and ϕ j , j = 1, . . . , d are any functions from the Fourier basis (meaning that they are of the form √ 2 cos(2kπX) or √ 2 sin(2kπX)), then A = 1 (because they form an orthogonal system) and B ≤ √ 2d.

A localized basis like the evenly spaced histogram basis of the unit interval

ϕ j (x) = √ d1 x ∈ (j -1)/d, j/d , j = 1, . . . , d,
will also be such that A = 1 and B = √ d. Similar computations could be made for other local bases, like wavelet bases.

Note that when χ is of order √ d, and κ and κ ′ of order 1, Theorem 3.1 means that the excess risk of the min-max truncated estimator f is upper bounded by Cd/n provided that n ≥ Cd 2 for a large enough constant C.

3.2.4.

Adaptive design planning. Let us discuss the case when X is some observed random variable whose distribution is only approximately known. Namely let us assume that (ϕ j ) d j=1 is some basis of functions in L 2 P with some known coefficient χ, where P is an approximation of the true distribution of X in the sense that the density of the true distribution P of X with respect to the distribution P is in the range (η -1 , η). In this situation, the coefficient χ satisfies the inequality χ ≤ η 3/2 χ. Indeed

E X∼P u, ϕ(X) 4 ≤ ηE X∼ P u, ϕ(X) 4 ≤ η χ 2 E X∼ P u, ϕ(X) 2 2 ≤ η 3 χ 2 E X∼P u, ϕ(X) 2 2 .
In the same way, κ ≤ η 7/2 κ. Indeed,

E sup u, ϕ(X) 4 ; E u, ϕ(X) 2 ≤ 1 ≤ η E sup u, ϕ(X) 4 ; E u, ϕ(X) 2 ≤ η ≤ η 3 E sup u, ϕ(X) 4 ; E u, ϕ(X) 2 ≤ 1 ≤ η 3 κ 2 E sup u, ϕ(X) 2 ; E u, ϕ(X) 2 ≤ 1 2 ≤ η 7 κ 2 E sup u, ϕ(X) 2 ; E u, ϕ(X) 2 ≤ 1 2 .
Let us conclude this section with some scenario for the case when X is a real-valued random variable. Let us consider the distribution function of P F (x) = P(X ≤ x).

Then, if P has no atoms, the distribution of F (X) would be uniform on (0, 1) if X were distributed according to P. In other words, P • F -1 = U, the uniform distribution on the unit interval. Starting from some suitable partial basis (ϕ j ) d j=1 of L 2 (0, 1), U like the ones discussed above, we can build a basis for our problem as

ϕ j (X) = ϕ j F (X) .
Moreover, if P is absolutely continuous with respect to P with density g, then P • F -1 is absolutely continuous with respect to P • F -1 = U, with density g • F -1 , and of course, the fact that g takes values in (η -1 , η) implies the same property for g • F -1 . Thus, if χ and κ are the coefficients corresponding to ϕ j (U ) when U is the uniform random variable on the unit interval, then the true coefficient χ (corresponding to ϕ j (X)) will be such that χ ≤ η 3/2 χ and κ ≤ η 7/2 κ.

3.3.

Computation of the estimator. For ease of description of the algorithm, we will write X for ϕ(X), which is equivalent to considering without loss of generality that the input space is R d and that the functions ϕ 1 , . . . ,ϕ d are the coordinate functions. Therefore, the function f θ maps an input x to θ, x . Let us introduce

L i (θ) = α θ, X i -Y i 2 .
For any subset of indices I ⊂ {1, . . . , n}, let us define

r I (θ) = λ θ 2 + 1 α|I| i∈I L i (θ).
We suggest the following heuristics to compute an approximation of arg min

θ∈Θ sup θ ′ ∈Θ D(θ, θ ′ ).
• Start from I 1 = {1, . . . , n} with the ordinary least squares estimate

θ 1 = arg min R d r I 1 .
• At step number k, compute

Q k = 1 |I k | i∈I k X i X T i .
• Consider the sets

J k,1 (η) = i ∈ I k : L i ( θ k )X T i Q -1 k X i 1 + 1 + L i ( θ k ) -1 2 < η ,
where Q -1 k is the (pseudo-)inverse of the matrix Q k .

• Let us define

θ k,1 (η) = arg min R d r J k,1 (η) , J k,2 (η) = i ∈ I k : L i θ k,1 (η) -L i θ k ≤ 1 , θ k,2 (η) = arg min R d r J k,2 (η) , (η k , ℓ k ) = arg min η∈R + ,ℓ∈{1,2} max j=1,...,k D θ k,ℓ (η), θ j , I k+1 = J k,ℓ k (η k ), θ k+1 = θ k,ℓ k (η k ).
• Stop when max j=1,...,k D( θ k+1 , θ j ) ≥ 0, and set θ = θ k as the final estimator of θ.

Note that there will be at most n steps, since I k+1 I k and in practice much less in this iterative scheme. Let us give some justification for this proposal.

Let us notice first that

D(θ + h, θ) = nαλ( θ + h 2 -θ 2 ) + n i=1 ψ α 2 h, X i θ, X i -Y i + h, X i 2 .
Hopefully, θ = arg min θ∈R d R(f θ ) + λ θ 2 is in some small neighbourhood of θ k already, according to the distance defined by Q ≃ Q k . So we may try to look for improvements of θ k by exploring neighbourhoods of θ k of increasing sizes with respect to some approximation of the relevant norm

θ 2 Q = E θ, X 2 .
Since the truncation function ψ is constant on (-∞, -1] and [1, +∞), the map θ → D(θ, θ k ) induces a decomposition of the parameter space into cells corresponding to different sets I of examples. Indeed, such a set I is associated to the set

C I of θ such that L i (θ) -L i ( θ k ) < 1 if and only if i ∈ I.
Although this may not be the case, we will do as if the map θ → D(θ, θ k ) restricted to the cell C I reached its minimum at some interior point of C I , and approximates this minimizer by the minimizer of r I .

The idea is to remove first the examples which will become inactive in the closest cells to the current estimate θ k . The cells for which the contribution of example number i is constant are delimited by at most four parallel hyperplanes.

It is easy to see that the square of the inverse of the distance of θ k to the closest of these hyperplanes is equal to

1 α X T i Q -1 k X i L i ( θ k ) 1 + 1 + 1 L i ( θ k ) 2 .
Indeed, this distance is the infimum of Q

1/2 k h , where h is a solution of h, X i 2 + 2 h, X i θ k , X i -Y i = 1 α .

It is computed by considering h of the form

h = ξ Q -1/2 k X i -1 Q -1
k X i and solving an equation of order two in ξ.

This explains the proposed choice of J k,1 (η). Then a first estimate θ k,1 (η) is computed on the basis of this reduced sample, and the sample is readjusted to J k,2 (η) by checking which constraints are really activated in the computation of D(θ k,1 (η), θ k ). The estimated parameter is then readjusted taking into account the readjusted sample (this could as a variant be iterated more than once). Now that we have some new candidates θ k,ℓ (η), we check the minimax property against them to elect I k+1 and θ k+1 . Since we did not check the minimax property against the whole parameter set Θ = R d , we have no theoretical warranty for this simplified algorithm. Nonetheless, similar computations to what we did could prove that we are close to solving min j=1,...,k R(f θ j ), since we checked the minimax property on the reduced parameter set { θ j , j = 1, . . . , k}. Thus the proposed heuristics is capable of improving on the performance of the ordinary least squares estimator, while being guaranteed not to degrade its performance significantly.

3.4. Synthetic experiments. In Section 3.4.1, we detail the different kinds of noises we work with. Then, Sections 3.4.2, 3.4.3 and 3.4.4 describe the three types of functional relationships between the input, the output and the noise involved in our experiments. A motivation for choosing these inputoutput distributions was the ability to compute exactly the excess risk, and thus to compare easily estimators. Section 3.4.5 provides details about the implementation, its computational efficiency and the main conclusions of the numerical experiments. • the standard Gaussian noise: W ∼ N (0, 1),

• a heavy-tailed noise defined by: W = sign(V )/|V | 1/q , with V ∼ N (0, 1) a standard Gaussian random variable and q = 2.01 (the real number q is taken strictly larger than 2 as for q = 2, the random variable W would not admit a finite second moment). • an asymmetric heavy-tailed noise defined by:

W = |V | -1/q if V > 0, -q q-1
otherwise, with q = 2.01 with V ∼ N (0, 1) a standard Gaussian random variable. • a mixture of a Dirac random variable with a low-variance Gaussian random variable defined by: with probability p, W = (1ρ)/p, and with probability 1p, W is drawn from

N - p(1 -ρ) 1 -p , ρ 1 -p - p(1 -ρ) (1 -p) 2 .
The parameter ρ ∈ [p, 1] characterizes the part of the variance of W explained by the Gaussian part of the mixture. Note that this noise admits exponential moments, but for n of order 1/p, the Dirac part of the mixture generates low signal-to-noise points.

Independent normalized covariates (INC(n, d))

. In INC(n, d), we consider ϕ(X) = X, and the input-output pair is such that

Y = θ * , X + σW,
where the components of X are independent standard normal distributions, θ * = (10, . . . , 10) T ∈ R d , and σ = 10.

Highly correlated covariates (HCC(n, d))

. In HCC(n, d), we consider ϕ(X) = X, and the input-output pair is such that

Y = θ * , X + σW,
where X is a multivariate centered normal Gaussian with covariance matrix Q obtained by drawing a (d, d)-matrix A of uniform random variables in [0, 1] and by computing Q = AA T , θ * = (10, . . . , 10) T ∈ R d , and σ = 10. So the only difference with the setting of Section 3.4.2 is the correlation between the covariates.

Trigonometric series (TS(n, d))

. Let X be a uniform random variable on [0, 1]. Let d be an even number. In TS(n, d), we consider ϕ(X) = cos(2πX), . . . , cos(dπX), sin(2πX), . . . , sin(dπX)

T , and the input-output pair is such that

Y = 20X 2 -10X - 5 3 + σW,
with σ = 10. One can check that this implies

θ * = 20 π 2 , . . . , 20 π 2 ( d 2 ) 2 , - 10 π , . . . , - 10 π( d 2 )
T ∈ R d .

Experiments.

Choice of the parameters and implementation details.. The min-max truncated algorithm has two parameters α and λ. In the subsequent experiments, we set the ridge parameter λ to the natural default choice for it: λ = 0. For the truncation parameter α, according to our analysis (see (3.9)), it roughly should be of order 1/σ 2 up to kurtosis coefficients. By using the ordinary least squares estimator, we roughly estimate this value, and test values of α in a geometric grid (of 8 points) around it (with ratio 3). Cross-validation can be used to select the final α. Nevertheless, it is computationally expensive and is significantly outperformed in our experiments by the following simple procedure: start with the smallest α in the geometric grid and increase it as long as θ = θ 1 , that is as long as we stop at the end of the first iteration and output the empirical risk minimizer.

To compute θ k,1 (η) or θ k,2 (η), one needs to determine a least squares estimate (for a modified sample). To reduce the computational burden, we do not want to test all possible values of η (note that there are at most n values leading to different estimates). Our experiments show that testing only three levels of η is sufficient. Precisely, we sort the quantity

L i ( θ k )X T i Q -1 k X i 1 + 1 + L i ( θ k ) -1 2
by decreasing order and consider η being the first, 5-th and 25-th value of the ordered list. Overall, in our experiments, the computational complexity is approximately fifty times larger than the one of computing the ordinary least squares estimator.

Results.. The tables and figures have been gathered in Appendix A. Tables 1 and2 give the results for the mixture noise. Tables 3, 4 and 5 provide the results for the heavy-tailed noise and the standard Gaussian noise. Each line of the tables has been obtained after 1000 generations of the training set. These results show that the min-max truncated estimator is often equal to the ordinary least squares estimator f (ols) , while it ensures impressive consistent improvements when it differs from f (ols) . In this latter case, the number of points that are not considered in f , i.e. the number of points with low signal-to-noise ratio, varies a lot from 1 to 150 and is often of order 30. Note that not only the points that we expect to be considered as outliers (i.e. very large output points) are erased, and that these points seem to be taken out by local groups: see Figures 1 and2 in which the erased points are marked by surrounding circles. Besides, the heavier the noise tail is (and also the larger the variance of the noise is), the more often the truncation modifies the initial ordinary least squares estimator, and the more improvements we get from the minmax truncated estimator, which also becomes much more robust than the ordinary least squares estimator (see the confidence intervals in the tables).

Finally, we have also tested more traditional methods in robust regression: namely, the M-estimators with Huber's loss, L 1 -loss, and Tukey's bisquare influence function, and also the least trimmed squares estimator, the Sestimator and the MM-estimator (see [START_REF] Rousseeuw | Robust regression by means of S-estimators[END_REF][START_REF] Yohai | High breakdown-point and high efficiency robust estimates for regression[END_REF] and references within). These methods rely on diminishing the influence of points having "unreasonably" large residuals. They were developed to handle training sets containing true outliers, i.e. points (X, Y ) not generated by the distribution P . This is not the case in our estimation framework. By overweighting points having reasonably small residuals, these methods are often biased even in settings where the noise is symmetric and the regression function f (reg) : x → E[Y |X = x] belongs to F lin (i.e., f (reg) = f * lin ), and also even when there is no noise (but f (reg) / ∈ f * lin ). The worst results were obtained by the L 1 -loss, since estimating the (conditional) median is here really different from estimating the (conditional) mean. The MM-estimator and the M-estimators with Huber's loss and Tukey's bisquare influence function give good results as long as the signal-to-noise ratio is low. When the signal-to-noise ratio is high, a lack of consistency drastically appears in part of our simulations, showing that these methods are thus not suited for our estimation framework.

The S-estimator is almost consistently improving on the ordinary least squares estimator (in our simulations). However, when the signal-to-noise ratio is low (that is, in the setting of the aforementioned simulations with σ = 10), the improvements are much less significant than the ones of the min-max truncated estimator.

4. Main ideas of the proofs. The goal of this section is to explain the key ingredients appearing in the proofs which both allow to obtain subexponential tails for the excess risk under a non-exponential moment assumption and get rid of the logarithmic factor in the excess risk bound.

4.1. Sub-exponential tails under a non-exponential moment assumption via truncation. Let us start with the idea allowing us to prove exponential inequalities under just a moment assumption (instead of the traditional exponential moment assumption). To understand it, we can consider the (apparently) simplistic 1-dimensional situation in which we have Θ = R and the marginal distribution of ϕ 1 (X) is the Dirac distribution at 1. In this case, the risk of the prediction function 2 , so that the least squares regression problem boils down to the estimation of the mean of the output variable. If we only assume that Y admits a finite second moment, say E Y 2 ≤ 1, it is not clear whether for any ε > 0, it is possible to find θ such that with probability at least 1 -2ε,

f θ is R(f θ ) = E (Y -θ) 2 = E (Y -EY ) 2 + (EY -θ)
(4.1) R(f θ) -R(f * ) = E(Y ) -θ 2 ≤ c log(ε -1 ) n ,
for some numerical constant c. Indeed, from Chebyshev's inequality, the trivial choice θ = 1 n n i=1 Y i just satisfies: with probability at least 1 -2ε,

R(f θ) -R(f * ) ≤ 1 nε ,
which is far from the objective (4.1) for small confidence levels (consider ε = exp(-√ n) for instance). The key idea is thus to average (soft) truncated values of the outputs. This is performed by taking

θ = 1 nλ n i=1 log 1 + λY i + λ 2 Y 2 i 2 , with λ = 2 log(ε -1 ) n . Since we have log E exp(nλ θ) = n log 1 + λE(Y ) + λ 2 2 E(Y 2 ) ≤ nλE(Y ) + n λ 2 2 ,
the exponential Chebyshev's inequality guarantees that with probability at least 1ε, we have nλ

θ -E(Y ) ≤ nλ 2 /2 + log(ε -1 ), hence θ -E(Y ) ≤ 2 log(ε -1 ) n .
Replacing Y by -Y in the previous argument, we obtain that with probability at least 1ε, we have

nλ E(Y ) + 1 nλ n i=1 log 1 -λY i + λ 2 Y 2 i 2 ≤ n λ 2 2 + log(ε -1 ). Since -log(1 + x + x 2 /2) ≤ log(1 -x + x 2 /2), this implies E(Y ) -θ ≤ 2 log(ε -1 ) n
. The two previous inequalities imply Inequality (4.1) (for c = 2), showing that sub-exponential tails are achievable even when we only assume that the random variable admits a finite second moment (see [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF] for more details on the robust estimation of the mean of a random variable).

Localized PAC-Bayesian inequalities to eliminate a logarithm factor.

Let us first recall that the Kullback-Leibler divergence between distributions ρ and µ defined on F is

(4.2) K(ρ, µ)    E f ∼ρ log dρ dµ (f ) if ρ ≪ µ, +∞ otherwise,
where dρ dµ denotes as usual the density of ρ w.r.t. µ. For any real-valued (measurable) function h defined on F such that exp[h(f )]π(df ) < +∞, we define the distribution π h on F by its density:

(4.3) dπ h dπ (f ) = exp[h(f )] exp[h(f ′ )]π(df ′ )
.

The analysis of statistical inference generally relies on upper bounding the supremum of an empirical process χ indexed by the functions in a model F. Concentration inequalities appear as a central tool to obtain these bounds. An alternative approach, called the PAC-Bayesian one, consists in using the entropic equality

(4.4) E exp sup ρ∈M ρ(df )χ(f ) -K(ρ, π ′ ) = π ′ (df )E exp χ(f ) .
where M is the set of probability distributions on F. Let ř : F → R be an observable process such that for any f ∈ F, we have

E exp χ(f ) ≤ 1 for χ(f ) = λ[R(f ) -ř(f )
] and some λ > 0. Then (4.4) leads to: for any ε > 0, with probability at least 1ε, for any distribution ρ on F, we have

(4.5) ρ(df )R(f ) ≤ ρ(df )ř(f ) + K(ρ, π ′ ) + log(ε -1 ) λ .
The lefthand-side quantity represents the expected risk with respect to the distribution ρ. To get the smallest upper bound on this quantity, a natural choice of the (posterior) distribution ρ is obtained by minimizing the righthand-side, that is by taking ρ = π ′ -λř (with the notation introduced in (4.3)). This distribution concentrates on functions f ∈ F for which ř(f ) is small. Without prior knowledge, one may want to choose a prior distribution π ′ = π which is rather "flat" (e.g., the one induced by the Lebesgue measure in the case of a model F defined by a bounded parameter set in some Euclidean space). Consequently the Kullback-Leibler divergence K(ρ, π ′ ), which should be seen as the complexity term, might be excessively large.

To overcome the lack of prior information and the resulting high complexity term, one can alternatively use a more "localized" prior distribution. Here we use Gaussian distributions centered at the function of interest (for instance, the function f * ), and with covariance matrix proportional to the inverse of the Gram matrix Q. The idea of using PAC-Bayesian inequalities with Gaussian prior and posterior distributions goes back to Langford and Shawe-Taylor [START_REF] Langford | PAC-bayes & margins[END_REF] in the context of linear classification.

The detailed proofs of Theorems 2.1, 2.2 and 3.1 can be found in the supplementary material [START_REF] Audibert | Supplement to "Robust linear least squares regression[END_REF].

SUPPLEMENTARY MATERIAL

To be filled by production: (http://lib.stat.cmu.edu/aoas/???/???). The supplementary material provides the proofs of Theorems 2.1, 2.2 and 3.1. )

ER( f (ols) ) -R(f * ) ER( f ) -R(f * ) ER[( f (ols) )| f = f (ols) ] -R(f * ) E[R( f )| f = f (ols) ] -R(f * ) INC(n=200,

CONTENTS 1 2 13 3. 2 . 3

 121323 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Our statistical task . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 Why should we be interested in this task . . . . . . . . . . . . . . . 5 Outline and contributions . . . . . . . . . . . . . . . . . . . . . . .6 Ridge regression and empirical risk minimization . . . . . . . . . . 6 3 A min-max estimator for robust estimation . . . . . . . . . . . . . 9 3.1 The min-max estimator and its theoretical guarantee . . . . . 9 3.2 The value of the uncentered kurtosis coefficients χ and κ . . . 11 3.2.1 Gaussian design . . . . . . . . . . . . . . . . . . . . . 12 3.2.2 Independent design . . . . . . . . . . . . . . . . . . . . Bounded design . . . . . . . . . . . . . . . . . . . . . . 14 3.2.4 Adaptive design planning . . . . . . . . . . . . . . . . 15 3.3 Computation of the estimator . . . . . . . . . . . . . . . . . . 16 3.4 Synthetic experiments . . . . . . . . . . . . . . . . . . . . . . 18 3.4.1 Noise distributions . . . . . . . . . . . . . . . . . . . . 18 3.4.2 Independent normalized covariates (INC(n, d)) . . . . 19 3.4.3 Highly correlated covariates (HCC(n, d)) . . . . . . . . 19 3.4.4 Trigonometric series (TS(n, d)) . . . . . . . . . . . . . 20 3.4.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . 20 4 Main ideas of the proofs . . . . . . . . . . . . . . . . . . . . . . . . 22 4.1 Sub-exponential tails under a non-exponential moment assumption via truncation . . . . . . . . . . . . . . . . . . . . . 22 4.2 Localized PAC-Bayesian inequalities to eliminate a logarithm factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 Supplementary Material . . . . . . . . . . . . . . . . . . . . . . . . . . 24 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 A Experimental results for the min-max truncated estimator (Section 3.3) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1. Introduction.

Theorem 3 . 1 .

 31 Let us assume that (3.1) and (3.2) hold. For some numerical constants c and c ′ , for n > cκχD, by taking

  Figures and tables are postponed to Appendix A. 3.4.1. Noise distributions. In our experiments, we consider different types of noise that are centered and with unit variance:

Fig 1 .Fig 2 .

 12 Fig 1.Circled points are the points of the training set generated several times from TS(1000, 10) (with the mixture noise with p = 0.005 and ρ = 0.4) that are not taken into account in the min-max truncated estimator (to the extent that the estimator would not change by removing simultaneously all these points). The min-max truncated estimator x → f (x) appears in dash-dot line, while x → E(Y |X = x) is in solid line. In these six simulations, it outperforms the ordinary least squares estimator.

Table 3

 3 Comparison of the min-max truncated estimator f with the ordinary least squares estimator f (ols) with the heavy-tailed noise (see Section 3.4.1).

		) = R( f (ols)
	nb of iterations	nb of iter. with R( f

) nb of iter. with R( f

) < R( f (ols)

APPENDIX A: EXPERIMENTAL RESULTS FOR THE MIN-MAX TRUNCATED ESTIMATOR (SECTION 3.3)

Table 1 Comparison of the min-max truncated estimator f with the ordinary least squares estimator f (ols) for the mixture noise (see Section 3.4.1) with ρ = 0.1 and p = 0.005. In parenthesis, the 95%-confidence intervals for the estimated quantities.

nb of iterations nb of iter. with R( f

) nb of iter. with R( f Comparison of the min-max truncated estimator f with the ordinary least squares estimator f (ols) for the mixture noise (see Section 3.4.1) with ρ = 0.4 and p = 0.005. In parenthesis, the 95%-confidence intervals for the estimated quantities.

nb of iterations nb of iter. with R( f ) nb of iter. with R( f


Jean-Yves Audibert
email: audibert@imagine.enpc.fr

Olivier Catoni
email: olivier.catoni@ens.fr

SUPPLEMENT TO "ROBUST LINEAR LEAST SQUARES REGRESSION"

Keywords: AMS 2000 subject classifications: 62J05, 62J07 Linear regression, Generalization error, Shrinkage, PAC-Bayesian theorems, Risk bounds, Robust statistics, Resistant estimators, Gibbs posterior distributions, Randomized estimators, Statistical learning theory

1. Proofs of Theorems 2.1 and 2.2. To shorten the formulae, we will write X for ϕ(X), which is equivalent to considering without loss of generality that the input space is R d and that the functions ϕ 1 , . . . ,ϕ d are the coordinate functions. Therefore, the function f θ maps an input x to θ, x . With a slight abuse of notation, R(θ) will denote the risk of this prediction function.

Let us first assume that the matrix Q λ = Q + λI is positive definite. This indeed does not restrict the generality of our study, even in the case when λ = 0, as we will discuss later (Remark 1.1).

Consider the change of coordinates

X = Q -1/2 λ X. Let us introduce R(θ) = E ( θ, X -Y ) 2 , so that R(Q 1/2 λ θ) = R(θ) = E ( θ, X -Y ) 2 . Let Θ = Q 1/2 λ θ; θ ∈ Θ . Consider r(θ) = 1 n n i=1 θ, X i -Y i 2 , (1.1) r(θ) = 1 n n i=1 θ, X i -Y i 2 , (1.2) θ 0 = arg min θ∈Θ R(θ) + λ Q -1/2 λ θ 2 , (1.3) θ ∈ arg min θ∈Θ r(θ) + λ θ 2 , (1.4) θ 1 = Q 1/2 λ θ ∈ arg min θ∈Θ r(θ) + λ Q -1/2 λ θ 2 .
(1.5) For α > 0, let us introduce the notation

W i (θ) = α θ, X i -Y i 2 -θ 0 , X i -Y i 2 , W (θ) = α θ, X -Y 2 -θ 0 , X -Y 2 .
For any θ 2 ∈ R d and β > 0, let us consider the Gaussian distribution centered at θ 2

ρ θ 2 (dθ) = β 2π d/2 exp - β 2 θ -θ 2 2 dθ.
Lemma 1.1. For any η > 0 and α > 0, with probability at least 1exp(-η), for any

θ 2 ∈ R d , -n log 1 -E W (θ) + E W (θ) 2 /2 ρ θ 2 (dθ) ≤ - n i=1 log 1 -W i (θ) + W i (θ) 2 /2 ρ θ 2 (dθ) + K(ρ θ 2 , ρ θ 0 ) + η,
where K(ρ θ 2 , ρ θ 0 ) is the Kullback-Leibler divergence function :

K(ρ θ 2 , ρ θ 0 ) = log dρ θ 2 dρ θ 0 (θ) ρ θ 2 (dθ).
Proof. Since E ρ θ 0 (dθ)

n i=1 1 -W i (θ) + W i (θ) 2 /2 1 -E W (θ) + E W (θ) 2 /2 ≤ 1, with probability at least 1 -exp(-η) log ρ θ 0 (dθ) n i=1 1 -W i (θ) + W i (θ) 2 /2 1 -E W (θ) + E W (θ) 2 /2 ≤ η.
We conclude the proof using the convex inequality (see [START_REF] Donsker | Asymptotic evaluation of certain Markov process expectations for large time, I[END_REF], [3, Proposition

1.4.2] or [1, page 159]) log ρ θ 0 (dθ) exp h(θ) ≥ ρ θ 2 (dθ)h(θ) -K(ρ θ 2 , ρ θ 0 ).
Let us compute some useful quantities

K(ρ θ 2 , ρ θ 0 ) = β 2 θ 2 -θ 0 2 , (1.6) ρ θ 2 (dθ) W (θ) = α ρ θ 2 (dθ) θ -θ 2 , X 2 + W (θ 2 ) = W (θ 2 ) + α X 2 β , (1.7) ρ θ 2 (dθ) θ -θ 2 , X 4 = 3 X 4 β 2 , (1.8) (1.9) ρ θ 2 (dθ) W (θ) 2 = α 2 ρ θ 2 (dθ) θ -θ 0 , X 2 θ + θ 0 , X -2Y 2 = α 2 ρ θ 2 (dθ) θ -θ 2 + θ 2 -θ 0 , X θ -θ 2 + θ 2 + θ 0 , X -2Y 2 = ρ θ 2 (dθ) α θ -θ 2 , X 2 + 2α θ -θ 2 , X θ 2 , X -Y + W (θ 2 ) 2 = ρ θ 2 (dθ) α 2 θ -θ 2 , X 4 + 4α 2 θ -θ 2 , X 2 θ 2 , X -Y 2 + W (θ 2 ) 2 + 2α θ -θ 2 , X 2 W (θ 2 ) = 3α 2 X 4 β 2 + 2α X 2 β 2α θ 2 , X -Y 2 + W (θ 2 ) + W (θ 2 ) 2 .
Using the fact that

2α θ 2 , X -Y 2 + W (θ 2 ) = 2α θ 0 , X -Y 2 + 3W (θ 2 ),
and that for any real numbers a and b, 6ab ≤ 9a 2 + b 2 , we get Lemma 1.2.

ρ θ 2 (dθ) W (θ) = W (θ 2 ) + α X 2 β ,
(1.10)

ρ θ 2 (dθ) W (θ) 2 = W (θ 2 ) 2 + 2α X 2 β 2α θ 0 , X -Y 2 + 3W (θ 2 ) + 3α 2 X 4 β 2 (1.11) ≤ 10W (θ 2 ) 2 + 4α 2 X 2 β θ 0 , X -Y 2 + 4α 2 X 4 β 2 , (1.12)
and the same holds true when W is replaced with W i and (X, Y ) with

(X i , Y i ).
Another important thing to realize is that

E X 2 = E Tr X X T = E Tr Q -1/2 λ XX T Q -1/2 λ = E Tr Q -1 λ XX T = Tr Q -1 λ E(XX T ) = Tr Q -1 λ (Q λ -λI) = d -λ Tr(Q -1 λ ) = D . (1.13)
We can weaken Lemma 1.1 (page 2) noticing that for any real number x,

x - x 2 2 ≤ -log 1 -x + x 2 2 = log 1 + x + x 2 /2 1 + x 4 /4 ≤ log 1 + x + x 2 2 ≤ x + x 2 2 .
We obtain with probability at least 1 -exp(-η)

nE W (θ 2 ) + nα β E X 2 -5nE W (θ 2 ) 2 -E 2nα 2 X 2 β θ 0 , X -Y 2 + 2nα 2 X 4 β 2 ≤ n i=1 W i (θ 2 ) + 5W i (θ 2 ) 2 + α X i 2 β + 2α 2 X i 2 β θ 0 , X i -Y 2 + 2α 2 X i 4 β 2 + β 2 θ 2 -θ 0 2 + η.
Noticing that for any real numbers a and b, 4ab ≤ a 2 + 4b 2 , we can then bound

α -2 W (θ 2 ) 2 = θ 2 -θ 0 , X 2 θ 2 + θ 0 , X -2Y 2 = θ 2 -θ 0 , X 2 θ 2 -θ 0 , X + 2 θ 0 , X -Y 2 = θ 2 -θ 0 , X 4 + 4 θ 2 -θ 0 , X 3 θ 0 , X -Y + 4 θ 2 -θ 0 , X 2 θ 0 , X -Y 2 ≤ 2 θ 2 -θ 0 , X 4 + 8 θ 2 -θ 0 , X 2 θ 0 , X -Y 2 .
Theorem 1.3. Let us put

D = 1 n n i=1 X i 2 (let us remind that D = E X 2 from (1.13)), B 1 = 2E X 2 θ 0 , X -Y 2 , B 1 = 2 n n i=1 X i 2 θ 0 , X i -Y i 2 , B 2 = 2E X 4 , B 2 = 2 n n i=1 X i 4 , B 3 = 40 sup E u, X 2 θ 0 , X -Y 2 : u ∈ R d , u = 1 , B 3 = sup 40 n n i=1 u, X i 2 θ 0 , X i -Y i 2 : u ∈ R d , u = 1 , B 4 = 10 sup E u, X 4 : u ∈ R d , u = 1 , B 4 = sup 10 n n i=1 u, X i 4 : u ∈ R d , u = 1 .
With probability at least 1 -exp(-η), for any

θ 2 ∈ R d , nE W (θ 2 ) -nα 2 (B 3 + B 3 ) + β 2 θ 2 -θ 0 2 -nα 2 (B 4 + B 4 ) θ 2 -θ 0 4 ≤ n i=1 W i (θ 2 ) + nα β ( D -D) + nα 2 β (B 1 + B 1 ) + nα 2 β 2 (B 2 + B 2 ) + η.
Let us now assume that θ 2 ∈ Θ and let us use the fact that Θ is a convex set and that θ 0 = arg min θ∈Θ R(θ)

+ λ Q -1/2 λ θ 2 . Introduce θ * = arg min θ∈R d R(θ) + λ Q -1/2 λ θ 2 . As we have R(θ) + λ Q -1/2 λ θ 2 = θ -θ * 2 + R(θ * ) + λ Q -1/2 λ θ * 2 ,
the vector θ 0 is uniquely defined as the projection of θ * on Θ for the Euclidean distance, and for any

θ 2 ∈ Θ (1.14) α -1 E W (θ 2 ) + λ Q -1/2 λ θ 2 2 -λ Q -1/2 λ θ 0 2 = R(θ 2 ) -R(θ 0 ) + λ Q -1/2 λ θ 2 2 -λ Q -1/2 λ θ 0 2 = θ 2 -θ * 2 -θ 0 -θ * 2 = θ 2 -θ 0 2 + 2 θ 2 -θ 0 , θ 0 -θ * ≥ θ 2 -θ 0 2 .
This and the inequality

α -1 n i=1 W i (θ 1 ) + nλ Q -1/2 λ θ 1 2 -nλ Q -1/2 λ θ 0 2 ≤ 0
leads to the following result.

Theorem 1.4. With probability at least 1 -exp(-η),

R( θ) + λ θ 2 -inf θ∈Θ R(θ) + λ θ 2 = α -1 E W (θ 1 ) + λ Q -1/2 λ θ 1 2 -λ Q -1/2 λ θ 0 2
is not greater than the smallest positive non degenerate root of the following polynomial equation as soon as it has one

1 -α(B 3 + B 3 ) + β 2nα x -α(B 4 + B 4 )x 2 = 1 β max( D -D, 0) + α β (B 1 + B 1 ) + α β 2 (B 2 + B 2 ) + η nα .
Proof. Let us remark first that when the polynomial appearing in the theorem has two distinct roots, they are of the same sign, due to the sign of its constant coefficient. Let Ω be the event of probability at least 1-exp(-η) described in Theorem 1.3 (page 5). For any realization of this event for which the polynomial described in Theorem 1.4 does not have two distinct positive roots, the statement of Theorem 1.4 is void, and therefore fulfilled. Let us consider now the case when the polynomial in question has two distinct positive roots x 1 < x 2 . Consider in this case the random (trivially nonempty) closed convex set

Θ = θ ∈ Θ : R(θ) + λ θ 2 ≤ inf θ ′ ∈Θ R(θ ′ ) + λ θ ′ 2 + x 1 +x 2 2 . Let θ 3 ∈ arg min θ∈ Θ r(θ) + λ θ 2 and θ 4 ∈ arg min θ∈Θ r(θ) + λ θ 2 . We see from Theorem 1.3 that (1.15) R(θ 3 ) + λ θ 3 2 < R(θ 0 ) + λ θ 0 2 + x 1 + x 2 2 ,
because it cannot be larger from the construction of Θ. On the other hand, since Θ ⊂ Θ, the line segment

[θ 3 , θ 4 ] is such that [θ 3 , θ 4 ]∩ Θ ⊂ arg min θ∈ Θ r(θ)+ λ θ 2 . We can therefore apply equation (1.15) to any point of [θ 3 , θ 4 ] ∩ Θ, which proves that [θ 3 , θ 4 ] ∩ Θ is an open subset of [θ 3 , θ 4 ].
But it is also a closed subset by construction, and therefore, as it is non empty and

[θ 3 , θ 4 ] is connected, it proves that [θ 3 , θ 4 ] ∩ Θ = [θ 3 , θ 4 ],
and thus that θ 4 ∈ Θ. This can be applied to any choice of θ 3 ∈ arg min θ∈ Θ r(θ) + λ θ 2 and θ 4 ∈ arg min θ∈Θ r(θ) + λ θ 2 , proving that arg min θ∈Θ r(θ) + λ θ 2 ⊂ arg min θ∈ Θ r(θ) + λ θ 2 and therefore that any θ 4 ∈ arg min θ∈Θ r(θ) + λ θ 2 is such that

R(θ 4 ) + λ θ 4 2 ≤ inf θ∈Θ R(θ) + λ θ 2 + x 1 .
because the values between x 1 and x 2 are excluded by Theorem 1.3.

The actual convergence speed of the least squares estimator θ on Θ will depend on the speed of convergence of the "empirical bounds" B k towards their expectations. We can rephrase the previous theorem in the following more practical way: Theorem 1.5. Let η 0 , η 1 , . . . , η 5 be positive real numbers. With probability at least

1 -P D > D + η 0 - 4 k=1 P B k -B k > η k -exp(-η 5 ), R( θ) + λ θ 2 -inf θ∈Θ R(θ) + λ θ 2 is smaller than the smallest non de- generate positive root of (1.16) 1 -α(2B 3 + η 3 ) + β 2nα x -α(2B 4 + η 4 )x 2 = η 0 β + α β (2B 1 + η 1 ) + α β 2 (2B 2 + η 2 ) + η 5 nα ,
where we can optimize the values of α > 0 and β > 0, since this equation has non random coefficients. For example, taking for simplicity

α = 1 8B 3 + 4η 3 , β = nα 2 ,
we obtain

x - 2B 4 + η 4 4B 3 + 2η 3 x 2 = 16η 0 (2B 3 + η 3 ) n + 8B 1 + 4η 1 n + 32(2B 3 + η 3 )(2B 2 + η 2 ) n 2 + 8η 5 (2B 3 + η 3 ) n .
1.1. Proof of Theorem 2.1. Let us now deduce Theorem 2.1 from Theorem 1.5. Let us first remark that with probability at least 1 -ε/2

D ≤ D + B 2 εn ,
because the variance of D is less than

B 2 2n
. For a given ε > 0, let us take

η 0 = B 2 εn , η 1 = B 1 , η 2 = B 2 , η 3 = B 3 and η 4 = B 4 . We get that R λ ( θ) -inf θ∈Θ R λ (θ) is smaller than the smallest positive non degenerate root of x - B 4 2B 3 x 2 = 48B 3 n B 2 nε + 12B 1 n + 288B 2 B 3 n 2 + 24 log(3/ε)B 3 n ,
with probability at least

1 - 5 ε 6 - 4 k=1 P B k > B k + η k .
According to the weak law of large numbers, there is n ε such that for any n ≥ n ε ,

4 k=1 P B k > B k + η k ≤ ε/6.
Thus, increasing n ε and the constants to absorb the second order terms, we see that for some n ε and any n ≥ n ε , with probability at least 1 -ε, the excess risk is less than the smallest positive root of

x - B 4 2B 3 x 2 = 13B 1 n + 24 log(3/ε)B 3 n .
Now, as soon as ac < 1/4, the smallest positive root of x -

ax 2 = c is 2c 1 + √ 1 -4ac
. This means that for n large enough, with probability at least

1 -ε, R λ ( θ) -inf θ R λ (θ) ≤ 15B 1 n + 25 log(3/ε)B 3 n ,
which is precisely the statement of Theorem 2.1, up to some change of notation.

1.2. Proof of Theorem 2.2. Let us now weaken Theorem 1.4 in order to make a more explicit non asymptotic result and obtain Theorem 2.2. From now on, we will assume that λ = 0. We start by giving bounds on the quantity defined in Theorem 1.3 in terms of

B = sup f ∈span{ϕ 1 ,...,ϕ d }-{0} f 2 ∞ /E[f (X)] 2 .
Since we have

X 2 = Q -1/2 λ X 2 ≤ dB, we get d = 1 n n i=1 X i 2 ≤ dB, B 1 = 2E X 2 θ 0 , X -Y 2 ≤ 2dB R(f * ), B 1 = 2 n n i=1 X i 2 θ 0 , X i -Y i 2 ≤ 2dB r(f * ), B 2 = 2E X 4 ≤ 2d 2 B 2 , B 2 = 2 n n i=1 X i 4 ≤ 2d 2 B 2 , B 3 = 40 sup E u, X 2 θ 0 , X -Y 2 : u ∈ R d , u = 1 ≤ 40B R(f * ), B 3 = sup 40 n n i=1 u, X i 2 θ 0 , X i -Y i 2 : u ∈ R d , u = 1 ≤ 40B r(f * ), B 4 = 10 sup E u, X 4 : u ∈ R d , u = 1 ≤ 10B 2 , B 4 = sup 10 n n i=1 u, X i 4 : u ∈ R d , u = 1 ≤ 10B 2 .
Let us put

a 0 = 2dB + 4dBα[R(f * ) + r(f * )] + η αn + 16B 2 d 2 αn 2 , a 1 = 3/4 -40αB[R(f * ) + r(f * )],
and

a 2 = 20αB 2 .
Theorem 1.4 applied with β = nα/2 implies that with probability at least 1 -η the excess risk R( f (erm) ) -R(f * ) is upper bounded by the smallest positive root of a 1 x -a 2 x 2 = a 0 as soon as a 2 1 > 4a 0 a 2 . In particular, setting ε = exp(-η) when (1.17) holds, we have

R( f (erm) ) -R(f * ) ≤ 2a 0 a 1 + a 2 1 -4a 0 a 2 ≤ 2a 0 a 1 .
We conclude that Theorem 1.6. For any α > 0 and ε > 0, with probability at least 1 -ε, if the inequality

(1.17) 80 (2 + 4α[R(f * ) + r(f * )])Bd + log(ε -1 ) n + 4Bd n 2 < 3 4B -40α[R(f * ) + r(f * )] 2
holds, then we have

(1.18) R( f (erm) )-R(f * ) ≤ J (2 + 4α[R(f * ) + r(f * )])Bd + log(ε -1 ) n + 4Bd n 2 , where J = 8/(3α -160α 2 B[R(f * ) + r(f * )])
Now, the Bienaymé-Chebyshev inequality implies

P r(f * ) -R(f * ) ≥ t ≤ E r(f * ) -R(f * ) 2 t 2 ≤ E[Y -f * (X)] 4 /nt 2 .
Under the finite moment assumption of Theorem 2.2, we obtain that for any ε ≥ 1/n, with probability at least 1 -ε,

r(f * ) < R(f * ) + E[Y -f * (X)] 4 .
From Theorem 1.6 and a union bound, by taking

α = 80B[2R(f * ) + E[Y -f * (X)] 4 -1
, we get that with probability 1 -2ε,

R( f (erm) ) -R(f * ) ≤ J 1 B 3Bd + log(ε -1 ) n + 4Bd n 2 , with J 1 = 640 2R(f * ) + E [Y -f * (X)] 4
. This concludes the proof of Theorem 2.2.

Remark 1.1. Let us indicate now how to handle the case when Q is degenerate. Let us consider the linear subspace S of R d spanned by the eigenvectors of Q corresponding to positive eigenvalues. Then almost surely Span{X i , i = 1, . . . , n} ⊂ S. Indeed for any θ in the kernel of Q, E θ, X 2 = 0 implies that θ, X = 0 almost surely, and considering a basis of the kernel, we see that X ∈ S almost surely, S being orthogonal to the kernel of Q. Thus we can restrict the problem to S, as soon as we choose

θ ∈ span X 1 , . . . , X n ∩ arg min θ n i=1 θ, X i -Y i 2 ,
or equivalently with the notation

X = (ϕ j (X i )) 1≤i≤n,1≤j≤d and Y = [Y j ] n j=1 , θ ∈ im X T ∩ arg min θ X θ -Y 2
This proves that the results of this section apply to this special choice of the empirical least squares estimator. Since we have R d = ker X ⊕im X T , this choice is unique. Finally, we also have that inequality (2.3) of the paper still holds by replacing d by rank(Q).

2. Proof of Theorem 3.1. We use the same notations as in Section 1. We write X for ϕ(X), therefore, the function f θ maps an input x to θ, x . We consider the change of coordinates

X = Q -1/2 λ X.
Thus, from (1.13), we have E X 2 = D. We will use

R(θ) = E ( θ, X -Y ) 2 , so that R(Q 1/2 θ) = E ( θ, X -Y ) 2 = R(f θ ). Let Θ = Q 1/2 λ θ; θ ∈ Θ ,
and consider

θ 0 = arg min θ∈Θ R(θ) + λ Q -1/2 λ θ 2 .
With these notations,

θ = Q -1/2 λ θ 0 , σ = E θ 0 , X -Y 2 , χ = sup u∈R d E u, X 4 1/2 E u, X 2 , κ = E X 4 1/2 E X 2 = E X 4 1/2 D , κ ′ = E θ 0 , X -Y 4 1/2 σ 2 ,
and

T = Θ = max θ,θ ′ ∈Θ θ -θ ′ .
For α > 0, we introduce

J i (θ) = θ, X i -Y i , J(θ) = θ, X -Y L i (θ) = α θ, X i -Y i 2 , L(θ) = α θ, X -Y 2 W i (θ) = L i (θ) -L i (θ 0 ), W (θ) = L(θ) -L(θ 0 ), and
r ′ (θ, θ ′ ) = λ Q -1/2 λ θ 2 -Q -1/2 λ θ ′ 2 + 1 nα n i=1 ψ L(θ) -L(θ ′ ) . Let θ = Q 1/2 λ θ ∈ Θ. We have (2.1) -r ′ (θ 0 , θ) = r ′ ( θ, θ 0 ) ≤ max θ 1 ∈Θ r ′ ( θ, θ 1 ) ≤ γ + max θ 1 ∈Θ r ′ (θ 0 , θ 1 ),
where the quantity γ = max

θ 1 ∈Θ r ′ ( θ, θ 1 ) -inf θ∈Θ max θ 1 ∈Θ
r ′ (θ, θ 1 ) can be made arbitrary small by a proper choice of the estimator. Using an upper bound r ′ (θ 0 , θ 1 ) that holds uniformly in θ 1 , we will control both left and right hand sides of (2.1).

To achieve this, we will upper bound

(2.2) r ′ (θ 0 , θ 1 ) = λ Q -1/2 λ θ 0 2 -Q -1/2 λ θ 1 2 + 1 nα n i=1 ψ -W i (θ 1 )
by the expectation of a distribution depending on θ 1 of a quantity that does not depend on θ 1 , and then use the PAC-Bayesian argument to control this expectation uniformly in θ 1 . The distribution depending on θ 1 should therefore be taken such that for any θ 1 ∈ Θ, its Kullback-Leibler divergence with respect to some fixed distribution is small (at least when θ 1 is close to θ 0 ). Let us start with the following result.

Lemma 2.1. Let f, g : R → R be two Lebesgue measurable functions such that f (x) ≤ g(x), x ∈ R. Let us assume that there exists h ∈ R such that x → g(x) + hx 2 /2 is convex. Then for any probability distribution µ on the real line,

f xµ(dx) ≤ g(x)µ(dx) + min sup f -inf f, h 2
Var(µ) .

Proof. Let us put x 0 = xµ(dx) The function

x → g(x) + h 2 (x -x 0 ) 2
is convex. Thus, by Jensen's inequality

f (x 0 ) ≤ g(x 0 ) ≤ µ(dx) g(x) + h 2 (x -x 0 ) 2 = g(x)µ(dx) + h 2 Var(µ).
On the other hand

f (x 0 ) ≤ sup f ≤ sup f + g(x) -inf f µ(dx) = g(x)µ(dx) + sup f -inf f.
The lemma is a combination of these two inequalities.

The above lemma will be used with f = g = ψ, where ψ is the increasing influence function

ψ(x) =            -log(2), x ≤ -1, log(1 + x + x 2 /2), -1 ≤ x ≤ 0, -log(1 -x + x 2 /2), 0 ≤ x ≤ 1, log(2),
x ≥ 1.

Since we have for any x ∈ R

-log 1 -x + x 2 2 = log 1 + x + x 2 2 1 + x 4 4 < log 1 + x + x 2 2 ,
the function ψ satisfies for any

x ∈ R * -log 1 -x + x 2 2 < ψ(x) < log 1 + x + x 2 2 .
Moreover

ψ ′ (x) = 1 -x 1 -x + x 2 2 , ψ ′′ (x) = x(x -2) 2 1 -x + x 2 2 2 ≥ -2, 0 ≤ x ≤ 1,
showing (by symmetry) that the function x → ψ(x) + 2x2 is convex on the real line.

For any θ ′ ∈ R d and β > 0, we consider the Gaussian distribution with mean θ ′ and covariance β -1 I:

ρ θ ′ (dθ) = β 2π d/2 exp - β 2 θ -θ ′ 2 dθ.
From Lemmas 1.2 and 2.1 (with µ the distribution of -W i (θ) + α X i 2 β when θ is drawn from ρ θ 1 and for a fixed pair (X i , Y i )), we can see that

ψ -W i (θ 1 ) = ψ ρ θ 1 (dθ) -W i (θ) + α X i 2 β ≤ ρ θ 1 (dθ)ψ -W i (θ) + α X i 2 β + min log(4), Var ρ θ 1 L i (θ) . Let us compute 1 α 2 Var ρ θ 1 L i (θ) = Var ρ θ 1 J 2 i (θ) -J 2 i (θ 1 ) = ρ θ 1 (dθ) J 2 i (θ) -J 2 i (θ 1 ) 2 - X i 4 β 2 = ρ θ 1 (dθ) θ -θ 1 , X i 2 + 2 θ -θ 1 , X i J i (θ 1 ) 2 - X i 4 β 2 = 2 X i 4 β 2 + 4L i (θ 1 ) X i 2 αβ . (2.3)
Let ξ ∈ (0, 1), and let us remark that

L i (θ 1 ) ≤ L i (θ) ξ + α θ -θ 1 , X i 2 1 -ξ .
We get min log(4),

Var ρ θ 1 L i (θ) = min log(4), 4α X i 2 L i (θ 1 ) β + 2α 2 X i 4 β 2 ≤ ρ θ 1 (dθ) min log(4), 4α X i 2 L i (θ) βξ + 2α 2 X i 4 β 2 + 4α 2 X i 2 θ -θ 1 , X i 2 β(1 -ξ) ≤ ρ θ 1 (dθ) min log(4), 4α X i 2 L i (θ) βξ + 2α 2 X i 4 β 2 + min log(4), 4α 2 X i 4 β 2 (1 -ξ) .
Let us now put a = 3 log(4) < 2.17, b = a + a 2 log(4) < 8.7 and let us remark that min log(4), x + min log(4), y

≤ log 1 + a min{log(4), x} + log(1 + ay)

≤ log 1 + ax + by , x, y ∈ R + . Thus min log(4), Var ρ θ 1 L i (θ) ≤ ρ θ 1 (dθ) log 1 + 4aα X i 2 L i (θ) βξ + 2α 2 X i 4 β 2 a + 2b 1 -ξ .
We can then remark that

ψ(x) + log(1 + y) = log exp[ψ(x)] + y exp[ψ(x)] ≤ log exp[ψ(x)] + 2y ≤ log 1 + x + x 2 2 + 2y , x ∈ R, y ∈ R + . Thus, putting c 0 = a + 2b 1 -ξ , we get (2.4) ψ -W i (θ 1 ) ≤ ρ θ 1 (dθ) log[A i (θ)],
with

A i (θ) = 1 -W i (θ) + α X i 2 β + 1 2 -W i (θ) + α X i 2 β 2 + 8aα X i 2 L i (θ) βξ + 4c 0 α 2 X i 4 β 2 .
Similarly, we define A(θ) by replacing (X i , Y i ) by (X, Y ). Since we have

E exp n i=1 log[A i (θ)] -n log[EA(θ)] = 1,
from the usual PAC-Bayesian argument, we have with probability at least 1 -ε, for any

θ 1 ∈ R d , ρ θ 1 (dθ) n i=1 log[A i (θ)] -n ρ θ 1 (dθ) log[A(θ)] ≤ K(ρ θ 1 , ρ θ 0 ) + log(ε -1 ) ≤ β θ 1 -θ 0 2 2 + log(ε -1 ).
From (2.2) and (2.4), with probability at least 1 -ε, for any

θ 1 ∈ R d , we get r ′ (θ 0 , θ 1 ) ≤ 1 α log 1 + E ρ θ 1 (dθ) -W (θ) + α X 2 β + 1 2 -W (θ) + α X 2 β 2 + 8aα X 2 L(θ) βξ + 4c 0 α 2 X 4 β 2 + β θ 1 -θ 0 2 2nα + log(ε -1 ) nα + λ Q -1/2 λ θ 0 2 -Q -1/2 λ θ 1 2 .
Moreover from (2.3) and α X 2 β = -L(θ 1 ) + ρ θ 1 (dθ)L(θ), we deduce that

ρ θ 1 (dθ) -W (θ) + α X 2 β 2 = Var ρ θ 1 L(θ) + W (θ 1 ) 2 = W (θ 1 ) 2 + 4αL(θ 1 ) X 2 β + 2α 2 X 4 β 2 . Proposition 2.2. With probability at least 1 -ε, for any θ 1 ∈ R d , r ′ (θ 0 , θ 1 ) ≤ 1 α log 1 + E -W (θ 1 ) + W (θ 1 ) 2 2 + 2 + 8a/ξ α X 2 L(θ 1 ) β + 1 + 8a/ξ + 4c 0 α 2 X 4 β 2 + β θ 1 -θ 0 2 2nα + log(ε -1 ) nα + λ Q -1/2 λ θ 0 2 -Q -1/2 λ θ 1 2 ≤ E J(θ 0 ) 2 -J(θ 1 ) 2 + 1 2α W (θ 1 ) 2 + (2 + 8a/ξ) X 2 L(θ 1 ) β + (1 + 8a/ξ + 4c 0 )α X 4 β 2 + β θ 1 -θ 0 2 2nα + log(ε -1 ) nα + λ Q -1/2 λ θ 0 2 -Q -1/2 λ θ 1 2 .
Using the triangular inequality and Cauchy-Schwarz's inequality, we get 1 α 2 E W (θ 1 ) 2 = E θ 1 -θ 0 , X 2 + 2 θ 1 -θ 0 , X J(θ 0 )

2 ≤ E θ 1 -θ 0 , X 4 1/2 + 2E θ 1 -θ 0 , X 4 1/4 E J(θ 0 ) 4 1/4 2 (2.5)

≤ χ θ 1 -θ 0 2 E θ 1 -θ 0 θ 1 -θ 0 , X 2 + 2 θ 1 -θ 0 σ κ ′ χ E θ 1 -θ 0 θ 1 -θ 0 , X 2 2
≤ χq max q max + λ θ 1 -θ 0 2 θ 1 -θ 0 χq max q max + λ + 2σ √ κ ′ 2 , and 1 α E X 2 L(θ 1 ) = E X θ 1 -θ 0 , X + X J(θ 0 )

2 ≤ E X 4 1/2 E θ 1 -θ 0 , X 4 1/4 + E J(θ 0 ) 4 1/4 2 (2.6)

≤ κD θ 1 -θ 0 χq max q max + λ + 2σ √ κ ′ 2 .
Let us put

R(θ) = R(θ) + λ Q -1/2 λ θ 2 , c 1 = 4(2 + 8a/ξ), c 2 = 4(1 + 8a/ξ + 4c 0 ), δ = c 1 κκ ′ Dσ 2 n + 2χ log(ε -1 ) n + c 2 κ 2 D 2 n 2 2 √ κ ′ σ + Θ √ χ 2 1 - 4c 1 κχD n .
We have proved the following result.

Proposition 2.3. With probability at least 1 -ε, for any

θ 1 ∈ R d , r ′ (θ 0 , θ 1 ) ≤ R(θ 0 ) -R(θ 1 ) + α 2 χ θ 1 -θ 0 2 2 √ κ ′ σ + θ 1 -θ 0 √ χ 2 + c 1 α 4β κD √ κ ′ σ + θ 1 -θ 0 √ χ 2 + c 2 ακ 2 D 2 4β 2 + β θ 1 -θ 0 2 2nα + log(ε -1 ) nα .
Let us assume from now on that θ 1 ∈ Θ, our convex bounded parameter set. In this case, as seen in (1.14), we have θ 0 -θ 1 2 ≤ R(θ 1 ) -R(θ 0 ). We can also use the fact that

√ κ ′ σ + θ 1 -θ 0 √ χ 2 ≤ 2κ ′ σ 2 + 2χ θ 1 -θ 0 2 .
We deduce from these remarks that with probability at least 1 -ε,

r ′ (θ 0 , θ 1 ) ≤ -1+ αχ 2 2 √ κ ′ σ+ Θ √ χ 2 + β 2nα + c 1 ακDχ 2β R(θ 1 )-R(θ 0 ) + c 1 ακDκ ′ σ 2 2β + c 2 ακ 2 D 2 4β 2 + log(ε -1 ) nα .
Let us assume that n > 4c 1 κχD and let us choose

β = nα 2 ,

Proofs of Theorems

2.1 and 2.2

Proof of Theorem

3.1

Plugging this into (2.1), we get

Computing the numerical values of the constants when ξ = 0.8 gives c 1 < 95 and c 2 < 1511.


Anne Bretagnolle

Hélène Liliane Lizzi

Hélène Mathian

Adrien Van Hamme

Animating the cities…Dynamic exploration

Keywords:

Introduction

The project consists in a generic visual tool for mapping and analyzing the evolution of a system of cities on a long-term period (2 centuries), in different countries. The project started in 2008 in the framework of a national grant1 . The first step aimed to build a dynamic urban database over two centuries, in its conceptual form as well as in its fill in. The conception of the visual tool for exploring the database represents the last part of this project and has been developed in 2010 [START_REF] Van Hamme | Animer les villes. Exploration dynamique de bases de données urbaines harmonisées dans le temps long[END_REF]. Three main objectives were challenged: Designing a generic tool, that could be possibly used in the future for several countries under study (United States, France, South Africa, and possibly China and India); Allowing a dynamic exploration of the different databases, which raises a range of difficulties as they are spatio-temporal and multi-level (cities themselves and systems of cities); Giving an interactive application, that could be used by a large range of Internet users, from pupils to researchers or urban planners. The resulting interface combines time, location and attributes, according to the triad model [START_REF] Peuquet | It's about time; A conceptual framework for the representation of temporal dynamics in geographic information systems[END_REF]. Attributes are characteristics of cities, such as population, surfaces, density, socio-economic indicators and will include also inter-urban accessibility. The application integrates interactivity in order to propose an exploratory and animated cartography of urban dynamics, according to an evolutive conception of the urban system (Pumain 1997). The aim is to stimulate the desire of exploration of a variety of trajectories (Dykes &al. 2004[START_REF] Andrienko | Exploratory Analysis of Spatial and Temporal Data -A Systematic Approach[END_REF]. This paper focuses on the two first implementations, United States and France.

1/ Background and Objectives: getting to see evolving cities on the long term period

Managing spatio-temporal databases Urban databases are particularly complicated when approached on the long term period. Cities perimeters evolve frequently, especially since the first revolution of transportation, which is the depart point of our study. Usually, when official urban statistics integrate historical data, they define them for the whole period according to their final delineation. For our part, we assume that urban delineation is an active component of the city evolution that needs to be taken into account. If the delineations are too restrictive and don't take into account the emergence of suburbs, at the end of the 19 th century, and the raise of suburban rings, in the second half of the 20 th century, the results of the statistical computations of socio-economic indicators may suffer from bias (such as an under-estimation of larger cities).

The originality of our approach lays on the evolving definition of the city that underlies the geospatial database as well as on the way the system of cities is shown in a dynamic and multilevel perspective. Indeed, the project consists in constructing harmonized urban databases, i.e. using time definitions of cities. For the United States, we started from a theoretical approach of an "evolving city model", and we used the "one hour time-budget", dedicated in different world region or time period to travel every day to work place, in order to construct spatial delineations of the cities at each period [START_REF] Anne | Measuring urbanisation in the United States, from the first colonial trading posts to the Metropolitan Areas (1790 to 2000)[END_REF]. In this time budget, we fit the municipalities of the 19 th century, the urban agglomeration with their suburbs of the first half of the 20 th century, and the functional urban area with their suburban rings of the second half of the 20 th century (Figure 1). Historical transportation networks were used to estimate the average euclidian distance that can be covered in the one-hour time budget, at each census data point (from 1790 to 2000, every ten years). Inside the resulting radius, all the building blocks were aggregated (municipalities or counties, depending on the period) and this aggregation form the "city" in our database. Different types of indicators are under the scope of the project Harmonie-cités and will be represented as attributes of cities on the future website: population and surface, economic activities, social categories, and accessibility on transportation networks. Currently, the website has been designed for the first thematic focus (population, density, surface).

2/ Approach and methods: adopting a multi-level perspective

Our purpose is to use the power of new technologies to conceive a visual tool integrating our conception of how cities are evolving in their main components. Our level of observation is the one of the city, as a whole. Cities are conceived as elements of an urban system, which mean that their evolution is constrained by their position in the system, either global in the urban hierarchy, either local in their spatial regional or national context. The application has thus to integrate this multilevel contextualisation of the cities.

Different points of view for visualizing cities

For the conception of the tool, attention was given to allow a dynamic and interactive exploration of cities at different scale levels. This is linked to one of the most important properties of cities : they are not only a concentration of people, enterprises and services, but they also are part of systems of cities through the multiplicity of relationships that they draw with other cities, at regional scale, national scale, macro-regional scale, and for some of them today, at the world scale. As in cartography where any database is accurate for one specific cartographic scale, the accuracy of one urban database is associated to one specific geographic level. In our case, the spatio-temporal reference level is the one of system of cities. This is not easy to report in the sense that people are generally used to focus on cities and not system of cities. The links that underlay the interactions inside the system are difficult to catch : flows data are rare, and generally the representations are just based on the design of transportation or information networks (air plane, railroad, Internet connections…). An ambitious goal of the project is to get the level "system of cities" more accessible to the user, enabling him to discover some relations between (the list is not exhaustive) :

• Accessibility of cities at regional or national scales and their economic dynamic (see the railroad "hub" of Chicago at the 19 th century, or Paris at the same period) • Regional locations and economic innovation cycles (the Rust Belt in the United States, the Sun Belt in France or in the United States) • Social dynamics and city sizes (concentration of "creative classes" and white collars in large metropolises for instance). The visualizing tool should allow this simultaneous exploration of the multi-level dynamics, going from cities to urban system through different significant levels of contexts that make sense for the exploration of the city dynamics.

The city level that we consider here consists first in aggregation of elementary local administrative units used by National Statistical Boards to collect census data (the building blocks being generally areas, districts, "arrondissements", "quartiers", municipalities etc.). Of course, this doesn't give an overall view of the dynamics "inside" the city, but it gives a synthetic representation of the huge transformation of urban spatial extent, city by city, since two centuries. For example, when following the evolution of Paris's building blocks, we start with one commune in 1831 and finish with more than one thousand communes in 2006. Depending on the cities, the building blocks consist in municipalities (US 1790(US -2000)), counties or sub-counties (US 1940(US -2000)), communes (France, 1831(France, -2006)). In addition to this way of considering city evolution, another important point is to allow a simultaneous representation of the city with some attributes (Figure 2) and a view of the city as a part of a larger urban system. Then, the challenge is to represent dynamically the co-evolution of the delineation and the characteristics of the resident population and activities constituting the city, as a good proxy for social and economic specialisation of the city, which is seen here as a whole, an element of a higher system. Different contextual levels can be considered, for systems of cities. Indeed, one city is generally involved in different kind of systems (regional, national, macro-regional…), which are interlocked. Even the notion of national system of cities is questionable: if cities are known to have more relationship with other cities inside national frontiers than with cities abroad, these national relations don't involve all the cities but only some of them, either the largest, either the closest inside regional contexts, either the ones which belong to the same economic specialisations [START_REF] Bretagnolle | The organisation of urban systems[END_REF]. In order to respect this variety of approaches, we have decided not to impose a unique definition of systems of cities but to give the user the possibility of choosing different perspectives. The system level can be :

• Regional (for example cities belonging to some administrative regions, in France, or to federal states, in United States) • National (all the French or US cities at one date) • Defined by the experts of the projects (for example in US at the 19 th century, when strong interdependences between city only occur only inside regional part of the country, like the North-Eastern coast, the South-Eastern one, the Californian one, etc.) • Semantic (for example choosing cities larger than a selected threshold, or belonging to the top ten, or other) or defined by the user himself (Figure 3). Using a multilevel approach allows to combine not only the vision of the object but also its location and relative evolution inside a larger system of cities, coupled with its representation as an evolving aggregation of building blocks. This combination is provided by the use of attributes that are not computed only for one city but also for the urban system level, as illustrated for the example of Toulouse (Figure 4). In that case, one can observe two different stages. First, a relative decline (until the half of the 20 th century), mainly due to the low industry activities and the low accessibility on national transportation networks. Then, a relative growth, partly linked to the coming of aeronautic industry and to the local amenities regarding the sun-belt context.

3/ Results: a generic, dynamic and interactive tool

We wanted the visualizing tool to be an interface with the spatio-temporal database, reflecting the way we conceptualized the city and the forces underlying their evolution and their coevolution patterns, under logics operating either at the level of the whole urban system either at the level of a more local urban systems. We present in this part the different conceptual choices that have been done relating time and dynamic trends, not only from a cartographical animation point of view but also from a semiological one.

Conceptualizing a generic structure

We have put a strong emphasis on the construction of a generic structure, as the application is currently illustrated by two urban systems but some other projects are developed by researchers of the laboratory (South Africa, India, China…). Concerning the database formats, we had to face the diversity of the projects : some urban databases have been constructed several years ago, by a diversity of researchers, and the database models are not the same. Consequently, we have chosen to go from urban databases to homogeneous cartographic databases, which constitute the inputs of the application and feed directly the cartographic module. The input is in a very standard format. It gives for each city, at each date, its characteristics and its geometry (generally a shape, but sometimes only a point because of a lack of data, for example in US database until 1940). The objects are also associated to a generic type of city : definition evolves throughout the two centuries, depending on the perimeter contained in the one-hour time budget. For example, New York or Paris's trajectory varies from municipality to agglomeration of municipalities then to functional area.

Dynamic: at the heart of the interface

The interface combines in a dynamic way three windows which offer to the user a multidimensional view of the spatio-temporal database. The windows are dynamically joined, two by two, and synchronized with the time line (Figure 5). The whole period is represented by the range of the timeline, and the time unit is the year. The continuous representation of the period (time line) should give the impression of continuity, and events are marked by red dots on the timeline. By this way, the time for animation is proportional to the real time. Cursor moves on the timeline and representations change only when passing over red dots: We chose not to interpolate the geometries between two recordings in order to give the closest representation of the database. Another reason for this choice is that interpolation was not possible for the former periods as we don't have the shapes of the cities but only their centre location (US, before 1940). Tools for animation are the classical ones (play, stop, pause...) and a cursor is used for localising the current year. It is also possible for the user to act on the time and to select the dates to be visualized. In that case, each user leads the exploration of the period, date to date with its own perception of the order, and the time becomes the user's time.

The multi-level approach

The three windows allow to visualize a view of the selected contextual urban system, from local to global (Figure 5, upper window), a view of the spatial evolution of a selected city (low left window) and the quantitative evolution of a selected indicator for the city and its context (low right window). Consequently, the interface combines the dynamic of an urban system as well as the focus on the dynamic of a selected city, revealed through its spatial extension and its thematic dimensions, in relation with its national or regional context. As described above, some contexts are predefined (given by experts as significant regional contexts) but the user may also choose one specific using the "zoom". The combination of the three linked windows gives an original powerful lighting of the spatio-temporal database.

Figure 5 : the interface of the project (Orléans, France)

At the level of the city, representing the succession of the perimeters of the selected city for such a long period (more than 20 census dates in US, more than 30 in France) is challenging. The solution adopted (Figure 6, example of Orléans, in France) has the advantage of showing very clearly and quite simply the dramatic spatial extension of cities. Each date is characterized by a particular shade inside a variation of scale colour, and the oldness of an area of the city is given by the resulting colour of the overlaying shades. The central areas are in dark red at the end of the simulation, whereas the recent suburban rings appear in light rose. The legend summaries the way the figure can be understood, detailing 30 different colours with referred dates being just impossible to do. Interactivity: thinking stimulation and self-discovery The interface has not been conceived as an animated "atlas", a series of urban images to be looked at passively. On the contrary, the architecture of the application invites the user to a self-discovery, allowing him to process by try and error. One can have a static approach by choosing one date or a dynamic one by choosing animation. The upper and largest window is dedicated to a systemic view, and even if one may zoom on the map, this window suggests that urban dynamic has to be explored at a national or regional scale because cities interact and co-evolve. The user may also select a city and follow in the bottom left window its spatial extension at the building block level. A third window completes this dynamic exploration with some representation of statistical measures. Thus it is possible to connect the evolution of the selected city with its spatial extension and to follow its position (Figure 7, red curve) relatively to the urban system, either national or local (blue curve), depending on the spatial coverage selected in the upper window.

Figure 7: Links between functionalities of the interface

The multiplicity of links between sub-windows are an invitation for the user to explore himself the dynamics of cities and the strong interdependences that exist between the trajectory of a city and the one of a regional or national system. One can discover progressively the potentialities of the website. At a very rough level, meaning without pointing on a precise spatial context or a precise date, it is the evolution of the whole system of cities that appears. For example, the U.S. interface allows the visualisation of the progressive space filling by cities, during the pioneer front of the 19th century (Figure 8). It also enlightens the remarkable persistency of the urban hierarchy through time (Figure 9), old Eastern cities as Boston or New York remaining constantly in the top 10 of the cities, even after the booming of the Western then Southern cities. At a more detailed level, the user may choose some spatial and temporal contexts that enlighten some specificity of the selected cities. For example, considering only the rust-belt region in order to analyse the economic or demographic trajectory of a city like Detroit, or considering the years 1950-1980 in order to follow the trajectories of sun-belt cities during the TIC innovation cycle.

Different levels of difficulties

The public that could be interested in this product is large, from teachers and pupils at schools to researchers and urban planners. We have planned an easy downloading of the main attributes tables, maps and graphs, and of the information sheets that are associated to the indicators. Thee difficulty levels will be available, easy (for pupils for example), average (one-page format) and "to go further" (the papers published by the researchers who have conceived the databases and the statistical exploitations).

4/ Perspectives and conclusion

The aim of the project was to conceive a thinking stimulation tool and not a simple animated atlas of urban maps. Even if we are still in an experimental phase of the conception, the visualisation tool is now ready to be tested and improved taking into account suggestions of users (validation phase). The website has already been tested with teachers and pupils at the "Salon de la Géomatique" of the International Festival of Saint-Dié (October 2010), and other tests will be necessary in the future. The application concerning South Africa should also be conceived in the next months. The second phase of the project will be longer, as it consists in working on the introduction of new time-harmonized data, for social categories and transportations inter-urban indicators.

The databases are still in construction. We hope that the strong association with research fields in urban studies will feed in the next years the visualization website and allow frequent updates and maintenance of the tool.

Technical Appendix:

The architecture of the website is principally based on the JavaScript libraries OpenLayers and HighChart, and a few GML layers (one per year, for each country). The statistical measures are also stored in a MySQL database to be able to make asynchronous calls, especially using PHP (in AJAX).

Figure 1 :

 1 Figure 1: From the evolution of the concept of city to the database model

Figure 2 :

 2 Figure 2: A city trajectory with its attributes (Toulouse, France)

Figure 3 :

 3 Figure 3: Systems of cities defined by "experts" at the same period, according to the railroad interdependencies.Left: National context (France, 1870), Right: Regional context (Laker country, US, 1870)

Figure

  Figure 4: A relative view of a city trajectory (Toulouse, France)

Figure 6 :

 6 Figure 6: step by step animation of Toulouse spatial growth (1954-2006, selected dates)

Figure 8 :

 8 Figure 8: Step by step animation of the evolution of urban patterns in U.S. (selected dates, 1810-1910)

Figure 9 :

 9 Figure 9: Step by step animation of the urban populations in US (selected dates, 1930-2000)

Action Nationale de la Recherche Corpus et Outils de la Recherche en SHS(2008)(2009)(2010)(2011). The project is called « Harmonie-cités : construction de bases de données harmonisées sur les populations, les activités et les réseaux des villes » (http://www.parisgeo.cnrs.fr/index.php?page=programmeDescription&id=3&langue=en).