Ilaria Giulini

M Pierre

Alquier Examinateur

M Gérard

Biau Examinateur

M Stéphane

Boucheron Examinateur

M Olivier

Catoni Directeur De Thèse

M Nicolò

Cesa-Bianchi Rapporteur

Mallat Examinateur

M Pascal

Massart Rapporteur

M Alexandre

Tsybakov Examinateur

Présentée et soutenue le 24 septembre 2015 devant le jury composé de

Results on the estimation of the covariance and Gram matrices follow from random matrix theory and the accuracy of the estimation is usually evaluated in terms of the operator norm.

. It consists of two steps: first the dataset is embedded in a space in which clusters are more evident and then clusters are separated using a standard algorithm, e.g. the k-means algorithm. Also in this case the number of clusters is assumed to be known in advance. Among all the definitions of the Laplacian matrix, we choose

In chapter 2 we investigate the problem of estimating the Gram matrix via the classical empirical estimator and we provide non-asymptotic dimension-free bounds for the approximation error |θ Gθ -θ Ḡθ| using as a tool the robust estimator introduced before.

The results presented above allow us to characterize the stability of principal component analysis independently of the dimension d.   .

We also provide a similar bound in terms of the Frobenius norm by slightly changing the definition of the estimator Q + .

Chapter 4: Spectral Clustering

In this last chapter we present a new algorithm for spectral clustering. The idea is to couple spectral clustering with some preliminary change of representation in a reproducing kernel

Introduction

This thesis focuses on obtaining generalization bounds for random samples in reproducing kernel Hilbert spaces. The approach consists in first obtaining non-asymptotic dimensionfree bounds in finite-dimensional spaces using some PAC-Bayesian inequalities and then in generalizing the results to separable Hilbert spaces. We will investigate questions arising in the field of unsupervised classification and more specifically spectral clustering. We will also make some preliminary experiments on the possible relevance of some new spectral clustering algorithms in image analysis.

The goal of image analysis is to recognize and classify patterns that may have been transformed by some set of transformations, such as translations, rotations or scaling, and that more generally are affected by the conditions in which the images have been taken. The challenge is to find a representation in which those transformations of a same pattern are grouped together. However this objective is difficult to achieve even in the case of translations, that represent the simplest kind of transformations. Indeed, translated copies of a single pattern span a high dimensional vector space, so that it is hard to delineate clusters.

What we envision in this thesis is to use kernel methods to characterize clusters as a set of directions in some potentially infinite-dimensional Hilbert space.

Kernel-based methods are a class of statistical learning methods used to detect general types of relations between data. Detecting linear structures in the data has been the subject of much research in statistics and machine learning, and standard techniques, such as linear support vector machines (SVMs) and principal component analysis (PCA [START_REF] Jolliffe | Principal component analysis[END_REF]), are efficient and well understood. However real-world data analysis problems often have a non-linear structure. Kernel methods provide a workaround. Indeed, they provide a way to deal with non-linear structures in the data with the efficiency of linear algorithms and they have been successfully used for supervised learning (e.g. SVMs) and for unsupervised learning (e.g. kernel-PCA [START_REF] Schölkopf | Advances in kernel methods. chapter Kernel Principal Component Analysis[END_REF]).

An early survey on support vector machines for pattern recognition can be found in Burges [START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF]. For an introduction to support vector machines and kernel methods we refer to Cristianini and Shawe-Taylor [START_REF] Cristianini | An Introduction to Support Vector Machines and Other Kernel-based Learning Methods[END_REF] and to Schölkopf and Smola [START_REF] Schölkopf | Learning with kernels : support vector machines, regularization, optimization, and beyond[END_REF]. For a survey on kernel methods we refer to Hofmann, Schölkopf and Smola [START_REF] Hofmann | Kernel methods in machine learning[END_REF] and to Shawe-Taylor and Cristianini [START_REF] Shawe | Kernel Methods for Pattern Analysis[END_REF].

The basic idea of kernel methods is that any type of data can be implicitly embedded, via a (non-linear) function called a feature map, in a higher-dimensional space, the feature space, in which traditional linear methods become efficient. In other words, they simplify the geometry of the problem but increase the dimension of the space which may now be extremely large (or even infinite).

Kernel-based methods consist of two parts: firstly, the embedding in the feature space is computed and secondly, a learning algorithm is designed to detect the linear structure of the embedded data. The interesting point is that the embedding can be done without ex-plicit knowledge of the feature map but by simply computing the inner products between the images of the points in the feature space, using a kernel function. This approach is called the kernel trick.

In these kernel-based methods an important role is played by kernel random matrices which are matrices of the form k(X i , X j ) where k is a positive semi-definite kernel and X 1 , . . . , X n is a random sample. We will mainly investigate the statistical framework where X 1 , . . . , X n is a sample of independent and identically distributed (i.i.d.) vectors drawn according to an unknown probability distribution P. The goal is to estimate the integral operator L k f (x) = k(x, y)f (y) dP(y) from the matrix k(X i , X j ), under suitable assumptions. This integral operator is related to the Gram operator defined in the feature space, that is a generalization, to infinite dimension, of the Gram matrix. We now present the Gram matrix in finite dimension, since most of the computations made in this thesis are done in finite dimension first, and then generalized to infinite dimension. To be able to go from finite dimension to infinite dimension, we will establish dimension free inequalities.

Let X ∈ R d be a random (column) vector distributed according to P. The Gram matrix G of X is the symmetric positive semi-definite d × d matrix

G = E XX ,
where E is the expectation with respect to P.

The study of the spectral properties of the Gram matrix is of interest in the case of a non-centered criterion and it coincides, in the case of centered data (i.e. E[X] = 0), with the study of the covariance matrix Σ = E X -E(X) X -E(X) .

The empirical Gram matrix is defined as

Ḡ = 1 n n i=1 X i X i
and it is obtained by replacing the distribution P in the definition of G with the sample distribution 1 n n i=1 δ X i , also called the empirical measure. It is known, by the law of large numbers, that the empirical Gram matrix converges to the true one as the sample size grows to infinity. Instead of the asymptotic regime, we consider a non-asymptotic approach, that provides, for any fixed sample size n, non asymptotic deviation inequalities holding true with probability close to one. random matrices with independent rows. More precisely he observes that the empirical Gram matrix can be written as 1 n A A, where A is the matrix whose i-th row is the vector X i and hence, by construction, A is a matrix with independent rows. A different approach is proposed by Tropp [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF]. He observes that the empirical Gram matrix can be expressed as a sum of independent random matrices and proves that the Bernstein inequality extends to matrices, providing an exponential concentration inequality for the operator norm of a sum of independent random matrices.

The empirical Gram matrix becomes less efficient when the data have a long tail distribution, to improve on this, we construct in chapter 1 a more robust estimator and provide non-asymptotic dimension-free bounds of its estimation error. We then extend these bounds to any separable Hilbert space, taking advantage of the fact that they are independent of the dimension. In chapter 2 we use such an estimator to deduce some dimension-independent results for the classical empirical estimator.

Bounds on the deviations of the empirical Gram operator from the true Gram operator in separable Hilbert spaces can be found in Koltchinskii and Lounici [START_REF] Koltchinskii | Asymptotics and concentration bounds for spectral projectors of sample covariance[END_REF] in the case of Gaussian random vectors. Similarly to our results, these bounds are dimension-free and they are characterized in terms of the operator norm of the Gram operator and of its trace (more precisely in terms on the effective rank which is the ratio between the trace of the Gram operator and its operator norm).

Many learning algorithms rely on the spectral properties of the covariance matrix, for example principal component analysis (PCA).

Principal component analysis is a classical dimensionality reduction method that transforms the original coordinates into a new reduced set of variables, called principal components, that correspond to the directions where the variance is maximal. Since this set of directions lies in the space generated by the eigenvectors associated with the largest eigenvalues of the covariance matrix of the sample, the dimensionality reduction is achieved by projecting the dataset into the space spanned by these eigenvectors, that in the following we call largest eigenvectors.

Asymptotic results regarding the PCA projectors are provided in Biau and Mas [START_REF] Biau | PCA-kernel estimation[END_REF]. Results on PCA in Hilbert spaces can be found in Koltchinskii and Lounici [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF], [START_REF] Koltchinskii | Normal approximation and concentration of spectral projectors of sample covariance[END_REF]. The authors study the problem of estimating the spectral projectors of the covariance operator by their empirical counterparts in the case of Gaussian centered random vectors, based on the bounds obtained in [START_REF] Koltchinskii | Asymptotics and concentration bounds for spectral projectors of sample covariance[END_REF], and in the setting where both the sample size n and the trace of the covariance operator are large.

In [START_REF] Shawe-Taylor | On the eigenspectrum of the gram matrix and its relationship to the operator eigenspectrum[END_REF], [START_REF] Shawe-Taylor | On the eigenspectrum of the gram matrix and the generalisation error of kernel pca[END_REF], Shawe-Taylor, Williams, Cristianini and Kandola present concentration bounds for sums of eigenvalues of a kernel random matrix and use these results to provide a bound on the performance of (kernel-)PCA. These studies are continued in Zwald, Bousquet and Blanchard [START_REF] Zwald | Statistical properties of kernel principal component analysis[END_REF]. The main idea of these works is to view the kernel random matrix as the empirical version of an underlying integral operator. A first study on the relationship between the spectral properties of a kernel matrix and the corresponding integral operator is done in Koltchinski and Giné [START_REF] Koltchinskii | Random matrix approximation of spectra of integral operators[END_REF] for the case of a symmetric square integrable kernel k. They prove that the ordered spectrum of the kernel matrix

K ij = 1
n k(X i , X j ) converges to the ordered spectrum of the kernel integral operator L k . In [START_REF] Rosasco | On learning with integral operators[END_REF], Rosasco, Belkin and De Vito study the connection between the spectral properties of the empirical kernel matrix K ij and those of the corresponding integral operator L k by introducing two extension operators on the reproducing kernel Hilbert space defined by k, that have the same spectrum (and related eigenfunctions) as K and L k respectively. The introduction of these extension operators defined on the same Hilbert space overcomes the difficulty of dealing with objects (K and L k ) operating in different spaces.

Several methods have been proposed in the literature in order to provide a more robust version of PCA, e.g. [START_REF] Candès | Robust principal component analysis?[END_REF], [START_REF] Minsker | Geometric median and robust estimation in banach spaces[END_REF]. In [START_REF] Candès | Robust principal component analysis?[END_REF], Candès, Li, Ma and Wright show that it is possible to recover the principal components of a data matrix in the case where the observations are contained in a low-dimensional space but arbitrarily corrupted by additive noise. An alternative approach is suggested by Minsker in [START_REF] Minsker | Geometric median and robust estimation in banach spaces[END_REF] where a robust estimator of the covariance matrix, based on the geometric median, is used to provide non-asymptotic dimension-independent results concerning PCA.

Chapter 3 of this thesis deals with PCA. More precisely, we study how the results presented in the previous chapters lead to a robust PCA without assuming any geometric structure on the data.

The problem of knowing spectral properties of kernel random matrices also arises in the study of spectral clustering. In von Luxburg, Belkin and Bousquet [START_REF] Ulrike Von Luxburg | Consistency of spectral clustering[END_REF] methods similar to the ones described above are used to prove consistency of spectral clustering and to show the superiority of normalized spectral clustering over unnormalized spectral clustering from a statistical point of view [START_REF] Ulrike Von Luxburg | Consistency of spectral clustering[END_REF].

Clustering is the task of grouping a set of objects into classes, called clusters, in such a way that objects in the same group are more similar to each other than to those in other groups. Spectral clustering techniques use the spectrum of some data-dependent matrices to perform clustering. These matrices can be either the affinity (or similarity) matrix [START_REF] Donath | Lower bounds for the partitioning of graphs[END_REF] or the Laplacian matrix [START_REF] Fiedler | A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory[END_REF].

Given X 1 , . . . , X n a set of points to cluster, the affinity matrix measures the similarity, in terms of the distance, between each pair of points and a common definition is

A i,j = exp - X i -X j 2 2σ 2 (1)
where σ is a free scale parameter. The Laplacian matrix is obtained by rescaling the affinity matrix by its row sums.

and we view L as the empirical version of the integral operator with kernel K(x, y) = K(x, y) ( K(x, z)dP(z)) 1/2 ( K(y, z)dP(z)) 1/2 .

Remark that, in the case when the affinity matrix A has the form described in equation [START_REF] Andén | Multiscale scattering for audio classification[END_REF], the kernel K is the Gaussian kernel

K(x, y) = exp - x -y 2 2σ 2 . ( 2
)
Since connections between empirical operators and their continuous counterparts have been studied in many works, e.g. [START_REF] Rosasco | On learning with integral operators[END_REF], [START_REF] Ulrike Von Luxburg | Consistency of spectral clustering[END_REF], the study of spectral clustering can start from the study of the kernel K.

Our generalization bounds hold when spectral clustering is performed in the feature space defined by a reproducing kernel, opening the possibility of a kernel spectral clustering algorithm as described in chapter 4.

Before we give a short overview of the topics and results covered in this thesis, we explain the meaning of PAC-Bayesian.

The PAC-Bayesian approach was first introduced by McAllester [START_REF] Mcallester | Pac-bayesian model averaging[END_REF] and then formalized in Seeger [START_REF] Seeger | Pac-bayesian generalisation error bounds for gaussian process classification[END_REF] and Catoni [START_REF] Catoni | Statistical learning theory and stochastic optimization[END_REF]. It provides deviation bounds (Probably Approximately Correct bounds) for the generalization error, using Bayesian prior and posterior parameter distributions. Those bounds are frequentist, they do not require any Bayesian hypotheses on the unknown sample distribution. The role of prior and posterior parameter distributions here is purely technical, they serve as tools to obtain inequalities holding true uniformly with respect to the parameter, but are not involved in describing the statistical assumptions made about the sample distribution.

Chapter 1 and Chapter 2: Estimation of the Gram matrix

In order to estimate the Gram matrix, we consider the related problem of estimating the quadratic form

θ Gθ = E[ θ, X 2 ], θ ∈ R d ,
that computes the energy in the direction θ. Remark that to recover G it is sufficient to use the polarization identity G i,j = e i Ge j = 1 4 (e i + e j ) G(e i + e j ) -(e i -e j ) G(e i -e j ) ,

where {e i } d i=1 is the canonical basis of R d . As shown in Catoni [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF], if the distribution of θ, X 2 has a heavy tail for some values of θ, the quality of the approximation provided by the classical empirical estimator can be improved, using some M -estimator with a suitable influence function and a scale parameter depending on the sample size.

We present a robust estimator of the Gram matrix G and we provide non-asymptotic dimension-free bounds on the approximation error using some PAC-Bayesian inequalities related to Gaussian perturbations.

To reduce the influence of the tail of the distribution of θ, X 2 we truncate the empirical estimator via a continuous bounded function ψ (that is close to the identity in a neighborhood of zero) and we center it by a parameter λ > 0 so that we get

r λ (θ) = 1 n n i=1 ψ θ, X i 2 -λ . Let α(θ) = sup{α ∈ R + | r λ (αθ) ≤ 0}
and remark that if ψ is the identity, then the empirical estimator is equal to λ/ α(θ) 2 .

Therefore it is natural to consider as an estimator of θ Gθ a quantity related to λ/ α(θ) 2 , for a suitable value of λ.

We first estimate a confidence region for θ Gθ. Starting from the empirical criterion r λ we perturb the true parameter θ with a Gaussian perturbation centered in θ and with covariance matrix β -1 I, where β > 0 is a free parameter. The use of a PAC-Bayesian inequality leads to upper and lower bounds for r λ that hold, with high probability, uniformly in θ (i.e. for any θ ∈ R d ) and for any choice of β. Choosing the perturbation that optimizes the result, we construct a uniform confidence interval for θ Gθ that we denote by

B -(θ) := sup (λ,β) Φ - λ α(θ) 2 ≤ θ Gθ ≤ inf (λ,β) Φ -1 + λ α(θ) 2 =: B + (θ), (3)
where Φ -and Φ -1 + are non-decreasing functions depending on λ, β and on the kurtosis

κ = sup θ∈R d E( θ,X 2 )>0 E θ, X 4 E θ, X 2 2 .
We define as an estimator of G the minimal (in terms of the Frobenius norm) symmetric matrix Q whose corresponding quadratic form belongs to the optimal confidence interval (3) for any θ in a finite δ-net of S d (the unit sphere of R d ). To be sure that Q is nonnegative, we consider as an estimator its positive part Q + .

The quadratic estimator θ Q + θ is, by construction, more stable than the empirical estimator θ Ḡθ and it can be explicitly computed using a convex algorithm. In chapter 1 we provide a bound on the approximation error |θ Gθ -θ Q + θ| that holds, uniformly in θ, under weak moment assumptions and that does not depend explicitly on the dimension d of the ambient space. The result states that, with probability at least 1 -2 , for any θ ∈ S d , |θ Gθ -θ Q + θ| ≤ 2 max{θ Gθ, σ}µ max{θ Gθ, σ} + 7δ Tr(G 2 ) + σ with σ > 0 a suitable small threshold and µ(t)

2.032(κ -1) n 0.73 Tr(G) t + 4.35 + log( -1 ) + 98.5 κ Tr(G) nt [START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF] where Tr(G) = E X 2 denotes the trace of the Gram matrix and n ≤ 10 20 (this bound on n is only used to replace log(log(n)) with a more informative and more striking numerical constant). The threshold σ can be chosen going to zero as the sample size grows to infinity.

We observe that the quantity κ -1 that appears in the first factor of µ corresponds to the variance term, and more precisely to Var θ, X 2 /(θ Gθ) 2 (at least for the values of θ where κ is reached).

Remark that if we use a known upper bound instead of the exact value of κ, the result also holds for this upper bound. Remark also that κ ≤ 3 in the case when the distribution P is Gaussian. The bound presented above does not depend explicitly on the dimension d, that has been replaced by the entropy term Tr(G)/max{θ Gθ, σ}.

As a consequence, since the bound is dimension-independent, it can be generalized to any infinite-dimensional Hilbert space. More precisely, let H be a separable Hilbert space and let P be an unknown probability distribution on H. In the infinite-dimensional setting, the analogous to the Gram matrix is the Gram operator G : H → H defined by

Gθ = θ, v H v dP(v).
We consider an increasing sequence of finite-dimensional subspaces (H k ) k such that H = k H k and using a continuity argument we define an optimal confidence region for Gθ, θ H as in equation ( 3). Given X 1 , . . . , X n ∈ H an i.i.d. sample drawn according to P, we define

V k = span{Π k X 1 , . . . , Π k X n },
where Π k is the orthogonal projector on H k , and we consider the operator

Q = G k • Π V k (5)
where Π V k is the orthogonal projector on V k and G k is a linear operator on V k such that G k θ, θ H belongs to the optimal confidence region for any θ in a finite δ-net of S H ∩ V k .

Denoting by Q + the positive part of Q, we bound the approximation error

| Gθ, θ H -Q + θ, θ H |,
uniformly on the unit sphere of H, with the only additional assumption that the trace of the Gram operator Tr(G) is finite.

At the end of chapter 1 we generalize these results to estimate the expectation of a symmetric random matrix and we consider the problem of estimating the covariance matrix in the case where the expectation of X is unknown. We show that we do not need to estimate E[X] but that we can divide instead the data points into groups of a given size q and view the dataset as a sample of size n/q of q-tuples. We use the PAC-Bayesian approach to construct a robust estimator of the covariance matrix and to provide non-asymptotic dimension-free bounds on the approximation error. Finally we present an alternative (robust) estimator, related to the estimated matrix Q, which is easier to implement and we provide some empirical results that show the stability of such an estimator with respect to the empirical one.

Chapter 3: Principal Component Analysis

We recall that the study of PCA is linked with the study of the eigenvectors corresponding to the largest eigenvalues of the covariance matrix (or of the Gram matrix). In the following we consider the case of the Gram matrix.

Denoting by p 1 , . . . , p d an orthonormal basis of eigenvectors of G, we observe that the i-th eigenvalue of the Gram matrix is λ i = p i Gp i (and the same holds for the robust estimator Q + ). Starting from this remark, we are able to prove that each eigenvalue in the ordered spectrum of the Gram matrix is well approximated by the corresponding eigenvalue in the ordered spectrum of Q + , under weak moment assumptions.

From now on we assume that the eigenvectors are ranked according to the decreasing order of their eigenvalues, so that λ 1 is the largest eigenvalue and p 1 the corresponding eigenvector. A method to determine the number of relevant components (which correspond to the largest eigenvectors of G) is based on the difference in magnitude between successive eigenvalues. Let us assume that a significant gap in the spectrum of the Gram matrix is present at λ r -λ r+1 .

Denoting by Π r (respectively Π r ) the orthogonal projector on the r largest eigenvectors of G (respectively Q + ), we provide a bound on the approximation error Π r -Π r ∞ , in terms of the operator norm, which only depends on the trace of the Gram matrix, on its largest eigenvalue λ 1 and on the inverse of the size of the eigengap, but does not depend on the dimension d of the ambient space. In particular, the result relates the behavior of the estimator Π r to the size of the spectral gap λ r -λ r+1 and a good quality of approximation is determined by a large eigengap.

In order to avoid this kind of requirement, we propose to replace the projection on the r largest eigenvectors of the Gram matrix by a smooth cut-off of its eigenvalues. More precisely we aim at estimating f (G) by f (Q + ), where f is a Lipschitz function which is one on the largest eigenvalues and is zero on the smallest ones. In the case where there exists a (sufficiently large) gap in the spectrum of the Gram matrix, the usual projection on the first eigenvectors coincide with a cut-off with a Lipschitz constant exactly equal to the inverse of the size of the eigengap.

The worst case reformulation of our bound on the approximation error in terms of the operator norm depends, as before, on the trace of the Gram matrix and on its largest eigenvalue λ 1 and it replaces the size of the eigengap by the Lipschitz constant. The result states that, with probability at least 1 -2 , for any 1/L-Lipschitz function f ,

f (G) -f (Q + ) ∞ ≤ L -1 B(λ 1 ) + 6B(λ 1 ) Tr(G) + 4B(λ 1 ) 2
where

B(λ 1 ) 4λ 1   2.032(κ -1) n 0.73 Tr(G) λ 1 + 4.35 + log( -1 ) + 98.5 κ Tr(G) nλ 1
Hilbert space in order to bring down the representation of classes to a lower-dimensional space.

In order to better explain our approach we briefly describe one of the most successful spectral clustering algorithms [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF]. Given a set of points to cluster and the number c of classes, the algorithm introduced by Ng, Jordan, Weiss in [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] computes the c largest eigenvectors of the Laplacian matrix L and put them in columns to form a new matrix X.

After renormalization, each row of X is treated as a vector of R c and points are clustered according to this new representation using classical methods as k-means.

Let P be an unknown probability distribution on a compact subset of some separable Hilbert space (of possibly infinite dimension). Our approach relies on viewing the Laplacian matrix as the empirical version of the integral operator with kernel

K(x, y) = K(x, y) ( K(x, z)dP(z)) 1/2 ( K(y, z)dP(z)) 1/2
and on replacing the mere projection on the c largest eigenvectors of L (computed in [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF]) by a power of the kernel operator defined by K. This iteration, justified by the analysis of spectral clustering in terms of Markov chains, performs a smooth truncation of its eigenvalues that leads to a natural dimensionality reduction. Therefore it makes it possible to propose an algorithm that automatically estimates the number of classes (when it is not known in advance).

The algorithm can be described as a change of representation induced by a change of kernel followed by a (greedy) classification. Let us now describe in more detail the change of representation underlying our proposal for spectral clustering.

We use the Laplacian kernel K to build the new kernel

K2m (x, y) = K(y, z 1 ) K(z 1 , z 2 ) . . . K(z 2m-1 , x) dP ⊗(2m-1) (z 1 , . . . , z 2m-1 ), m > 0.
Denoting by H the reproducing kernel Hilbert space defined by K and by φ the corresponding feature map, the kernel K2m defines a new reproducing kernel Hilbert space with feature map G m-1/2 • φ, where G is the Gram operator on H, so that

K2m (x, y) = G m-1/2 φ(x), G m-1/2 φ(y) H .
In order to stabilize the representation induced by K2m we project this feature map on the unit sphere. This is equivalent to consider the normalized kernel

K m (x, y) = K2m (x, x) -1/2 K2m (x, y) K2m (y, y) -1/2 .
The change of kernel described above is ideal since it requires the knowledge of the law P.

In order to obtain a practical algorithm we construct, from a random sample of points in H, an estimated version of this ideal change of representation and we prove the convergence of this empirical algorithm by non-asymptotic bounds that are deduced from the bounds obtained for Gram operators in Hilbert spaces.

We first introduce a robust estimator K of K that can be written as

K(x, y) = φ(x), φ(y) H ,
for a suitable (observable) feature map φ, and then we replace it in the definition of K2m to obtain the new kernel

H2m (x, y) = G 2m-1 φ(x), φ(y) H
where G is the Gram operator defined by

Gv = v, φ(z) H φ(z) dP(z), v ∈ H.
Remark that the Gram operator G is still not observable since it depends on P. According to the results presented for the estimation of Gram operators in Hilbert spaces, we consider as an estimator of G the operator Q defined as the positive part of the operator introduced in equation ( 5), so that we get the estimated kernel

H 2m (x, y) = Q 2m-1 φ(x), φ(y) H .
The accuracy of this estimation is given in terms of two suitable quantities δ 1 , δ 2 > 0, that are linked with the function µ defined in equation ( 4) and depend on the sample size n as 1/ √ n. More precisely, with probability at least 1 -4 , for any x, y ∈ supp(P),

| H 2m (x, y) -K2m (x, y)| ≤ 2m (2δ 1 + δ 2 ) (1 + 2δ 1 + δ 2 ) 2m-2 (1 -δ 1 ) 4m + φ(x) H φ(y) H .
Coupling spectral clustering with a preliminary change of representation in a reproducing kernel Hilbert space can be seen as the analogous for unsupervised learning of the support vector machine (SVM) approach to supervised classification. Indeed, first, the kernel trick increases the dimension of the representation to simplify its geometry and then, spectral clustering decreases the size of the representation again. While SVMs compute a separating hyperplane that can be seen as a classification rule bearing on a representation of dimension one (the normal direction to the hyperplane), we show on some examples that the new representation induced by the kernel K m sends clusters to the neighborhood of an orthonormal set of c vectors (where c is the number of classes), that are therefore the vertices of a (regular) simplex, making subsequent classification a trivial task.

We apply this strategy to image analysis and we suggest with a small example that it is possible to learn transformation invariant representations from datasets that contain small successive transformations of the same pattern.

Chapter 1

The Gram Matrix

Our first goal is to estimate the Gram matrix from an i.i.d. sample. Based on some PAC-Bayesian inequalities we introduce a robust estimator leading to dimension free bounds on the estimation error. In particular, the dimension of the ambient space is replaced by an entropy term that depends on the trace of the Gram matrix.

Introduction

Let X ∈ R d be a random vector distributed according to the unknown probability measure

P ∈ M 1 + (R d ). The Gram matrix of X G = E XX = x x dP(x)
is a symmetric positive semi-definite d × d random matrix. The aim is to estimate the Gram matrix from a given sample of n independent and identically distributed (i.i.d.) vectors X 1 , . . . , X n ∈ R d drawn according to P. In relation to this problem, we consider the quadratic form

N (θ) = θ, x 2 dP(x), θ ∈ R d ,
which computes the energy in the direction θ, where θ, x denotes the standard inner product in R d . It can be seen as the quadratic form associated with the Gram matrix G. Indeed, according to the polarization identity

ξ Gθ = 1 4 [N (ξ + θ) -N (ξ -θ)] ,
we can recover the Gram matrix from the quadratic form N using the formula

G i,j = e i Ge j = 1 4 [N (e i + e j ) -N (e i -e j )] ,
where

{e i } d i=1 is the canonical basis of R d . The classical empirical estimator N (θ) = 1 n n i=1 θ, X i 2 , θ ∈ R d , (1.1)
is obtained by replacing, in the definition of N , the distribution P by the sample distribution 1 n n i=1 δ X i and it can be seen as the quadratic form associated with the empirical

Gram matrix Ḡ = 1 n n i=1 X i X i .
By the law of large numbers, the empirical Gram matrix converges to G almost surely as n goes to infinity. However, if θ, X 2 has a heavy tail distribution for at least some values of θ, the quality of approximation provided by the empirical mean estimator can be improved [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF].

Using a PAC-Bayesian approach we introduce in section 1.2 a new robust estimator of the Gram matrix and we provide non-asymptotic dimension-free bounds on the approximation error under weak moment assumptions. Since this result does not depend explicitly on the dimension d of the ambient space, it can be generalized to any infinite-dimensional Hilbert space, with the only assumption that the trace of the Gram matrix is finite, as shown in section 1.3.

In section 1.4 we generalize the results to estimate the expectation of a symmetric random matrix, while in section 1.5 we consider the problem of estimating the covariance matrix in the case of unknown expectation. Finally in section 1.6 we propose some empirical results that compare the performance of our robust estimator to that of the classical empirical one.

Estimate of the Gram matrix 1.2.1 Preliminaries

Our goal is to estimate, for any θ ∈ R d , the quadratic form

N (θ) = θ, x 2 dP(x)
from an i.i.d. sample X 1 , . . . , X n ∈ R d drawn according to the unknown probability distribution P ∈ M 1 + (R d ). In order to construct a robust estimator of N we consider a truncated version of the empirical estimator N . For any λ > 0 and for any fixed θ ∈ R d , we define

r λ (θ) = 1 n n i=1 ψ θ, X i 2 -λ ,
where λ is a centering parameter and where ψ : R → R is defined as

ψ(t) =        log(2) if t ≥ 1 -log 1 -t + t 2 2 if 0 ≤ t ≤ 1 -ψ(-t) if t ≤ 0. (1.2)
The function ψ is symmetric non-decreasing and bounded, it satisfies

-log 1 -t + t 2 2 ≤ ψ(t) ≤ log 1 + t + t 2 2 , t ∈ R,
and its role is to reduce the influence of the tail of the distribution of θ, X 2 . The introduction of the function ψ is similar to the use of an influence function in robust statistics, although here we do not use ψ to modify the definition of the mean and we use a scale parameter depending on the sample size. The scale parameter is somehow hidden in our formulation, in fact its role is played by the norm θ . We do not introduce an explicit scale parameter because it would have been redundant with the ability to change θ and because we will obtain results holding uniformly for any θ ∈ R d .

We consider the vector α(θ)θ, in the direction of θ, where the multiplicative factor α(θ) is defined by

α(θ) = sup {α ∈ R + | r λ (αθ) ≤ 0} . (1.3)
Since the function α → r λ (αθ) is continuous, r λ ( α(θ)θ) = 0 as soon as α(θ) < +∞.

Considering the fact that the influence function ψ is close to the identity in a neighborhood of zero, we observe that

0 = r λ ( α(θ)θ) = 1 n n i=1 ψ α(θ) 2 θ, X i 2 -λ 1 n n i=1 α(θ) 2 θ, X i 2 -λ = α(θ) 2 N (θ) -λ.
This suggests to consider as an estimator of the quadratic form N (θ) the quantity N (θ) = λ/ α(θ) 2 , for a suitable value λ > 0.

In section 1.2.2 we present a non-asymptotic dimension-free bound on the approximation error that holds under weak moment assumptions ( Proposition 1.19). However, since this estimator N is (unfortunately) no more a quadratic form, in section 1.2.3 we construct a quadratic estimator for N .

A PAC-Bayesian approach

We use a PAC-Bayesian approach based on Gaussian perturbations of the parameter θ to first construct a confidence region for N (θ) that is uniform in θ and then to define and study the robust estimator N .

Starting from the empirical criterion r λ , for any θ ∈ R d , we perturb the parameter θ with the Gaussian perturbation π θ ∼ N (θ, β -1 I) of mean θ and covariance matrix β -1 I. The parameter β > 0 is a free real parameter that can be seen as the inverse of the variance of the perturbation in each direction and will be determined later. The introduction of the perturbation π θ is a technical tool to analyze the average value of N in a neighborhood of θ.

We first introduce some technical lemmas that are useful to derive an upper bound for the empirical criterion r λ .

Lemma 1.1. We have

θ , x 2 dπ θ (θ ) = θ, x 2 + x 2 β .
Proof. Let W ∈ R d be a random variable distributed according to π θ . Since W is a Gaussian random variable with mean θ and covariance matrix β -1 I, for any x ∈ R d , the random variable W, x is a one-dimensional Gaussian random variable with mean θ, x and variance x (β -1 I)x = x 2 β . Consequently, we have

θ , x 2 dπ θ (θ ) = E W, x 2 = E W, x 2 + Var W, x = θ, x 2 + x 2 β ,
which concludes the proof.

As a consequence, we get

ψ θ, x 2 -λ = ψ θ , x 2 - x 2 β -λ dπ θ (θ ) .
What we would like to do is to pull the expectation with respect to π θ out of the influence function ψ, with a minimal loss of accuracy. To this purpose, we introduce the function

χ(z) =        ψ(z) z ≤ z 1 ψ(z 1 ) + p 1 (z -z 1 ) -(z -z 1 ) 2 /8 z 1 ≤ z ≤ z 1 + 4p 1 ψ(z 1 ) + 2p 2 1 z ≥ z 1 + 4p 1 (1.4)
where z 1 ∈ [0, 1] is such that ψ (z 1 ) = -1/4 and p 1 is defined by the condition p 1 = ψ (z 1 ). Using the explicit expression of the first and second derivative of ψ,

ψ (z) = 1 -z 1 -z + z 2 /2 and ψ (z) = -z + z 2 /2 (1 -z + z 2 /2) 2 ,
we can compute z 1 and p 1 . We find

z 1 = 1 -4 √ 2 -5, ψ(z 1 ) = -log 2( √ 2 -1) ,
which implies

p 1 = ψ (z 1 ) = 4 √ 2 -5 2( √ 2 -1)
,

and sup χ = ψ(z 1 ) + 2p 2 1 = 1 + 2 √ 2 2 -log 2( √ 2 -1) .
Further we observe that, for any z ∈ R,

ψ(z) ≤ χ(z).
Indeed, the inequality is trivial for

z ≤ z 1 , since χ(z) = ψ(z). For z ∈ [z 1 , z 1 + 4p 1 ],
performing a Taylor expansion at z 1 , we obtain that

ψ(z) = ψ(z 1 ) + p 1 (z -z 1 ) - 1 8 (z -z 1 ) 2 + z z 1 ψ (u) 2 (z -u) 2 du ≤ ψ(z 1 ) + p 1 (z -z 1 ) - 1 8 (z -z 1 ) 2 = χ(z), since ψ (u) ≤ 0 for u ∈ [0, 1[. Finally we observe that, for any z ≥ z 1 + 4p 1 , χ(z) = ψ(z 1 ) + 2p 2 1 > log(2) ≥ ψ(z).
On the other hand, for any z ∈ R,

χ(z) ≤ log 1 + z + z 2 /2 . (1.5)
Indeed, for z ≤ z 1 , we have already seen that the inequality is satisfied since χ(z) = ψ(z). Moreover we observe that the function

f (z) = log 1 + z + z 2 /2
is such that f (z 1 ) ≥ χ(z 1 ) and also f (z 1 ) ≥ χ (z 1 ). Performing a Taylor expansion at z 1 , we get

f (z) = f (z 1 ) + f (z 1 )(z -z 1 ) + z z 1 f (u)(z -u) 2 du ≥ χ(z 1 ) + χ (z 1 )(z -z 1 ) + inf f (z -z 1 ) 2 2 .
Since for any t

∈ [z 1 , z 1 + 4p 1 ], inf f = f ( √ 3 -1) = -1/4 = χ (t), we deduce that f (z) ≥ χ(z 1 ) + p 1 (z -z 1 ) - 1 8 (z -z 1 ) 2 = χ(z).
In particular, f (z 1 + 4p 1 ) ≥ χ(z 1 + 4p 1 ). Recalling that f is an increasing function while χ is constant on the interval [z 1 + 4p 1 , +∞[, we conclude that, for all z ∈ R, equation (1.5) is satisfied.

We are now going to prove that the quantity ψ θ, x 2 -λ can be bounded by the expectation of the function χ up to a loss term. Before, we introduce a lemma that allows us to pull the expectation out of the function χ.

Lemma 1.2. Let Θ be a measurable space. For any ρ ∈ M 1 + (Θ) and any h ∈ L 1 ρ (Θ),

χ h dρ ≤ χ(h) dρ + 1 8 Var h dρ ,
where by definition

Var h dρ = h(θ) -h dρ 2 dρ(θ) ∈ R ∪ {+∞}.
Proof. Performing a Taylor expansion of the function χ at z = h dρ we see that

χ h(θ) = χ(z) + h(θ) -z χ (z) + h(θ) z h(θ) -u χ (u) du ≥ χ(z) + h(θ) -z χ (z) + inf χ h(θ) -z 2 2 .
Remarking that inf χ = -1/4 and integrating the inequality with respect to ρ we obtain We conclude by taking the minimum of these two bounds.

χ h(θ) dρ(θ) ≥ χ h dρ - 1 8 h(θ) -h(θ)dρ 2 dρ(θ) = χ h dρ - 1
If we apply Lemma 1.3 to our problem we obtain that

ψ θ, x 2 -λ = ψ θ , x 2 - x 2 β -λ dπ θ (θ ) ≤ χ θ , x 2 - x 2 β -λ dπ θ (θ ) + min log(4), 1 8 Var θ , x 2 dπ θ (θ ) .
Defining m = θ, x , σ = x √ β and W ∼ N (0, σ 2 ) a centered Gaussian random variable, we observe that the variance of the random variable (m+W ) 2 is equal to Var θ , x 2 dπ θ (θ ) . Hence

Var θ , x 2 dπ θ (θ ) = Var (m + W ) 2 = Var W 2 + 2mW = E (W 2 + 2mW ) 2 -E W 2 + 2mW 2 = E W 4 + 4m 2 σ 2 -σ 4 = 4m 2 σ 2 + 2σ 4
where in the last line we have used the fact that E bigl(W 4 = 3σ 4 . We conclude that

Var θ , x 2 dπ θ (θ ) = 4 θ, x 2 x 2 β + 2 x 4 β 2 .
We have thus proved that

ψ θ, x 2 -λ ≤ χ θ , x 2 - x 2 β -λ dπ θ (θ ) + min log(4), θ, x 2 x 2 2β + x 4 4β 2 .
Lemma 1.4. For any positive real numbers a, b and c, and any Gaussian random variable

W ∼ N (0, σ 2 ), min a, bm 2 + c ≤ E min 2a, 2b m + W ) 2 + 2c .
Proof. We first claim that for any positive real numbers a, b, c,

min a, bm 2 + c ≤ min a, b(m + W ) 2 + c + min a, b(m -W ) 2 + c , since bm 2 + c ≤ max b(m + W ) 2 + c, b(m -W ) 2 + c .
Taking the expectation with respect to W of this inequality and remarking that W and -W have the same probability distribution we conclude that min a, bm

2 + c ≤ 2E min a, b(m + W ) 2 + c = E min 2a, 2b(m + W ) 2 + 2c .
Applying this result to our problem with a = log(4), b = x 2 /(2β), c = x 4 /(4β 2 ), m = θ, x , and m + W ∼ θ , x when θ ∼ π θ , we get

ψ θ, x 2 -λ ≤ χ θ , x 2 - x 2 β -λ dπ θ (θ ) + min 4 log(2), θ , x 2 x 2 β + x 4 2β 2 dπ θ (θ ).
Lemma 1.5. For any positive constants b, y and any z ∈ R,

χ(z) + min{b, y} ≤ log 1 + z + z 2 2 + y exp(sup χ) exp(b) -1 b .
Proof. For any positive real constants a, b, Applying this inequality to a = exp χ(z) and reminding that χ(z) ≤ log 1 + z + z 2 /2 , we conclude the proof.

According to Lemma 1.5, choosing b = 4 log(2), z = θ , x 2 -x 2 /β -λ and y = θ , x 2 x 2 /β + x 4 /2β 2 , we get

ψ θ, x 2 -λ ≤ log 1 + θ x 2 - x 2 β -λ + 1 2 θ , x 2 - x 2 β -λ 2 + c x 2 β θ , x 2 + x 2 2β dπ θ (θ ),
where

c = 15 4 log(2) exp (sup χ) = 15 8 log(2)( √ 2 -1) exp 1 + 2 √ 2 2 ≤ 44.3. (1.6)
In terms of the empirical criterion r λ we have proved the following result.

Proposition 1.6. With the above notation,

r λ (θ) ≤ 1 n n i=1 log 1 + θ , X i 2 - X i 2 β -λ + 1 2 θ , X i 2 - X i 2 β -λ 2 + c X i 2 β θ , X i 2 + X i 2 2β dπ θ (θ ).
We are now ready to use the following general purpose PAC-Bayesian inequality.

Proposition 1.7. Let ν ∈ M 1 + (R d ) be a prior probability distribution on R d and let f : R d × R d → [a, +∞] be a measurable function with a > -1. For any ∈ ]0, 1[, with probability at least 1 -, for any posterior distribution

ρ ∈ M 1 + (R d ), 1 n n i=1 f (X i , θ ) dρ(θ ) ≤ log E exp f (X, θ ) dρ(θ ) + K(ρ, ν) + log -1 n (1.7)
where

K(ρ, ν) =      log dρ dν dρ if ρ ν +∞ otherwise
is the Kullback divergence of ρ with respect to ν. By convention, a non measurable event is said to happen with probability at least 1 -when it includes a measurable event of probability non-smaller than 1 -.

Before proving the result we make some comments. Intuitively the above PAC-Bayesian inequality replaces the empirical mean with the expectation with respect to P. Doing this we loose in accuracy and this lost is evaluated by an entropy term which depends on the Kullbach divergence and on the confidence level set by the parameter .

In the following we consider as prior distribution ν = π 0 ∼ N (0, β -1 I). Since the result holds for any ρ ∈ M 1 + (R d ), if we consider as a family of posterior distributions the Gaussian perturbations π θ ∼ N (θ, β -1 I) | θ ∈ R d , β > 0 we get a bound that is uniform in θ and we can choose β to optimize the bound. To conclude we observe that with the above choice of prior and posteriors the Kullback divergence is

K(π θ , π 0 ) = β θ 2 2 .
Proof. We divide the proof into two parts.

Step 1. We prove that given ν ∈ M 1 + (R d ) a prior probability distribution and f : R d × R d → [a, +∞] a measurable function, there exists a measurable function F (X 1 , . . . , X n ) such that exp sup

ρ∈M 1 + (R d ) K(ρ,ν)<+∞ n 1 n n i=1 f (X i , θ ) -log E exp(f (X, θ )) × × 1 E exp f (X, θ ) < +∞ dρ(θ ) -K(ρ, ν) ≤ F (X 1 , . . . , X n ) and E F (X 1 , . . . , X n ) ≤ 1.
To give a meaning to this statement, we define for any measurable function h :

R d → R h(θ ) dρ(θ ) =              -∞ when min{h, 0} dρ = -∞ +∞ when min{h, 0} dρ > -∞ and max{h, 0} dρ = +∞ h dρ when h ∈ L 1 ρ (R d ).
We also introduce, for

θ ∈ R d , M (θ ) = 1 n n i=1 f (X i , θ ) and B(θ ) = log E exp(f (X, θ )) .
For any

ρ ∈ M 1 + (R d ) such that K(ρ, ν) < +∞, by definition of the Kullbach divergence, exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dρ(θ ) -K(ρ, ν) = exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ -log dρ dν (θ ) dρ(θ ) .
We put

h ρ = n(M -B)1[B < +∞] -log dρ dν .
In the case when

h ρ ∈ L 1 ρ (R d ), by the Jensen inequality, exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dρ(θ ) -K(ρ, ν) ≤ exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ -log dρ dν (θ ) dρ(θ ) = exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ -log dρ dν (θ ) dρ dν (θ ) 1 dρ dν (θ ) > 0 dν(θ ) = exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ 1 dρ dν (θ ) > 0 dν(θ ) ≤ exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dν(θ ).
Else, when min{h ρ , 0} dρ = -∞, then, with our conventions about ill-defined integrals, By using the definition of h ρ we conclude that

exp h ρ dρ = 0 ≤ exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dν(θ ) ∈ [0, +∞].
+∞ = exp h ρ dρ = exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ 1 dρ dν (θ ) > 0 dν(θ ) ≤ exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dν(θ ).
We have then proved that

exp h ρ dρ ≤ exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dν(θ ).
Hence, taking the supremum in ρ, we obtain that exp sup

ρ∈M 1 + (R d ) K(ρ,ν)<+∞ n M (θ ) -B(θ ) 1 B(θ ) < +∞ dρ(θ ) -K(ρ, ν) ≤ exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dν(θ ) def = F (X 1 , . . . , X n ).
The function F defined above is measurable with respect to (X 1 , . . . , X n ) and its expectation satisfies, from the Fubini theorem for non-negative functions,

E F (X 1 , . . . , X n ) = E exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ dν(θ ).
To conclude we need to prove that E exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ = 1. When B(θ ) = +∞, this is obvious. When B(θ ) < +∞, this means, by definition of B(θ ), that E exp f (X, θ ) < +∞. In this case,

E exp nM (θ ) = E exp n i=1 f (X i , θ ) = n i=1 E exp f (X i , θ ) = exp nB(θ ) , since X 1 , . . . , X n is an i.i.d. sample, so that E exp n M (θ ) -B(θ ) 1 B(θ ) < +∞ = 1.
Step 2. By the Markov inequality, since E F (X 1 , . . . , X n ) ≤ 1, we have

P F (X 1 , . . . , X n ) ≥ -1 ≤ .
Hence, with probability at least 1 -, log F (X 1 , . . . , X n ) ≤ log( -1 ).

The event log F (X 1 , . . . , X n ) ≤ log( -1 ) is a measurable event included in the, non necessarily measurable, event

A = sup ρ∈M 1 + (R d ), K(ρ,ν)<+∞ n M (θ ) -B(θ ) 1 B(θ ) < +∞ dρ(θ ) -K(ρ, ν) ≤ log( -1 ) .
Let us consider the event studied in the proposition, defined as

B = ∀ ρ ∈ M 1 + (R d ) s.t. K(ρ, ν) < +∞, M (θ ) dρ(θ ) ≤ B(θ ) dρ(θ ) + K(ρ, ν) + log( -1 ) n .
We observe that the inequality imposed on ρ in event B is trivially satisfied when B(θ ) dρ(θ ) = +∞. Thus, we have to check the inequality only when B(θ ) dρ(θ ) < +∞.

When A is satisfied, then for any

ρ ∈ M 1 + (R d ) such that K(ρ, ν) is finite, B(θ ) dρ(θ ) < +∞ and 1 B(θ ) = ∞ dρ(θ ) = 0, we have M (θ ) dρ(θ ) = M (θ )1 B(θ ) < +∞ dρ(θ ) ≤ B(θ )1 B(θ ) < +∞ dρ(θ ) + K(ρ, ν) + log( -1 ) n = B(θ ) dρ(θ ) + K(ρ, ν) + log( -1 ) n ,
proving that event B is also satisfied. In summary, we have proved that log F (X 1 , . . . , X n ) ≤ log( -1 ) ⊂ A ⊂ B.

Moreover, since the event log F (X 1 , . . . , X n ) ≤ log( -1 ) is measurable and is realized with probability at least 1 -, we conclude that, with probability at least 1 -, the event B is also realized, which proves the proposition.

We now come back to our problem and we apply the PAC-Bayesian inequality (1.7) to Proposition 1.6 in order to obtain a uniform bound on the empirical criterion. As already said, the prior is ν = π 0 ∼ N (0, β -1 I) and the family of posterior distributions is

π θ ∼ N (θ, β -1 I) | θ ∈ R d , β > 0 .
In what follows, the parameter can take any value in the interval ]0, 1[.

Proposition 1.8. With probability at least 1 -, for any θ ∈ R d , r λ (θ) ≤ E t(X, θ ) + 1 2 t(X, θ ) 2 + c X 2 β θ , X 2 + X 2 2β dπ θ (θ ) + β θ 2 2n + log( -1 ) n , where t(x, θ ) = θ , x 2 - x 2 β -λ and c = 15 8 log(2)( √ 2-1) exp 1+2 √ 2 2 .
Proof. We apply the PAC-Bayesian inequality (1.7) to Proposition 1.6 with

f (X i , θ ) = log 1 + t(X i , θ ) + 1 2 t(X i , θ ) 2 + c X i 2 β θ , X i + X i 2 2β . (1.8)
We obtain that, with probability at least 1 -, for any

θ ∈ R d , r λ (θ) ≤ 1 n n i=1 f (X i , θ )dπ θ (θ ) ≤ log E exp(f (X, θ )) dπ θ (θ ) + K(π θ , π 0 ) + log( -1 ) n ,
where

K(π θ , π 0 ) = β θ 2 2 . Observing that log E exp f (X, θ ) = log 1 + E t(X, θ ) + 1 2 t(X, θ ) 2 + c X 2 β θ , X 2 + X 2 2β ≤ E t(X, θ ) + 1 2 t(X, θ ) + c X 2 β θ , X 2 + X 2 2β ,
since log(1 + t) ≤ t, we conclude the proof.

In the following we introduce a sequence of results that provide more explicit bounds on the empirical criterion.

We first integrate explicitly with respect to the Gaussian distribution π θ .

Proposition 1.9. With probability at least 1 -, for any θ ∈ R d ,

r λ (θ) ≤ E θ, X 2 -λ + 1 2 θ, X 2 -λ 2 + 4 θ, X 2 X 2 β + 2 X 4 β 2 + c X 2 β θ, X 2 + 3 X 2 2β + β θ 2 2n + log( -1 ) n .
Proof. The proof is based on the identities

θ , X 2 dπ θ (θ ) = θ, X 2 + X 2 β , Var θ , X 2 dπ θ (θ ) = 4 θ, X 2 X 2 β + 2 X 4 β ,
that we already proved before.

Let us introduce

s 4 = E X 4 1/4 , (1.9) κ = sup θ∈R d E( θ,X 2 )>0 E θ, X 4 E θ, X 2 2 ,
(1.10) assuming that these two quantities are finite. Using the Cauchy-Schwarz inequality we rewrite the bound in Proposition 1.9 as follows.

Proposition 1.10. With probability at least 1 -, for any

θ ∈ R d , r λ (θ) ≤ κ 2 N (θ) -λ 2 + 1 + (κ -1)λ + (2 + c)κ 1/2 s 2 4 β N (θ) -λ + (κ -1)λ 2 2 + (2 + c)κ 1/2 s 2 4 λ β + (2 + 3c)s 4 4 2β 2 + β θ 2 2n + log( -1 ) n . (1.11)
Proof. The computation is based on the two inequalities

E X, θ 4 ≤ κN (θ) 2 , E θ, X 2 X 2 ≤ κ 1/2 s 2 4 N (θ),
and some elementary grouping of factors.

We now consider a finite set Λ ⊂ R + \ {0} 2 of possible values of the couple of parameters (λ, β) that will be determined later. Let |Λ| denote the number of elements of Λ.

Proposition 1.11. For any choice of (λ, β) ∈ Λ, we put

ξ = κλ 2 , µ = λ(κ -1) + (2 + c)κ 1/2 s 2 4 β , γ = λ 2 (κ -1) + (2 + c)κ 1/2 s 2 4 β + (2 + 3c)s 4 4 2β 2 λ + log(|Λ|/ ) nλ , δ = β 2nλ , (1.12)
where c =

15 8 log(2)( √ 2-1) exp 1+2 √ 2 2
, s 4 is defined by equation (1.9), and κ is defined by equation (1.10).

With probability at least 1 -, for any

θ ∈ R d , any (λ, β) ∈ Λ, r λ (θ) λ ≤ ξ N (θ) λ -1 2 + (1 + µ) N (θ) λ -1 + γ + δ θ 2 . (1.13)
Moreover, with probability at least 1 -, for any θ ∈ R d , any (λ, β) ∈ Λ,

r λ (θ) λ ≥ -ξ N (θ) λ -1 2 + (1 -µ) N (θ) λ -1 -γ -δ θ 2 . (1.14)
Proof. The first inequality is just a matter of rewriting the previous proposition with more compact notation and taking a union bound with respect to the allowed values of (λ, β) ∈ Λ. To obtain the second inequality, we observe that

-r λ (θ) = 1 n n i=1 ψ λ -θ, X i 2 .
Proceeding as previously done, with the necessary signs updates, we obtain that with probability at least 1 -, for any θ ∈ R d , any (λ, β) ∈ Λ,

-r λ (θ) ≤ E -t(X, θ ) + 1 2 t(X, θ ) 2 + c X 2 β θ , X 2 + X 2 2β dπ θ (θ ) + β θ 2 2n + log(|Λ|/ ) n .
Explicitly integrating with the respect to the Gaussian distribution π θ and rewriting the bound with compact notation, we conclude the proof.

From the above Proposition 1.11 we deduce a confidence interval for the quadratic form N (θ) in terms of λ/ α(θ) 2 , where α(θ) has been defined in equation (1.3) on page 17. More precisely the following proposition holds.

Proposition 1.12. With probability at least 1 -, for any θ ∈ R d , any (λ, β) ∈ Λ,

N (θ) ≤ Φ -1 θ,+ λ α(θ) 2 , (1.15)
where Φ θ,+ is the non-decreasing function defined as

Φ θ,+ (t) = t 1 + γ + δλ θ 2 /t 1 -µ -γ -2δλ θ 2 /t -1 1 ξ + µ + γ + 2δλ θ 2 /t < 1 .

Moreover, with probability at least 1 -, for any

θ ∈ R d , any (λ, β) ∈ Λ, Φ θ,- λ α(θ) 2 ≤ N (θ), (1.16)
where

Φ θ,-(t) = t 1 - γ + δλ θ 2 /t 1 + µ -γ -δλ θ 2 /t + 1 ξ -µ + γ + δλ θ 2 /t < 1 = t 1 - γ + δλ θ 2 /t 1 + µ -γ -δλ θ 2 /t 1 ξ -µ + 2γ + 2δλ θ 2 /t < 1
is an non-decreasing function.

Proof. We first observe that the polynomial

-ξ N (θ) λ -1 2 + (1 -µ) N (θ) λ -1 -γ -δ θ 2
can be written as

-ξ N (θ) λ -1 2 + 1 -µ - λδ θ 2 N (θ) N (θ) λ -1 -γ - λδ θ 2 N (θ) .
In order to simplify notation, we introduce the functions τ (θ) = λδ θ 2 N (θ) and

p θ (z) = -ξz 2 + 1 -µ -τ (θ) z -γ -τ (θ), z ∈ R.
We observe that τ (αθ) = τ (θ), and consequently p αθ (z) = p θ (z) for any α ∈ R + . We consider the case when p θ (1) > 0, meaning that

ξ + µ + γ + 2τ (θ) < 1.
(1.17)

In this case, the second degree polynomial p θ has two distinct real roots, z -1 and z +1 , where

z σ = 1 -µ -τ (θ) + σ 1 -µ -τ (θ) 2 -4ξ γ + τ (θ) 2ξ , σ ∈ {1, -1}.
We claim that

z -1 < γ + τ (θ) 1 -µ -γ -2τ (θ) < z +1 ,
which directly follows once we have shown that

p θ γ + τ (θ) 1 -µ -γ -2τ (θ) > 0. (1.18)
Observing that equation (1.17) can also be written as -ξ > -1 -µ -γ -2τ (θ) , we get

p θ γ + τ (θ) 1 -µ -γ -2τ (θ) = -ξ γ + τ (θ) 1 -µ -γ -2τ (θ) 2 + 1 -µ -τ (θ) γ + τ (θ) 1 -µ -γ -2τ (θ) -γ -τ (θ) > -γ + τ (θ) 2 1 -µ -γ -2τ (θ) + 1 -µ -τ (θ) γ + τ (θ) 1 -µ -γ -2τ (θ) -γ -τ (θ) = 0,
proving hence equation (1.18). We recall that by definition of α given in equation (1.3) on page 17, r λ (αθ) ≤ 0, for any α ≤ α(θ). Thus, according to equation (1.14), with probability at least 1 -, for any α ∈ [0, α(θ)],

p θ α 2 N (θ) λ -1 ≤ r λ (αθ) λ ≤ 0.
This means that

-1, α(θ) 2 N (θ) λ -1 ∩ z -1 , z +1 = ∅. Since p θ (1) > 0 > p θ (0), then z -1 ≥ 0 > -1. It follows that α(θ) 2 N (θ)/λ -1 ≤ z -1 .
This proves that, with probability at least 1 -, for any θ ∈ R d satisfying equation (1.17),

N (θ) ≤ λ α(θ) 2 1 + z -1 ≤ λ α(θ) 2 1 + γ + τ (θ) 1 -µ -γ -2τ (θ)
.

We observe that the last inequality

N (θ) ≤ λ α(θ) 2 1 + γ + τ (θ) 1 -µ -γ -2τ (θ)
can be written, when equation (1.17) is satisfied, as

Φ θ,+ N (θ) ≤ λ α(θ) 2 .
Moreover, this inequality is trivially true when condition (1.17) is not satisfied, because its left-hand side is equal to zero and its right-hand side is non-negative. Thus, we have proved that with probability at least 1 -, for any θ

∈ R d , Φ θ,+ N (θ) ≤ λ α(θ) 2 .
Proving the second part of the proposition requires a new argument and not a mere update of signs in the proof of the first part. Although it may seen at first sight that we are just aiming at a reverse inequality, the situation is more subtle than that. Let us first remark that in the case when

ξ -µ + γ + δ α(θ) 2 θ 2 < 1, (1.19)
is not satisfied, equation (1.16) is trivially satisfied because Φ θ,- λ α(θ) 2 = 0. Now, in the case when equation (1.19) is true, α(θ) < +∞, so that r λ α(θ)θ = 0, and therefore, according to equation (1.13), with probability at least 1 -,

0 ≤ q α(θ)θ α(θ) 2 N (θ) λ -1 , where q θ (z) = ξz 2 + (1 + µ)z + γ + δ θ 2 .
Since condition (1.19) can also be written as q α(θ)θ (-1) < 0, it implies that the second order polynomial q α(θ)θ has two roots and that α(θ) 2 N (θ) λ -1 is on the right of its largest root, which is larger than -1. On the other hand, we observe that, under condition (1.19), putting τ (θ) = δ α(θ) 2 θ 2 , we get

q α(θ)θ - γ + τ (θ) 1 + µ -γ -τ (θ) = ξ γ + τ (θ) 1 + µ -γ -τ (θ) 2 -(1+µ) γ + τ (θ) 1 + µ -γ -τ (θ) +γ+ τ (θ) < γ + τ (θ) 2 1 + µ -γ -τ (θ) - (1 + µ) γ + τ (θ) 1 + µ -γ + τ (θ) + γ + τ (θ) = 0.
Therefore, with probability at least 1

-, if condition (1.19) is satisfied, α(θ) 2 N (θ) λ -1 ≥ - γ + τ (θ) 1 + µ -γ -τ (θ) ,
which rewrites as equation (1.16) under condition (1.19). Hence, we conclude that with probability at least 1 -, for any θ ∈ R d , equation (1.16) holds.

To simplify notation, we omit the dependence on θ of the functions defined above, so that, from now on, we write Φ -and Φ + instead of Φ θ,-and Φ θ,+ .

Now that we have an observable non-asymptotic confidence interval for N (θ), we are going to define and study in the following an estimator for N (θ). Let Λ ⊂ R + \ {0} 2 be a finite set as before. We consider, for some energy level σ ∈ R + to be chosen later, the set

Γ = (λ, β, t) ∈ Λ × R + | ξ + µ + γ + 2 δλ max{t, σ} < 1
and the bound

B λ,β (t) =      γ + λδ/ max{t, σ} 1 -µ -γ -2λδ/ max{t, σ} (λ, β, t) ∈ Γ +∞ otherwise.
(1.20)

We consider the family of estimators

N λ (θ) = λ α(θ) 2 , θ ∈ R d , ( 1.21)
and we observe that, since α is homogeneous (of degree -1) in θ,

N λ (θ) = θ 2 N λ (θ/ θ ) .
Before defining our estimator and presenting some uniform bounds on the estimation error, we introduce a technical lemma.

Lemma 1.13. With probability at least 1 -, for any

θ ∈ R d , any (λ, β) ∈ Λ, max{N (θ), σ θ 2 } 1 + B λ,β θ -2 N (θ) -1 ≤ max{ N λ (θ), σ θ 2 }. (1.22) With probability at least 1 -, for any θ ∈ R d , any (λ, β) ∈ Λ, max{ N λ (θ), σ θ 2 } 1 -B λ,β θ -2 N λ (θ) ≤ max{N (θ), σ θ 2 }. With probability at least 1 -, for any θ ∈ R d , any (λ, β) ∈ Λ, max{N (θ), σ θ 2 } ≤ max{ N λ (θ), σ θ 2 } 1 + B λ,β θ -2 N λ (θ) .
Proof. We may assume, without loss of generality, that θ = 1. We observe that, for any

z, t, σ ∈ R + , if Φ + (z) ≤ t then Φ + (z) ≤ max{t, σ}.
Moreover, since by definition of Φ + , we have Φ + (σ) ≤ σ, in particular Φ + (σ) ≤ max{t, σ}.

We have then proved that, if Φ + (z) ≤ t, then Φ + max{z, σ} ≤ max{t, σ}.

The same reasoning applies to Φ

-proving that, if Φ -(z) ≤ t, then Φ -max{z, σ} ≤ max{t, σ}.
We also observe that, by definition of B λ,β , for any (λ, β) ∈ Λ,

Φ + max{z, σ} = max{z, σ} 1 + B λ,β (z) -1 (1.23) Φ -max{z, σ} ≥ max{z, σ} 1 -B λ,β (z) . (1.24)
Applying these results to Proposition 1.12 we obtain the first two inequalities.

To obtain the third inequality, we recall that, by definition of inverse function,

Φ -1 + (t) = sup{y | Φ + (y) ≤ t}. Denoting ȳ = sup{y | Φ + (y) ≤ t}, we get Φ -1 + max{t, σ} = max{ȳ, σ} ≤ max{t, σ} 1 + B λ,β (t) , (1.25)
since t ≤ ȳ and B λ,β is a non-increasing function. To conclude, it is sufficient to apply equation (1.25) to inequality

N (θ) ≤ Φ -1 + λ α(θ) 2 . Let us put ( λ, β) = arg min (λ,β)∈Λ B λ,β θ -2 N λ (θ) .
We define our robust estimator N as

N (θ) = N λ (θ) (1.26)
where N λ is defined in equation (1.21).

Next proposition provides a (uniform) bound on the approximation error.

Proposition 1.14. With probability at least 1 -2 , for any θ ∈ R d , max N (θ), σ θ 2 max N (θ), σ θ 2 -1 ≤ inf (λ,β)∈Λ B λ,β θ -2 N (θ) 1 + B λ,β θ -2 N (θ)
.

Proof. We may assume that θ = 1 without loss of generality, from homogeneity considerations. By Lemma 1.13 with the choice of parameters ( λ, β), we get that with probability at least 1 -2 ,

max{N (θ), σ} max{ N (θ), σ} -1 ≤ B λ, β N (θ) .
Since, by definition, ( λ, β) are the values which minimize

B λ,β N λ (θ) , B λ, β N (θ) = inf (λ,β)∈Λ B λ,β N λ (θ) = inf (λ,β)∈Λ B λ,β max{ N λ (θ), σ} . Equation (1.22) concludes the proof. Indeed, inf (λ,β)∈Λ B λ,β max{ N λ (θ), σ} ≤ inf (λ,β)∈Λ B λ,β max{N (θ), σ} (1 + B λ,β (N (θ))) -1 ≤ inf (λ,β)∈Λ B λ,β N (θ) (1 + B λ,β (N (θ))) -1 .
In the following we are going to present some simplified and more explicit bounds on the approximation error.

We introduce the subset Γ of Γ defined as

Γ = (λ, β, t) ∈ Λ × R + | ξ + µ + γ + 4δλ/ max{t, σ} < 1, µ + γ + 2δλ/ max{t, σ} ≤ 1/2,
and 2γ + δλ/ max{t, σ} ≤ 1/2 and the function

B λ,β (t) =      γ + λδ/ max{t, σ} 1 -µ -2γ -4λδ/ max{t, σ} (λ, β, t) ∈ Γ +∞ otherwise.
Lemma 1.15. For any (λ, β) ∈ Λ and any t ∈ R + ,

B λ,β t 1 + B λ,β (t) ≤ B λ,β (t), B λ,β (t) 1 -B λ,β (t) ≤ B λ,β (t).
Proof. We first observe that when (λ, β, t) ∈ Γ then B λ,β (t) = +∞ and hence the two inequalities are trivial. We now assume (λ, β, t) ∈ Γ and we put τ = λδ/ max{t, σ}. We prove the second inequality first. Since Γ ⊂ Γ, we have

B λ,β (t) 1 -B λ,β (t) = γ + τ 1 -µ -2γ -3τ ≤ B λ,β (t).
In order to prove the first inequality, we first check that λ, β, t 1+B λ,β (t) ∈ Γ. We first observe that, since

max {t/[1 + B λ,β (t)], σ} ≥ max{t, σ}/[1 + B λ,β (t)], then ξ + µ + γ + 2δλ/ max t 1 + B λ,β (t) , σ ≤ ξ + µ + γ + 2[1 + B λ,β (t)] δλ max{t, σ} = ξ + µ + γ + 2τ + 2τ B λ,β (t). Moreover, as (λ, β, t) ∈ Γ , we get B λ,β (t) = γ + τ 1 -µ -γ -2τ ≤ 1, so that ξ + µ + γ + 2δλ/ max{t/[1 + B λ,β (t)], σ} ≤ ξ + µ + γ + 4τ < 1, which proves that, indeed, λ, β, t 1+B λ,β (t) ∈ Γ. Therefore B λ,β t 1 + B λ,β (t) ≤ γ + τ (1 + B λ,β (t)) 1 -µ -γ -2τ [1 + τ B λ,β (t)] = (γ + τ ) (1 + τ /(1 -µ -γ -2τ )) 1 -µ -γ -2τ -2τ B λ,β (t) ,
where in the last line we have used the definition of B λ,β . Observing that

1 -µ -γ -2τ -2τ B λ,β (t) = (1 -µ -γ -2τ ) 2 -2τ (γ + τ ) 1 -µ -γ -2τ , we obtain B λ,β t 1 + B λ,β (t) ≤ (γ + τ ) (1 -µ -γ -τ ) (1 -µ -γ -2τ ) 2 -2τ (γ + τ ) = (γ + τ ) (1 -µ -γ -τ ) (1 -µ -γ -τ ) 2 + τ 2 -2τ (1 -µ -γ -τ ) -2τ 2 -2γτ = γ + τ 1 -µ -γ -τ -2τ -(τ 2 + 2γτ ) / (1 -µ -γ -τ ) .
Considering that

τ 2 + 2γτ / (1 -µ -γ -τ ) ≤ τ, since when (λ, β, t) ∈ Γ , it is true that 1 -µ -γ -τ ≥ 1/2 and 2γ + τ ≤ 1/2, we conclude that B λ,β t 1 + B λ,β (t) ≤ γ + τ 1 -µ -γ -4τ = B λ,β (t).
Combining the above lemma with Proposition 1.14, we obtain the following result.

Proposition 1.16. With probability at least

1 -2 , for any θ ∈ R d , max{N (θ), σ θ 2 } max{ N (θ), σ θ 2 } -1 ≤ inf (λ,β)∈Λ B λ,β θ -2 N (θ) .
We now compute an explicit upper bound for the right-hand side of the inequality presented in the above proposition.

We first observe that, for any

θ ∈ R d , θ -2 N (θ) = E θ -2 θ, X 2 ≤ E X 2 ≤ E X 4 1/2 = s 2 4 ,
where s 4 is defined in equation (1.9) on page 26. Let a > 0 and

K = 1 + a -1 log n 72(2 + c)κ 1/2 ,
where we recall that c =

15 8 log(2)( √ 2-1) exp 1+2 √ 2 2
and κ is defined in equation (1.10). The choice of K will be justified by equation (1.29). We fix a geometric grid for θ -2 N (θ)/s 2 4 of the form exp (-ja) , 0 ≤ j < K, and we define the finite set Λ of all possible values of the couple (λ, β) as

Λ = (λ j , β j ) | 0 ≤ j < K , ( 1.27)
where

λ j = 2 n(κ -1) (2 + 3c) 4(2 + c)κ 1/2 exp(-ja) + log(K/ ) β j = 2(2 + c)κ 1/2 s 4 4 n exp -(j -1/2)a .
Proposition 1.17. Let σ > 0 be a threshold such that σ ≤ s 2 4 . We define the explicit bound

ζ(t) = 2(κ -1) (2 + 3c)s 2 4 4(2 + c)κ 1/2 t + log(K/ ) cosh(a/4) + 2(2 + c)κ 1/2 s 2 4 t cosh(a/2)
and

B * (t) =      n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ}) 6 + (κ -1) -1 ζ(max{t, σ}) ≤ √ n +∞ otherwise.
For any t ∈ R + , inf

(λ,β)∈Λ B λ,β (t) ≤ B * min{t, s 2 4 } . (1.28)
As a consequence, with probability at least 1 -2 , for any

θ ∈ R d , max{N (θ), σ θ 2 } max{ N (θ), σ θ 2 } -1 ≤ B * θ -2 N (θ) .
Before proving the result we recall explicitly some definitions previously introduced. We have defined the set

Γ = (λ, β, t) ∈ Λ × R + | (4κ -3)λ 2 + 2(2 + c)κ 1/2 s 2 4 β + (2 + 3c)s 4 4 2β 2 λ + log(K/ ) nλ + β n max{t, σ} < 1 ,
and the bound

B λ,β (t) = (κ -1)λ 2 + (2 + c)κ 1/2 s 2 4 β + (2 + 3c)s 4 4 2β 2 λ + log(K/ ) nλ + β 2n max{t, σ} 1 - 3(κ -1)λ 2 - 2(2 + c)κ 1/2 s 2 4 β - (2 + 3c)s 4 4 2β 2 λ - log(K/ ) nλ - β n max{t, σ} when (λ, β, t) ∈ Γ, and B λ,β (t) = +∞ otherwise. The definition of N is given in equa- tion (1.26).
We remark that the function ζ is non-increasing and observable, so that in the case when the threshold σ is defined by the equation

6 + (1 -κ) -1 ζ(σ) = √ n,
then the condition in the definition of B * is satisfied for all t ∈ R + . Moreover the function B * : R + → R + is non-increasing and in this case

B * (t) ≤ B * (0) ≤ 1 2 , t ∈ R + .
We also observe that it is not restrictive to assume σ ≤ s 2 4 , since we know that N (θ) ≤ s 2 4 θ 2 . Moreover, we will see later that if we choose the threshold σ such that

8ζ(σ) = √ n,
then it decays to zero at speed 1/n as the sample size grows to infinity.

Proof. We recall that the function B λ,β is non-increasing so that B λ,β (t) ≤ B λ,β min{t, s 2 4 } . Moreover, since max min{t, s 2 4 }, σ = min max{t, σ}, s 2 4 , it is sufficient to prove equation (1.28) for max{t, σ} ∈ [0, s 2 4 ]. As equation (1.28) is trivial when B * (t) = +∞, we may assume that B * (t) < +∞, so that 6ζ(max{t, σ}) ≤ √ n. In particular, by considering only the second factor in the definition of ζ, we obtain that

2(2 + c)κ 1/2 s 2 4 max{t, σ} ≤ 2(2 + c)κ 1/2 s 2 4 max{t, σ} cosh(a/2) ≤ √ n 6 , which implies max{t, σ} s 2 4 ≥ 72(2 + c)κ 1/2 n ≥ exp (-a(K -1)) . Therefore we have log max{t, σ} s 2 4 ∈ -a(K -1), 0 (1.29)
which justifies the definition of K. Because of equation (1.29), there exists  ∈ {0, . . . , K -1}, for which log max{t, σ} s 2

4

+ a ≤ a/2.

Equivalently, it means that

exp - + 1 2 a ≤ max{t, σ} s 2 4 ≤ exp - - 1 2 a .
(1.30)

We now recall that by equation (1.12) on page 27

γ + δλ/ max{t, σ} = λ 2 (κ -1) + (2 + c)κ 1/2 s 2 4 β + (2 + 3c)s 4 4 2β 2 λ + log(K/ ) nλ + β 2n max{t, σ}
and we observe that (λ * , β * ) defined as

λ * = λ * (t) = 2 n(κ -1) (2 + 3c)s 2 4 4(2 + c)κ 1/2 max{t, σ} + log(K/ ) (1.31) β * = 2(2 + c)κ 1/2 s 2 4 max{t, σ}n (1.32)
are the desired values that optimize γ + δλ/ max{t, σ}. We also remark that, by equation (1.30),

β  exp(-a/2) ≤β * ≤ β  (1.33) λ  exp(-a/4) ≤λ * ≤ λ  exp(a/4). (1.34) Thus, evaluating γ + δλ/ max{t, σ} in (λ  , β  ) ∈ Λ, we obtain that γ  + δ  λ  / max{t, σ} = λ * (κ -1) 2 λ  λ * + (2 + c)κ 1/2 s 2 4 β * β * β  + (2 + 3c)s 4 4 2 β 2  λ * λ * λ  + log(K/ ) nλ * λ * λ  + β * 2n max{t, σ} β  β * ≤ λ * (κ -1) 2 λ  λ * + 1 nλ * (2 + 3c)s 2 4 4(2 + c)κ 1/2 max{t, σ} + log(K/ ) λ * λ  + (2 + c)κ 1/2 s 2 4 2n max{t, σ} β * β  + β  β * ≤ 2(κ -1) n (2 + 3c)s 2 4 4(2 + c)κ 1/2 max{t, σ} + log(K/ ) cosh log λ  λ * + 2(2 + c)κ 1/2 s 2 4 n max{t, σ} cosh log β  β * .
Using equation (1.33) and equation (1.34), we get

γ  + δ  λ  / max{t, σ} ≤ 2(κ -1) n (2 + 3c)s 2 4 4(2 + c)κ 1/2 max{t, σ} + log(K/ ) cosh a 4 + 2(2 + c)κ 1/2 s 2 4 n max{t, σ} cosh a 2 .
We also observe that

µ  + γ  + 4δ  λ  / max{t, σ} ≤ 4 γ  + δ  λ  / max{t, σ} ≤ 4n -1/2 ζ(t),
since by definition µ  ≤ 2γ  . In the same way, observing that

ξ  = κλ  2 ≤ γ  1 + 1 κ -1 we obtain ξ  + µ  + γ  + 4δ  λ  / max{t, σ} < 4 + (κ -1) -1 γ  + 4δ  λ  / max{t, σ} ≤ 6 + (κ -1) -1 n -1/2 ζ(max{t, σ})
and similarly,

2 µ  + γ  + 2δ  λ  / max{t, σ} ≤ 6n -1/2 ζ(max{t, σ}), 2 2γ  + δ  λ  / max{t, σ} ≤ 4n -1/2 ζ(max{t, σ}).
This implies that, whenever B * (t) < +∞, then (λ  , β  , t) ∈ Γ . We have then proved that inf

(λ,β)∈Λ B λ,β (t) ≤ B λ  , β  (t) ≤ n -1/2 ζ(max{t, σ}) 1 -4n -1/2 ζ(max{t, σ}) = B * (t).
We now need to prove the last part of the proposition. We first apply equation (1.28) to Proposition 1.16 and we obtain that, with probability at least 1 -2 , for any

θ ∈ R d , max{N (θ), σ θ 2 } max{ N (θ), σ θ 2 } -1 ≤ B * min θ -2 N (θ), s 2 4 .
Since θ -2 N (θ) ≤ s 2 4 , we conclude the proof.

Before introducing some comments on the bound presented in Proposition 1.17 we make more explicit its formulation by replacing the quantity s 2 4 by the trace of the Gram matrix.

Lemma 1.18. We have

s 2 4 ≤ √ κ Tr(G),
where Tr(G) = x 2 dP(x) denotes the trace of G.

Proof. By definition

s 4 4 = x 4 dP(x) = 1≤j≤d 1≤k≤d x 2 j x 2 k dP(x)
and by the Cauchy-Schwarz inequality we obtain

s 4 4 ≤ 1≤j≤d 1≤k≤d
x 4 j dP(x)

1 2 x 4 k dP(x) 1 2
.

Recalling that, by definition, κ = sup

θ∈R d θ, x 4 dP(x) ( θ, x 2 dP(x)) 2 , we get s 4 4 ≤ κ 1≤j≤d 1≤k≤d x 2 j dP(x) x 2 k dP(x) = κ   d j=1 x 2 j dP(x)   2 = κ x 2 dP(x) 2 ,
which concludes the proof.

We now apply the above lemma to Proposition 1.17. We recall that, given a > 0,

K = 1 + a -1 log n 72(2 + c)κ 1/2
where κ is defined in equation (1.10) on page 26 and c =

15 8 log(2)( √ 2-1) exp 1+2 √ 2 2 . Proposition 1.19. Let us fix a threshold σ ∈ R + such that σ ≤ s 2 4 . With probability at least 1 -2 , for any θ ∈ R d , max{N (θ), σ θ 2 } max{ N (θ), σ θ 2 } -1 ≤ B * θ -2 N (θ) ,
where

B * (t) =      n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ}) 6 + (κ -1) -1 ζ(max{t, σ}) ≤ √ n +∞ otherwise and ζ(t) = 2(κ -1) (2 + 3c) Tr(G) 4(2 + c)t + log(K/ ) cosh(a/4) + 2(2 + c)κ Tr(G) t cosh(a/2).
Let us provide some numerical evaluation of the constants. Let us choose a = 1/2. We can bound the logarithmic factor log log(n) hidden in log(K) with a relatively small constant assuming any reasonable upper bound for the number n of observations. For instance, if we assume that n ≤ 10 20 , and use the fact that κ ≥ 1 we can bound K with the constant 77, so that log(K) ≤ 4.35. In this case, ζ is bounded by

ζ(t) ≤ 2.032(κ -1) 0.73 Tr(G) t + 4.35 + log( -1 ) + 98.5 κ Tr(G) t , t ∈ R + .
We observe that the quantity κ -1 corresponds to a variance factor, and more precisely to

Var[ θ, X 2 ] N (θ) -2 .
In the special case when X is a Gaussian random vector we have κ = 3. Moreover, the bound B * does not depend explicitly on the dimension d and in particular the dimension is replaced by the entropy term Tr(G)/max θ -2 N (θ), σ . The remarkable fact is that, since the result is dimension-independent, it can be generalized to any infinite-dimensional Hilbert space, with the only assumption that the trace of the Gram matrix Tr(G) is finite.

For further details we refer to section 1.3.

Let us also make clear that we do not need to know the exact values of κ and Tr(G) to compute the estimator and evaluate the bound. We only need to know upper bounds for these two quantities. If we use those upper bounds in the definition of the estimator, the estimation error will hold true with κ and Tr(G) replaced with their upper bounds.

In order to have a meaningful (finite) bound we can choose for example the threshold σ such that 8

ζ(σ) = √ n,
so that B * (t) < +∞ for any t ∈ R + , assuming that we work with κ ≥ 3/2. With this choice the threshold decays to zero as the sample size grows to infinity. More precisely, using the inequality (

√ a + √ b) 2 ≤ 2(a + b), we get in this case that σ ≤ 100 κ Tr(G) n/128 -4.35 -log -1 .
(1.35)

A quadratic estimator of N

We have already said that our estimator N defined in equation (1.26) is unfortunately not a quadratic form. In this section we are going to deduce a quadratic estimator of the quadratic form N.

Let Λ ⊂ R + \ {0} 2 be the finite set defined by equation (1.27). Proposition 1.12 on page 28 provides a confidence region for N (θ). Define

B + (θ) = min (λ,β)∈Λ Φ -1 θ,+ N λ (θ) and B -(θ) = max (λ,β)∈Λ Φ θ,-N λ (θ) (1.36)
where we recall that N λ (θ) = λ/ α(θ) 2 and α(θ) is defined in equation (1.3) on page 17.

According to Proposition 1.12 we get the following result.

Proposition 1.20. With probability at least 1 -2 , for any

θ ∈ R d , B -(θ) ≤ N (θ) ≤ B + (θ).
In the following we are going to study the properties of a symmetric matrix Q that satisfies

Tr(Q 2 ) ≤ Tr(G 2 ) and B -(θ) ≤ θ Qθ ≤ B + (θ), θ ∈ Θ δ , (1.37)
where Θ δ is a finite δ-net of the unit sphere

S d = θ ∈ R d , θ = 1 , meaning that sup θ∈S d min ξ∈Θ δ θ -ξ ≤ δ.
Such a quadratic form exists with probability at least 1 -2 according to Proposition 1.20 and it can be computed using a convex optimization algorithm described in section 1.2.4.

Since, for any

θ ∈ S d , there is ξ ∈ Θ δ such that θ -ξ ≤ δ, we have θ Qθ -ξ Qξ = θ + ξ Q θ -ξ ≤ θ + ξ Q ∞ θ -ξ ≤ 2δ Tr(Q 2 ) ≤ 2δ
Φ -• Φ + θ Qθ -η ≤ N (θ) + η, Φ -• Φ + N (θ) -η ≤ θ Qθ + η,
where Φ -and Φ + are the values taken by Φ θ,-and Φ θ,+ , when θ ∈ S d .

We recall that the functions Φ θ,-and Φ θ,+ depend on θ through θ only.

Proof. According to Proposition 1.20, with probability at least 1 -2 , for any (λ, β) ∈ Λ,

θ Qθ ≤ Φ -1 + N λ (ξ) + η ≤ Φ -1 + • Φ -1 -N (ξ) + η.
Since equation (1.38) holds true also for N , we conclude that

θ Qθ ≤ Φ -1 + • Φ -1 -N (θ) + η + η.
In the same way

θ Qθ ≥ Φ -N λ (ξ) -η ≥ Φ -• Φ + N (ξ) -η ≥ Φ -• Φ + N (θ) -η -η.
We need some technical lemmas to take advantage of this result.

Lemma 1.22. Let us put

F (t) = max{t, σ}B * (min{t, s 2 4 }),
where B * is defined in Proposition 1.17. In the case when

κ ≥ 3/2 and 8ζ(σ) ≤ √ n, the function t → F (t) is non-decreasing for any t ∈ R + . Proof. If σ ≥ s 2 4 , then B * (min{t, s 2 4 }) = B * (σ), so that F (t) = max{t, σ}B * (σ) is obviously non-decreasing. Otherwise, σ ≤ s 2 4 , so that ζ max min{t, s 2 4 }, σ = ζ min max{t, σ}, s 2 4 .
Therefore the function F is of the form

F (t) = c ug(u) (1 -g(u))
,

where u = max{t, σ}, g(u) = a 1 /u + a 2 + a 3 /u, g(σ) ≤ 1/2
, and the constants c, a 1 , a 2 , and a 3 are positive. Let h(u) = a 1 /u + a 3 /u and observe that

g (u) = - 1 2u a 1 /u a 1 /u + a 2 1/2 + a 3 /u ≥ - 1 2u a 1 /u + a 3 /u = h (u)
and that g(u) ≥ h(u). Therefore h(u) ≤ g(u) ≤ 1/2, for any u ≥ σ, and

∂ ∂u log ug(u) 1 -g(u) = 1 u + g (u) g(u) 1 -g(u) ≥ 1 u + h (u) h(u) 1 -h(u) = 1 u - 1 2u 1 -h(u) ≥ 0,
showing that F is non-decreasing.

In the following we always assume that κ ≥ 3/2.

Lemma 1.23. For any (a, b) ∈ R 2 such that, for any (λ, β) ∈ Λ, Φ -• Φ + (a -η) ≤ b + η, and Φ -• Φ + (b -η) ≤ a + η,
(σ) ≤ σ, max Φ + (b -η), σ ≥ max Φ + (b -η), Φ + (σ) = Φ + max{(b -η), σ} = max{b -η, σ} 1 + B λ,β (b -
-η, σ} ≥ max Φ + (b -η), σ > max{a + η, σ}.
In this case

max{a + η, σ} ≥ max Φ -• Φ + (b -η), σ ≥ max Φ -• Φ + (b -η), Φ -(σ) ≥ Φ -max{Φ + (b -η), σ} ≥ max Φ + (b -η), σ 1 -B λ,β max Φ + (b -η), σ
according to equation (1.24). Moreover, continuing the above chain of inequalities

max{a + η, σ} ≥ max Φ + (b -η), Φ + (σ) 1 -B λ,β (max{a + η, σ}) = Φ + max{b -η, σ} 1 -B λ,β (a + η) ≥ max{b -η, σ} 1 -B λ,β (a + η) 1 + B λ,β (max{b -η, σ}) ≥ max{b -η, σ} 1 -B λ,β (a + η) 1 + B λ,β (a + η) . Therefore max{b -η, σ} ≤ max{a + η, σ} 1 + B λ,β (a + η) 1 -B λ,β (a + η) = max{a + η, σ} 1 + 2B λ,β (a + η) 1 -B λ,β (a + η) ≤ max{a + η, σ} 1 + 2 B λ,β (a + η)
according to Lemma 1.15. This concludes the proof of Step 1.

Step 2 Taking the infimum in (λ, β) ∈ Λ in equation (1.41), according to equation (1.28), we obtain that max{b -η, σ} ≤ max{a + η, σ} 1 + 2B * min{a + η, s 2 4 }) .

We can then use the fact that t → max{t, σ}B * (min{t, s 2 4 }) is non-decreasing (proved in Lemma 1.22) to deduce that max{b -η, σ} ≤ max{a + η, σ} + 2 max{b + η, σ}B * (min{b + σ, s 2 4 }), since there is nothing to prove when already max{b + η, σ} ≤ max{a + η, σ}. Remark that max{a + η, σ} ≤ max{a + η, σ + η} ≤ max{a, σ} + η and that in the same way max{b -η, σ} ≥ max{b, σ} -η. This proves that max{b, σ} -max{a, σ} ≤ 2 max{a + η, σ}B * min{a + η, s 2 4 } + 2η and max{b, σ} -max{a, σ} ≤ 2 max{b + η, σ}B * min{b + η, s 2 4 } + 2η. By symmetry, we can then exchange a and b to prove the same bounds for max{a, σ} -max{b, σ}, and therefore also for the absolute value of this quantity, which ends the proof of the lemma. Lemma 1.24. For any (a, b) ∈ R 2 such that, for any (λ, β) ∈ Λ,

Φ -• Φ + (a -η) ≤ b + η,
and

Φ -• Φ + (b -η) ≤ a + η,
and any threshold σ ∈ R + such that 8ζ(σ) ≤ √ n and σ ≤ s 2 4 , we have max{a, σ} -max{b, σ} ≤ 2 max{a, σ}B * min{a, s 2 4 } + 5η/2 max{a, σ} -max{b, σ} ≤ 2 max{b, σ}B * min{b, s 2 4 } + 5η/2. Proof. This is a consequence of the previous lemma, of the fact that B * (min{t, s 2 4 }) ≤ 1/4, and of the fact that max{a + η, σ} ≤ max{a, σ} + η.

We deduce from Lemma 1.21 and Lemma 1.24 the analogous for the quadratic estimator θ Qθ of the dimension-free bound presented in Proposition 1.17 for N (θ). There is still some improvement to bring: at this stage we are not sure that Q is nonnegative. Nevertheless we know that, for any

θ ∈ S d , there is ξ ∈ Θ δ such that θ -ξ ≤ δ, so that θ Qθ ≥ ξ Qξ -η ≥ B -(ξ) -η ≥ -η.
Decomposing Q into its positive and negative parts and writing

Q = Q + -Q -, we deduce that Q -∞ = sup θ∈S d θ Q -θ = -inf θ∈S d θ Qθ ≤ η,
where we recall that η = 2δ Tr(G 2 ). Therefore, for any

θ ∈ S d , max{θ Qθ, σ} -max{θ Q + θ, σ} ≤ θ Qθ -θ Q + θ = θ Q -θ ≤ η.
We have proved the following result.

max{θ Q + θ, σ} -max{N (θ), σ} ≤ 2 max N (θ), σ B * N (θ) + 7δ Tr(G 2 ), max{θ Q + θ, σ} -max{N (θ), σ} ≤ 2 max θ Q + θ, σ B * min{θ Q + θ, s 2 4 + 7δ Tr(G 2 ),
where B * is defined in Proposition 1.17.

Let us recall that we can obtain a bound stated in terms of Tr(G), using the inequality s 2 4 ≤ κ 1/2 Tr(G) (proved in Lemma 1.18) that implies

ζ(t) ≤ 2(κ -1) (2 + 3c) Tr(G) 4(2 + c)t + log(K/ ) cosh(a/4)+ 2(2 + c)κ Tr(G) t cosh(a/2). Remark that, for any a, b ∈ R + , a -b ≤ max{a, σ} -max{b, σ} + σ, so that |a -b| ≤ max{a, σ} -max{b, σ} + σ. Corollary 1.27. Assume that 8ζ(σ) ≤ √ n, σ ≤ s 2 4 and that κ ≥ 3/2. With probability at least 1 -2 , for any θ ∈ S d , θ Q + θ -N (θ) ≤ 2 max{N (θ), σ}B * N (θ) + 7δ Tr(G 2 ) + σ θ Q + θ -N (θ) ≤ 2 max{θ Q + θ, σ}B * min θ Q + θ, s 2 4 + 7δ Tr(G 2 ) + σ.
We conclude the section by providing two more bounds on the approximation error |N (θ)θ Q + θ| that depend on the largest eigenvalue of N and on the largest eigenvalue of Q + respectively. Applying Lemma 1. [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF]

N (θ) -θ Q + θ ≤ 2 max{ G ∞ , σ}B * ( G ∞ ) + 7δ G F + σ, N (θ) -θ Q + θ ≤ 2 max{ Q + ∞ , σ}B * min Q + ∞ , s 2 4 + 7δ G F + σ,
where G F = Tr(G 2 ) is the Frobenius norm of the Gram matrix.

How to compute Q

In this section we provide an explicit construction of the quadratic form Q defined in equation (1.37) on page 40, inspired by [START_REF] Catoni | Estimating the gram matrix and least square regression through pac-bayes bounds[END_REF]. Later, in section 1.6 we will present a simplified construction that leads to an estimator which is not exactly the same as Q but presents the same kind of behavior.

From a theoretical point of view, to obtain a quadratic estimator of N it is sufficient to choose the quadratic form associated with any non-negative symmetric matrix Q satisfying

B -(θ) ≤ θ Qθ ≤ B + (θ), θ ∈ S d , ( 1
Q = θ∈Θ ξ + (θ) -ξ -(θ) θθ ,
where

ξ + , ξ -= arg max ξ + ,ξ -∈(R 2 + ) Θ - 1 2 (θ,θ )∈Θ 2 [ξ + (θ) -ξ -(θ)] ξ + (θ ) -ξ -(θ ) θ, θ 2 + θ∈Θ [ξ + (θ)B -(θ) -ξ -(θ)B + (θ)] ,
is such that

Q = arg min H F | H = H, B -(θ) ≤ θ Hθ ≤ B + (θ), θ ∈ Θ . Moreover Q F ≤ G F ≤ Tr(G), and
|θ 2 Qθ 2 -θ 1 Qθ 1 | ≤ θ 1 + θ 2 G F θ 1 -θ 2 , θ 1 , θ 2 ∈ R d . (1.44)
Proof. We define

H = arg min H F | H = H, B -(θ) ≤ θ Hθ ≤ B + (θ), θ ∈ Θ
and we observe that it is the solution of the (minimax) optimization problem

V = inf H =H sup ξ + ,ξ - V (H, ξ + , ξ -),
where

V (H, ξ + , ξ -) = 1 2 H 2 F + θ∈Θ ξ + (θ) B -(θ) -θ Hθ + ξ -(θ) θ Hθ -B + (θ) .
On the event of probability 1 -2 described in Proposition 1.20 where

B -(θ) ≤ θ Gθ ≤ B + (θ),
the constraints are satisfied when H = G and therefore, since the constraints are linear, by the Slater condition, inf

H sup ξ + ,ξ - V (H, ξ + , ξ -) = sup ξ + ,ξ - inf H V (H, ξ + , ξ -).
We can compute explicitly the solution of inf

H V (H, ξ + , ξ -)
and we obtain inf

H V (H, ξ + , ξ -) = V (Q, ξ + , ξ -). It follows that H = Q, since sup ξ + ,ξ -inf H V (H, ξ + , ξ -) = V (Q, ξ + , ξ -).
Moreover, for any θ 1 , θ 2 ∈ R d , we have

|θ 2 Qθ 2 -θ 1 Qθ 1 | ≤ (θ 1 + θ 2 ) Q (θ 2 -θ 1 ) ≤ θ 1 + θ 2 Q ∞ θ 2 -θ 1 .
We conclude observing that

Q ∞ ≤ Q F ≤ G F . The inequality G F ≤ Tr(G)
is a consequence of the fact that the eigenvalues of G are non-negative, so that

G 2 F = Tr(G 2 ) ≤ Tr(G) 2 .

Gram operators in Hilbert spaces

We have already underlined the fact that the results presented in the previous section do not explicitly depend on the dimension d of the ambient space and that, for this reason, they can be generalized to any separable Hilbert space. In this section we study how these results extend to infinite-dimensional Hilbert spaces.

Let H be a separable Hilbert space and let P ∈ M 1 + (H) be an unknown probability distribution on H. We consider the Gram operator G : H → H defined by

Gθ = θ, v H v dP(v)
and we assume Tr(G) = E X 2 H < +∞, where X ∈ H denotes a random vector of law P. In analogy to the previous section we denote by N the quadratic form associated with G so that

N (θ) = Gθ, θ H = θ, v 2 H dP(v), θ ∈ H.
We consider (H k ) k an increasing sequence of subspaces of H of finite dimension such that ∪ k H k = H and we endow each space H k with the probability P k which arises from the disintegration of P. We denote by N k the quadratic form in H k associated with the probability distribution P k and we observe that, denoting by Π k the orthogonal projector on H k , for any θ ∈ H, we have

N k (Π k θ) = N (Π k θ) .
In the following we consider i.i.d. samples of size n in H drawn according to P. According to Proposition 1.12 on page 28, the event

A k = ∀θ ∈ H k , ∀(λ, β) ∈ Λ, Φ θ,- λ α(θ) 2 ≤ N (θ) ≤ Φ -1 θ,+ λ α(θ) 2 is such that P ⊗n A k ≥ 1 -2 . Since A k+1 ⊂ A k , by the continuity of measure, P ⊗n k∈N A k ≥ 1 -2 .
This means that with probability at least 1 -2 , for any θ ∈ k H k and any (λ,

β) ∈ Λ, Φ θ,- λ α(θ) 2 ≤ N (θ) ≤ Φ -1 θ,+ λ α(θ) 2 .
Consequently, since N (θ) = lim k→+∞ N Π k (θ) , for any θ ∈ H, the following result holds. Proposition 1.30. With probability at least 1 -2 , for any θ ∈ H,

B -(θ) ≤ N (θ) ≤ B + (θ)
where

B -(θ) = lim sup k→+∞ sup (λ,β)∈Λ Φ Π k (θ),- λ α(Π k (θ)) 2 , B + (θ) = lim inf k→+∞ inf (λ,β)∈Λ Φ -1 Π k (θ),+ λ α(Π k (θ)) 2 .
If we do not want to go to the limit, we can use the explicit bound

N (θ) -N (Π k (θ)) = θ + Π k (θ), G θ -Π k (θ) H | ≤ 2 θ H G ∞ θ -Π k (θ) H ≤ 2 θ H Tr(G) θ -Π k (θ) H = 2 θ H E X 2 H θ -Π k (θ) H .
This bound depends on θ -Π k θ H . We will see in the following another bound that goes uniformly to zero for any θ ∈ S H when k tends to infinity. In the same way, proceeding as already done in section 1.2 we state the analogous of Proposition 1.17.

Proposition 1.31. Let κ ≥ sup θ∈H E( θ,X 2 H )>0 E θ, X 4 H E θ, X 2 H 2 and s 4 ≥ E X 4 H 1/4
be known constants. Let a > 0 and

K = 1 + a -1 log n 72(2 + c)κ 1/2 with c = 15 8 log(2)( √ 2-1) exp 1+2 √ 2 2
. Define

ζ(t) = 2(κ -1) (2 + 3c)s 2 4 4(2 + c)κ 1/2 t + log(K/ ) cosh(a/4) + 2(2 + c)κ 1/2 s 2 4 t cosh(a/2)
and consider, according to equation (1.26) on page 32, the estimator

N (θ) = N λ (θ), θ ∈ k H k .
For any θ ∈ H, define N (θ) by choosing any accumulation point of the sequence N Π k (θ) . Define the bound

B * (t) = n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ}) ,
where σ ∈]0, s 2 4 ] is some energy level such that

6 + (κ -1) -1 ζ(max{t, σ}) ≤ √ n.
With probability at least 1 -2 , for any θ ∈ H,

max{N (θ), σ θ 2 H } max{ N (θ), σ θ 2 H } -1 ≤ B * θ -2 H N (θ) .
Proof. This is a consequence of the fact that lim We recall that if we choose the threshold σ such that 8ζ(σ) = √ n, then σ goes to zero as the sample size grows to infinity as shown in equation (1.35) on page 39.

Let X 1 , . . . , X n ∈ H be an i.i.d. sample drawn according to P. Define V = span X 1 , . . . , X n and

V k = span Π k X 1 , . . . , Π k X n = Π k (V ).
Let Θ δ be a δ-net of S H ∩ V k (where S H denotes the unit sphere in H). Remark that Θ δ is finite because dim(V k ) ≤ n < +∞. We can compute as in Lemma 1.29 on page 45 a linear operator

G k : V k → V k such that Tr( G 2 k ) ≤ Tr(G 2 ) and Φ -N λ (θ) ≤ G k θ, θ H ≤ Φ -1 + N λ (θ) , θ ∈ Θ δ .
Consider the estimator of G defined as

Q = G k • Π V k , (1.45)
where Π V k is the orthogonal projector on V k . For any θ ∈ S H ,

θ, Qθ H = Π V k θ, QΠ V k θ H ≤ Π V k θ 2 H ξ, Qξ H + η , where ξ ∈ Θ δ is the closest point in Θ δ to Π V k θ -1 H Π V k θ and where η = 2δ Tr(G 2 ). Since ξ ∈ H k , with probability at least 1 -, for any (λ, β) ∈ Λ, ξ, Qξ H ≤ Φ -1 + N λ (ξ) = Φ -1 + N λ ξ + Π V k θ -1 H Π k -Π V k θ .
Let us now remark that for any a ∈ [0, 1], we have Φ + (at) ≤ aΦ + (t), so that aΦ -1

+ (t) ≤ Φ -1 + (at). Let us also remark that N λ (aθ) = a 2 N λ (θ). Therefore θ, Qθ H ≤ Π V k θ 2 H Φ -1 + N λ Π V k θ -1 H Π V k θ H ξ + Π k -Π V k θ + η ≤ Π V k θ 2 H Φ -1 + • Φ -1 -N Π V k θ -1 H Π V k θ H ξ + Π k -Π V k θ + η ≤ Φ -1 + • Φ -1 -N Π V k θ H ξ + Π k -Π V k θ + η ≤ Φ -1 + • Φ -1 -N Π k θ + η + η. Indeed, Π V k θ H ξ + Π k -Π V k θ -Π k θ ≤ δ,
and this is a difference of two vectors belonging to the unit ball. In the same way

θ, Qθ H ≥ Π V k θ 2 H ξ, Qξ H -η ≥ Π V k θ 2 H Φ -N λ Π V k θ -1 H Π V k θ H ξ + Π k -Π V k θ -η ≥ Π V k θ 2 H Φ -• Φ + N Π k θ -1 H Π V k θ H ξ + Π k -Π V k θ -η ≥ Φ -• Φ + N Π k θ -η -η.
Let us decompose Q in its positive and negative parts and write

Q = Q + -Q -.
max θ, Q + θ H , σ -max Π k θ, GΠ k θ H , σ ≤ 2 max Π k θ, GΠ k θ H , σ B * Π k θ, GΠ k θ H + 7δ Tr(G 2 ) max θ, Q + θ H , σ -max Π k θ, GΠ k θ H , σ ≤ 2 max θ, Q + θ H , σ B * min{ θ, Q + θ H , s 2 4 } + 7δ Tr(G 2 ).
Let us consider {p i } +∞ i=1 an orthonormal basis of eigenvectors of G such that the corresponding sequence of eigenvalues {λ i , i = 1, . . . , +∞} is non-increasing. For any θ ∈ S H , we have

θ, Gθ H -Π k θ, GΠ k θ H = +∞ i=1 Π k θ, p i 2 H -θ, p i 2 H λ i = +∞ i=1 θ, Π k p i 2 H -θ, p i 2 H λ i = +∞ i=1 θ, Π k p i + p i H θ, p i -Π k p i H λ i ≤ inf m=1,...,+∞ m-1 i=1 2λ i p i -Π k p i H + λ m . Indeed, +∞ i=m Π k θ, p i 2 H ≤ 1, so that +∞ i=m Π k θ, p i 2 H λ i ≤ sup i=m,...,+∞ λ i = λ m ,
and in the same way

+∞ i=m θ, p i 2 H λ i ≤ λ m .
We deduce the following result.

Proposition 1.33. Consider some threshold σ ∈ R + such that σ ≤ s 2 4 and 8ζ(σ) ≤ √ n. Assume that κ ≥ 3/2. With probability at least 1 -2 , for any θ ∈ S H , for any k,

max θ, Q + θ H , σ -max θ, Gθ H , σ ≤ 2 max θ, Gθ H , σ B * θ, Gθ H + 7δ Tr(G 2 ) + 3ν k max θ, Q + θ H , σ -max θ, Gθ H , σ ≤ 2 max θ, Q + θ H , σ B * min{ θ, Q + θ H , s 2 4 } + 7δ Tr(G 2 ) + 2ν k ,
where B * is defined in Proposition 1.31 and

ν k = inf m=1,...,+∞ m-1 i=1 λ i p i -Π k p i H + λ m /2 ≤ inf m=1,...,+∞ m-1 i=1 λ i p i -Π k p i H + Tr(G)/(2m) -→ k→+∞ 0.
Remark that we can use this result to bound θ, (G -Q + )θ H , using the inequality

θ, (G -Q + )θ H ≤ max θ, Q + θ H , σ -max θ, Gθ H , σ + σ.
Let us mention to conclude this section another way to extend the results from finite dimension to separable Hilbert spaces. This second method consists in defining directly a Gaussian perturbation of the parameter in the Hilbert space. Indeed, taking a basis, we can identify any separable Hilbert space H with the sequence space 2 ⊂ R N . We can take as our prior parameter distribution π 0 the law on R N of a sequence of independent onedimensional Gaussian random variables with mean 0 and variance 1/β. For any θ ∈ 2 , we can then define the posterior π θ as i∈N N (θ i , 1/β) ∈ M 1 + (R N ). One can prove that

K(π θ , π 0 ) = +∞ i=1 K N ( θ i , β -1 ), N (0, β -1 ) = +∞ i=1 β 2 θ 2 i = β θ 2 H 2 .
Moreover, when θ ∼ π θ ,

W (x) = lim n→+∞ n i=0 θ i x i , x ∈ 2 ,
exists π θ -almost surely and is a Gaussian process indexed by 2 H, with mean

E W (x) = θ, x H , x ∈ H,
and covariance

E W (x)W (y) -E W (x) E W (y) = 1 β x, y H , x, y ∈ H.
Using these properties, we can directly generalize to 2 H the PAC-Bayesian bounds stated in R d in the previous section.

Symmetric random matrices 1.4.1 Preliminaries

In this section we generalize the previous results to estimate the expectation of a symmetric random matrix. Indeed the Gram matrix can be viewed as the expectation of the symmetric positive semi-definite random matrix XX , where X ∈ R d is a random vector.

Let A ∈ M d (R) be a symmetric random matrix of size d. As already observed for the Gram matrix, the expectation of A can be recovered via the polarization identity from the quadratic form

N A (θ) = E θ Aθ , θ ∈ R d ,
where the expectation is taken with respect to the unknown probability distribution of A on the space of symmetric matrices of size d. We decompose A = U DU in its positive and negative parts

A = A + -A -,
where A + and A -are symmetric positive semi-definite random matrices defined by

A + = U D + U and A -= U D -U ,
and D + (respectively D -) is the diagonal matrix whose entries are the positive (respectively negative) parts of those of D. According to the previous decomposition, the quadratic form N A rewrites as

N A (θ) = E θ Aθ = E θ A + θ -E θ A -θ = N A + (θ) -N A -(θ).
Thus in the following we will consider the case where A ∈ M d (R) is a symmetric positive semi-definite random matrix of size d.

Symmetric positive semi-definite random matrices

Let A ∈ M d (R) be a symmetric positive semi-definite random matrix of size d and let P be the (unknown) probability distribution of A on the space of symmetric positive semi-definite matrices of size d. Our goal is to estimate the quadratic form

N (θ) = E θ Aθ , θ ∈ R d ,
where E is the expectation with respect to P, from an i.i.d. sample A 1 , . . . , A n ∈ M d (R) of symmetric positive semi-definite matrices drawn according to P.

Since A can be decomposed as

A = U DU = U D 1/2 U 2 ,
where D ∈ M d (R) is the diagonal matrix of the eigenvalues of A and U ∈ M d (R) is the orthogonal matrix of its eigenvectors, the quadratic form N rewrites as

N (θ) = E A 1/2 θ 2
where

A 1/2 = U D 1/2 U is the square root of A.
The construction of the robust estimator of the quadratic form N follows the one already done in the case of the Gram matrix with the necessary adjustments.

For any λ > 0 and for any θ ∈ R d , we define

r λ (θ) = 1 n n i=1 ψ A 1/2 i θ 2 -λ ,
where the influence function ψ is defined as in equation (1.2) on page 16.

We use a PAC-Bayesian approach linked with Gaussian perturbations to define and study a robust estimator of N .

We perturb θ with the Gaussian perturbation π θ ∼ N (θ, β -1 I) of mean θ and covariance matrix β -1 I, where β > 0 is a free real parameter to be determined. The following result holds.

Lemma 1.34. We have

A 1/2 θ 2 dπ θ (θ ) = A 1/2 θ 2 + Tr(A) β .
Proof. Let W ∈ R d be a Gaussian random variable with mean A 1/2 θ and covariance matrix β -1 A. We have

E A 1/2 θ 2 dπ θ (θ ) = E W 2 = d i=1 E W, e i 2
where {e i } d i=1 is the canonical basis of R d . Since W, e i is a one-dimensional Gaussian random variable with mean A 1/2 θ, e i and variance β -1 e i Ae i , we conclude that

E A 1/2 θ 2 dπ θ (θ ) = d i=1 Var W, e i + E W, e i 2 = 1 β Tr(A) + A 1/2 θ 2 .
According to Lemma 1.34 the empirical criterion r λ rewrites as

r λ (θ) = 1 n n i=1 ψ A 1/2 i θ 2 - Tr(A i ) β -λ dπ θ (θ ) .
With the help of Lemma 1.3 on page 20, we pull the expectation outside the influence function and we get

ψ A 1/2 θ 2 -λ ≤ χ A 1/2 θ 2 - Tr(A) β -λ dπ θ (θ ) + min log(4), 1 8 Var A 1/2 θ 2 dπ θ (θ )
where the function χ is defined in equation (1.4) on page 18 and is such that

χ(z) ≤ log 1 + z + z 2 /2 , z ∈ R.
The following lemma provides a bound on the variance.

Lemma 1.35. We have

Var A 1/2 θ 2 dπ θ (θ ) = 4 β Aθ 2 + 2 β 2 Tr(A 2 ).
Proof. Let us decompose A into A = U DU , where U U = I and D = diag(λ 1 , . . . , λ d ).

Remark that U θ has the same distribution as U θ + W , where W ∼ N (0, β -1 I) so that

Var A 1/2 θ 2 dπ θ (θ ) = Var d i=1 (U θ) i + W i 2 λ i = d i=1 λ 2 i Var (U θ) i + W i 2 = d i=1 2 β 2 + 4 β U θ 2 i λ 2 i = 2 β 2 Tr(A 2 ) + 4 β Aθ 2 .
As a consequence, we obtain that

ψ A 1/2 θ 2 -λ ≤ χ A 1/2 θ 2 - Tr(A) β -λ dπ θ (θ ) + min log(4), 1 2β Aθ 2 + Tr(A 2 ) 4β 2 .
We now use Lemma 1.4 on page 20 with m = Aθ , a = log(4), b = 1/(2β) and c = Tr(A 2 )/(4β 2 ) to obtain

ψ A 1/2 θ 2 -λ ≤ χ A 1/2 θ 2 - Tr(A) β -λ dπ θ (θ ) + min 4 log(2), 1 β Aθ 2 + Tr(A 2 ) 2β 2 dπ θ (θ )
and we conclude, by Lemma 1.5, that

ψ A 1/2 θ 2 -λ ≤ log 1 + A 1/2 θ 2 - Tr(A) β -λ + 1 2 A 1/2 θ 2 - Tr(A) β -λ 2 + c β Aθ 2 + Tr(A 2 ) 2β dπ θ (θ ),
where c =

15 8 log(2)( √ 2-1) exp 1+2 √ 2 2
.

In terms of the empirical criterion r λ we have proved that

r λ (θ) ≤ 1 n n i=1 log 1 + A 1/2 i θ 2 - Tr(A i ) β -λ + 1 2 A 1/2 i θ 2 - Tr(A i ) β -λ 2 + c β A i θ 2 + Tr(A 2 i ) 2β dπ θ (θ ).
We use the PAC-Bayesian inequality presented in Proposition 1.7 on page 22 where the family of posterior distribution is

π θ ∼ N (θ, β -1 I) | θ ∈ R d , β > 0
and the prior distribution is π 0 ∼ N (0, β -1 I), to provide an upper bound for r λ which is uniform in θ.

Proposition 1.36. Define t(A, θ ) = A 1/2 θ 2 - Tr(A) β -λ. With probability at least 1 -, for any θ ∈ R d , r λ (θ) ≤ E t(A, θ ) + 1 2 t(A, θ ) 2 + c β Aθ 2 + Tr(A 2 ) 2β dπ θ (θ ) + β θ 2 2n + log( -1 ) n = E A 1/2 θ 2 -λ + 1 2 A 1/2 θ 2 -λ 2 + (c + 2) Aθ 2 β + (2 + 3c) Tr(A 2 ) 2β 2 + β θ 2 2n + log( -1 ) n .
Proof. The proof is essentially the same of the one of Proposition 1.8 on page 25.

In analogy to the definitions given in equation (1.10) on page 26, we introduce κ = sup

θ∈R d E( A 1/2 θ 2 )>0 E A 1/2 θ 4 E A 1/2 θ 2 2 .
Using the Cauchy-Schwarz inequality, remark that

E Aθ 2 ≤ E A ∞ A 1/2 θ 2 ≤ E A 2 ∞ 1/2 E A 1/2 θ 4 1/2 ≤ E A 2 ∞ 1/2 κ 1/2 E A 1/2 θ 2 . that rewrites as E Aθ 2 ≤ E A 2 ∞ 1/2 κ 1/2 N (θ). (1.46)
We then deduce from Proposition 1.36 the following result.

Proposition 1.37. With probability at least 1 -, for any θ ∈ R d ,

r λ (θ) ≤ κ 2 N (θ) -λ 2 + 1 + (κ -1)λ + (2 + c)κ 1/2 E A 2 ∞ 1/2 β N (θ) -λ + (κ -1)λ 2 2 + (2 + c)κ 1/2 E A 2 ∞ 1/2 λ β + (2 + 3c)E Tr(A 2 ) 2β 2 + β θ 2 2n + log( -1 ) n .
Remark that if we consider instead of equation (1.46) any upper bound of the form

E Aθ 2 ≤ f E(A) N (θ), Proposition 1.37 holds replacing E A 2 ∞ 1/2 κ 1/2 with f E(A)
. Similarly we can replace E Tr(A 2 ) with an upper bound.

We also observe that Proposition 1.37 is the analogous in the case of a symmetric positive semi-definite random matrix of Proposition 1.10 on page 27. Hence according to section 1.2.2 we define a robust estimator as follows.

Let Λ ⊂ R + \ {0} 2 be a finite set of possible values of the couple of parameters (λ, β).

We consider the family of estimators

N λ (θ) = λ α(θ) 2 , (λ, β) ∈ Λ, where α(θ) = sup {α ∈ R + | r λ (αθ) ≤ 0} and we define ( λ, β) = arg min (λ,β)∈Λ B λ,β θ -2 N λ (θ) ,
where B λ,β is defined in equation (1.20) on page 31. We consider as an estimator

N (θ) = N λ (θ), θ ∈ R d .
(1.47)

We now directly state the analogous of Proposition 1.19 on page 38, choosing the grid

Λ = (λ j , β j ) | 0 ≤ j < K ,
where a > 0 and

K = 1 + a -1 log n 72(2 + c)κ 1/2 , λ j = 2 (κ -1)n (2 + 3c)E Tr(A 2 ) 4(2 + c)κ 1/2 E A 2 ∞ exp(ja) + log K/ , β j = 2(2 + c)κ 1/2 E A 2 ∞ exp -(j -1/2))a . Proposition 1.38. Let σ ∈ R + be some energy level such that σ ≤ E A 2 ∞ 1/2 . With probability at least 1 -2 , for any θ ∈ R d , max{N (θ), σ θ 2 } max{ N (θ), σ θ 2 } -1 ≤ B * θ -2 N (θ) ,
where B * is defined as

B * (t) =      n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ}) 6 + (κ -1) -1 ζ(max{t, σ}) ≤ √ n +∞ otherwise and ζ(t) = 2(κ -1) (2 + 3c)E[Tr(A 2 )] 4(2 + c)κ 1/2 E A 2 ∞ 1/2 t + log(K/ ) cosh(a/4) + 2(2 + c)κ 1/2 E A 2 ∞ 1/2 t cosh(a/2).
We conclude this section by making some comments on the above result.

As already observed in the case of the Gram matrix, the quantity κ -1 corresponds to a variance term, and more precisely to

Var A 1/2 θ 2 N (θ) -2 .
Moreover the bound does not depend explicitly on the dimension d and in particular the dimension is replaced by the entropy terms

E Tr(A 2 ) κ 1/2 E A 2 ∞ 1/2 max { θ -2 N (θ), σ} and E A 2 ∞ 1/2 κ 1/2 max { θ -2 N (θ), σ} .
Let us remark that they satisfy

E A 2 ∞ 1/2 κ 1/2 ≤ E Tr(A 2 ) κ 1/2 E A 2 ∞ 1/2 ≤ E Tr(A) . (1.48) Indeed, A 2 ∞ ≤ Tr(A 2 ), E Tr(A 2 ) ≤ E A ∞ Tr(A) ≤ E A 2 ∞ 1/2 E Tr(A) 2 1/2 ,
and, denoting by

{e i } d i=1 an orthonormal basis of R d , E Tr(A) 2 = 1≤i≤d, 1≤j≤d E A 1/2 e i 2 A 1/2 e j 2 ≤ 1≤i≤d, 1≤j≤d E A 1/2 e i 4 1/2 E A 1/2 e j 4 1/2 ≤ κ 1≤i≤d, 1≤j≤d E A 1/2 e i 2 E A 1/2 e j 2 = κE Tr(A) 2 .
As a consequence, ζ can be bounded as shown

ζ(t) ≤ 2(κ -1) (2 + 3c)E[Tr(A)] 4(2 + c)t + log(K/ ) cosh(a/4) + 2(2 + c)κ E[Tr(A)] t cosh(a/2).
Let us remark also that in the case when rank(A) = 1 (as it is for the Gram matrix) then A 2 ∞ = Tr(A 2 ), so that the two entropy terms are the same. As the entropy terms are dominated by E Tr(A) , the result can be generalized to the case where A is a random symmetric positive semi-definite operator in an infinite-dimensional Hilbert space with the only additional assumption that E[Tr(A)] < +∞. For more detail we refer to section 1.3.

We also recall that if we choose the threshold σ such that 8ζ(σ) = √ n, so that B * (t) < +∞ for any t ∈ R + , then σ goes to zero as the sample size grows to infinity as shown in equation (1.35).

We conclude by recalling that the estimator N is not a quadratic form but that it is possible to deduce a positive semi-definite quadratic estimator of N as explained in section 1.2.3 on page 39.

Covariance matrix

Let X ∈ R d be a random vector distributed according to the unknown probability measure

P ∈ M 1 + (R d ). The covariance matrix of X Σ = E X -E(X) X -E(X)
is a symmetric positive semi-definite random matrix. Our goal is to estimate, uniformly in θ, the quadratic form associated with Σ

N (θ) = E θ, X -E(X) 2 , θ ∈ R d ,
from an i.i.d. sample X 1 , . . . , X n ∈ R d drawn according to P. We cannot use the results we have proved for the Gram matrix, since the quadratic form N depends on the unknown quantity E(X). However we can find a workaround, using the results of the previous section about random symmetric matrices. Indeed, the quadratic form N can be written as

N (θ) = 1 2 E θ, X -X 2
where X is an independent copy of X. More generally, given q ∈ N, we may consider q independent copies X (1) , . . . , X (q) of X and the random matrix

A = 1 q(q -1) 1≤j<k≤q X (j) -X (k) X (j) -X (k) .
Remark that

N (θ) = 1 q(q -1) E 1≤j<k≤q θ, X (j) -X (k) 2 = E θ A θ .
Consider the sample A 1 , . . . , A n/q of independent copies of A defined as

A i = 1 q(q -1) (i-1)q<j<k≤iq (X j -X k )(X j -X k ) .
We can use the results of the previous section to define a robust estimator of N (θ). We will discuss later how to choose q.

Let us introduce

κ = sup θ∈R d , E( A 1/2 θ 2 )>0 E A 1/2 4 E A 1/2 θ 2 2 ,
and κ = sup

θ∈R d , E( θ,X-E(x) 2 )>0 E θ, X -E(X) 4 E θ, X -E(X) 2 2 .
Lemma 1.39. The two kurtosis coefficients introduced above are related by the relation κ ≤ 1 + τ q (κ)/q, where τ q (κ) = κ -1 + 2 q -1 .

Proof. Replacing X with X -E(X) we may assume during the proof that E(X) = 0.

Recalling the definition of covariance, we have

E A 1/2 θ 4 = E 1 q(q -1) 1≤j<k≤q θ, X (j) -X (k) 2 2 = 1 q 2 (q -1) 2 1≤j<k≤q 1≤s<t≤q E θ, X (j) -X (k) 2 θ, X (s) -X (t) 2 = 1 q 2 (q -1) 2 1≤j<k≤q 1≤s<t≤q E θ, X (j) -X (k) 2 E θ, X (s) -X (t) 2 + 1≤j<k≤q E θ, X (j) -X (k) 4 -E θ, X (j) -X (k) 2 2 + 1≤j<k≤q 1≤s<t≤q |{j,k}∩{s,t}|=1 E θ, X (j) -X (k) 2 θ, X (s) -X (t) 2 -E θ, X (j) -X (k) 2 E θ, X (s) -X (t) 2 = 1 4 E θ, X (2) -X (1) 2 2 + 1 2q(q -1) E θ, X (2) -X (1) 4 -E θ, X (2) -X (1) 2 2 + q -2 q(q -1) E θ, X (1) -X (2) 2 θ, X (1) -X (3) 2 -E θ, X (1) -X (2) 2 2
Define W j = θ, X (j) .

E (W 1 -W 2 ) 2 2 = 4N (θ) 2 , E (W 1 -W 2 ) 4 = E W 4 1 + 6E W 2 1 E W 2 2 + E W 4 2 = 2E W 4 1 + 6E W 2 1 2 ≤ (2κ + 6)N (θ) 2 , E (W 1 -W 2 ) 2 (W 1 -W 3 ) 2 = E W 4 1 + 3E W 2 2 ) 2 ≤ (κ + 3)N (θ) 2 . Therefore E A 1/2 θ 4 ≤ 1 + (q -2)(κ -1) q(q -1) + (κ + 1) q(q -1) N (θ) 2 = 1 + τ q (κ) q N (θ) 2 , so that κ ≤ 1 + τ q (κ)/q, since E A 1/2 θ 2 = N (θ).
Remark that

E Tr(A) = Tr E(A) = Tr(Σ) = E X -E(X) 2
max{N (θ), σ θ 2 } max{ N (θ), σ θ 2 } -1 ≤ B * θ -2 N (θ) ,
where

B * (t) =        q n/q 1/2 ζ max{t, σ}) 1 -4 q n/q 1/2 ζ(max{t, σ}) , if 6 + q/τ q (κ) ζ max{t, σ} ≤ q n/q 1/2 , +∞, otherwise, ζ(t) = 2τ q (κ) (2 + 3c) Tr(Σ) 4(2 + c)t + log K/ cosh(a/4) + 2(2 + c) q + τ q (κ) Tr(Σ) t cosh(a/2), K = 1 + a -1 log n/q 72(2 + c) 1 + τ q (κ)/q 1/2 , τ q (κ) = κ -1 + 2 q -1 and c = 15 8 log(2)( √ 2-1) exp 1+2 √ 2 2 .
Proof. The proof follows from Proposition 1.38 replacing κ by κ and n by n/q . Indeed, using equation (1.48), we get

ζ(t) ≤ 2(κ -1) (2 + 3c)E[Tr(A)] 4(2 + c)t + log(K/ ) cosh(a/4) + 2(2 + c)κ E[Tr(A)] t cosh(a/2).
Since by Lemma 1.39 we see that κ -1 ≤ τ q (κ)/q, we conclude the proof.

As stated at the beginning of this proposition, we have used here rather crude entropy bounds. We can improve the constants by evaluating more carefully E Aθ 2 and E Tr(A 2 ) .

Lemma 1.41. We have

E Aθ 2 ≤ 1 - q -2 q(q -1) Σ ∞ N (θ) + 1 q κ + 1 q -1 Tr(Σ)N (θ) (1.49) E Tr A 2 ≤ 1 - q -2 q(q -1) Tr Σ 2 + 1 q κ + 1 q -1 Tr(Σ) 2 . (1.50) Proof. Replacing X with X -E(X) we may assume that E(X) = 0. Recall that E X 4 ≤ κE X 2 2 = κ Tr(Σ) 2
and E X (1) , X (2) 2 = Tr(Σ 2 ). We observe that

E Aθ 2 = E 1 q 2 (q -1) 2 1≤j<k≤q 1≤s<t≤q
θ, X (j) -X (k) X (j) -X (k) , X (s) -X (t) X (s) -X (t) , θ and E θ, X (1) -X (2) X (1) -X (2) , X (3) -X (4) X (3) -X (4) , θ = 4E θ, X (1) X (1) , X (2) X (2) , θ = 4 Σθ 2 ≤ 4 Σ ∞ N (θ), E θ, X (1) -X (2) X (1) -X (2) , X (1) -X (3) X (1) -X (3) , θ = E θ, X (1) 2 X (1) 2 + 3E θ, X (1) X (1) , X (2) X (2) , θ

≤ κ Tr(Σ)N (θ) + 3 Σ ∞ N (θ),
E θ, X (1) -X (2) X (1) -X (2) , X (1) -X (2) X (1) -X (2) , θ (1) X (1) , X (2) X (2) , θ

= 2E θ, X (1) 2 X (1) 2 + 2E θ, X (1) 2 E X (1) 2 + 4E θ, X
≤ 2(κ + 1) Tr(Σ)N (θ) + 4 Σ ∞ N (θ)
Therefore,

E Aθ 2 ≤ 1 - q -2 q(q -1) Σ ∞ N (θ) + 1 q κ + 1 q -1 Tr(Σ)N (θ).
which proves the first inequality. In the same way,

E Tr(A 2 ) = E 1 q 2 (q -1) 2 1≤j<k≤q 1≤s<t≤q X (j) -X (k) , X (s) -X (t) 2
, and E X (1) -X (2) , X (3) -X (4) 2 = 4E X (1) , X (2) 2 = 4 Tr Σ 2 , E X (1) -X (2) , X (1) -

X (3) 2 = E X (1) 4 + 3E X (1) , X (2) 2 ≤ κ Tr(Σ) 2 + 3 Tr Σ 2 , E X (1) -X (2) , X (1) -X (2) 2 = 2E X (1) 4 + 2E X (1) 2 2 + 4E X (1) , X (2) 2 ≤ 2(κ + 1) Tr(Σ) 2 + 4 Tr Σ 2 ,
which concludes the proof.

Using these tighter bounds, we can improve ζ to

ζ(t) = 2τ q (κ) (2 + 3c) 1 -q-2 q(q-1) Tr(Σ 2 ) + 1 q κ + 1 q-1 Tr(Σ) 2 4(2 + c) 1 -q-2 q(q-1) Σ ∞ + 1 q κ + 1 q-1 Tr(Σ) t + log K/ × cosh(a/4) + 2(2 + c) 1 -q-2 q(q-1) q Σ ∞ + κ + 1 q-1 Tr(Σ) t cosh(a/2).
Therefore, in the case when, q Σ ∞ ≤ Tr(Σ), we can take

ζ(t) = 2τ q (κ) (2 + 3c) Tr(Σ) 4(2 + c)t + log(K/ ) cosh(a/4) + 2(2 + c) κ + 1 + 2 q(q-1) Tr(Σ) t cosh(a/2).
If we compare this results with the bound obtained in Proposition 1.19 on page 38 for the Gram matrix estimator, we see that the first appearance of κ in the definition of ζ has been replaced with

τ q (κ) + 1 = κ + 2 q -1 ,
and that the second appearance of κ has been replaced with

κ + 1 + 2 q(q -1)
.

Thus, when Σ ∞ ≤ Tr(Σ)/2, that is not a very strong hypothesis, we can take at least q = 2, and obtain an improved bound for the estimation of Σ that is not much larger than the bound for the estimation of the centered Gram matrix, that requires the knowledge of E(X), since the difference between the two bounds is just a matter of replacing κ with κ+2.

Empirical results

In this section we present some empirical results that show the performance of our robust estimator Q.

We use in these experiments a simplified construction that does not lead exactly to the estimator Q, for which we have proved theoretical results in the previous sections, but should nonetheless exhibit the same kind of behaviour.

We use an iterative scheme based on the polarization formula to estimate the coefficients of the Gram matrix in an orthonormal basis of eigenvectors and then update this basis iteratively to a basis of eigenvectors of the current estimate.

Let X 1 , . . . , X n ∈ R d be a sample drawn according to the probability distribution P and let λ > 0. Let p ∈ R n and define S(p, λ) as the solution of

n i=1 ψ λ S(p, λ) -1 p 2 i -1 = 0.
In practice we compute S(p, λ) using the Newton algorithm. We observe that, when p i = θ, X i and λ is suitably chosen, S(p, λ) is an approximation of the estimator N (θ) of the quadratic form N (θ).

Define S(p) as the solution obtained when the parameter λ is set to

λ = m 1 v 2 n log( -1 ) 1 - 2 n log( -1 )
where

m = 1 n n i=1 p 2 i , v = 1 n -1 n i=1 p 2 i -m 2 and = 0.1.
According to [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF], this value of the scale parameter λ should be close to optimal for the estimation of a single expectation from an empirical sample distribution.

Let λ1 ≥ • • • ≥ λd ≥ 0 be the eigenvalues of the empirical Gram matrix Ḡ, that will be our starting point, and let u 1 , . . . , u d be a corresponding orthonormal basis of eigenvectors. We decompose the empirical Gram matrix as

Ḡ = U DU
where U is the orthogonal matrix whose columns are the eigenvectors of Ḡ and D is the diagonal matrix D = diag( λ1 , . . . , λd ). We observe that, by the polarization formula,

u i Gu j = 1 4 [N (u i + u j ) -N (u i -u j )] , i, j = 1, . . . , d,
where N (u i + σu j ), with σ ∈ {+1, -1}, is approximated by

S u i + σu j , X 2 , 1 ≤ ≤ n .
Taking notation, for any n × d matrix W , we define C(W ) as the d × d matrix of entries

C(W ) i,j = 1 4 S W ,i + W ,j 2 | 1 ≤ ≤ n -S W ,i -W ,j 2 | 1 ≤ ≤ n .
Let Y be the matrix whose i-th row is the vector X i , so that

(Y U ) ,k = u k , X , 1 ≤ ≤ n, 1 ≤ k ≤ d.
We update the Gram matrix estimate to

Q 0 = U C(Y U )U .
Then we iterate the update scheme, decomposing Q 0 as

Q 0 = O 0 D 0 O 0 ,
where O 0 O 0 = O 0 O 0 = I and D 0 is a diagonal matrix and computing

Q 1 = O 0 C(Y O 0 )O 0 .
The inductive update step is more generally the following. Assuming we have constructed Q k , we decompose it as

Q k = O k D k O k where O k O k = O k O k = I
and D k is a diagonal matrix and we define the new updated estimator of G as

Q k+1 = O k C(Y O k )O k .
We stop this iterative estimation scheme when Q k -Q k-1 F falls under a given threshold.

In the following numerical experiment we more simply performed four updates. We take as our robust estimator of G the last update Q k . We now present an example of the performance of this estimator, for some i.i.d. sample of size n = 100 in R 10 drawn according to the Gaussian mixture distribution P = (1 -α) N (0, M 1 ) + α N (0, 16 I),

where α = 0.05 and

M 1 =         2 1 0 1 1 0 0.01 0 . . . 0 0.01        
.

The Gram matrix of P is equal to

G = (1 -α)M 1 + 16 α I =         2
        .
This example illustrates a favorable situation where the performance of the robust estimator is particularly striking when compared to the empirical Gram matrix. As it can be seen on the expression of the sample distribution as well as on the configuration plots below, this is a situation of intermittent high variance : the sample is a mixture of a rare high variance signal and a frequent low variance more structured signal.

We tested the algorithm on 500 different samples drawn according to the Gaussian mixture distribution defined above. Random sample configurations are presented in fig. 1.1.

Figure 1.2 shows that the robust estimator Q significantly improves the error in terms of square of the Frobenius norm when compared to the empirical estimator Ḡ. The red solid line represents the empirical quantile function of the errors of the robust estimator, whereas the blue dotted line represents the quantiles of Ḡ -G 2 F . This quantile function is obtained by sorting the 500 empirical errors in increasing order.

The mean of the square distances Q -G 2 F on 500 trials is 5.6 ± 0.4, where the indicated mean estimator and confidence interval is the non-asymptotic confidence interval given by Proposition 2.4 of [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF] at confidence level 0.99. In the case of the empirical estimator, the mean is 15.5 ± 2. The empirical standard deviations are respectively 2 and 10. So we see that in this case the robust estimator reliably decreases the error by a factor larger than 2 and also produces errors with a much smaller standard deviation from sample to sample. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -5 0 5 10 -5 0 5 n1 = 6, n2 = 94 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -8 -6 -4 -2 0 2

= 6, n2 = 94 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -6 -4 -2 0 2 4

-2 0 2 n1 = 6, n2 = 94 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -6 -4 -2 0 where C is an absolute constant, then

Ḡ -G ∞ ≤ δ G ∞ .
Typically the result is used with m = O(d G ∞ ), so that the conclusion of Proposition 2.1 holds with sample size n ≥ Cδ -2 d log(d/ ).

This means that we need to choose the sample size n = O(d log(d)) to approximate the Gram matrix by the empirical Gram matrix.

A different approach is used in Tropp [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF], where a bound on the operator norm of the approximation error is deduced from the matrix Bernstein inequality. The author observes that the empirical Gram matrix is an unbiased estimator of the Gram matrix, i.e. E[ Ḡ] = G, which can be expressed as a sum of independent random matrices. Thus, a strategy to bound the approximation error is to look whether the random matrix Ḡ deviates from its mean. It is known that in the case of a random variable Z expressed as a sum of independent random variables the Bernstein inequality provides a bound on the probability that Z takes values far from its mean. We recall the Bernstein inequality.

Proposition 2.2. (Bernstein inequality)

Let S 1 , . . . , S n ∈ R be independent random variables such that, for any k = 1, . . . , n,

S k -E(S k ) ≤ R. Let Z = n k=1 S k and σ 2 = E Z -E(Z) 2 .
Then, for t ≥ 0,

P Z -E(Z) ≥ t ≤ 2 exp - t 2 /2 σ 2 + Rt/3 .
What is proved in [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF] is that the scalar Bernstein inequality extends directly to matrices, that is, it is possible to obtain an exponential concentration inequality for the operator norm of a sum of independent random matrices.

Proposition 2.3. (Matrix Bernstein inequality, [34]). Let

S 1 , . . . , S n ∈ M d 1 ,d 2 (R) be independent random matrices of size d 1 × d 2 , such that, for any k = 1, . . . , n, S k -E(S k ) ∞ ≤ R. Let Z = n k=1 S k and σ 2 = max E (Z -E(Z))(Z -E(Z)) ∞ , E (Z -E(Z)) (Z -E(Z)) ∞ .
Then, for t ≥ 0,

P [ Z -E(Z) ∞ ≥ t] ≤ (d 1 + d 2 ) exp - t 2 /2 σ 2 + Rt/3 .
We now describe how to apply the matrix Bernstein inequality to study of the Gram matrix.

Let X ∈ R d be a random vector such that X 2 ≤ B, and let X 1 , . . . , X n ∈ R d be an i.i.d. sample drawn according to P. We write the approximation error as

Ḡ -G = n k=1 S k where S k = 1 n X k X k -G . We have E[S k ] = 0 and S k ∞ ≤ 1 n X k 2 + E[ X 2 ] ≤ 2B n .
It remains to bound the variance σ 2 . We observe that

σ 2 = E n k=1 S k 2 ∞ = n k=1 E S 2 k ∞
where the last identity depends on the fact that S k are independent zero-mean matrices. Since

E S 2 k = 1 n 2 E X k X k -G 2 = 1 n 2 E X k 2 X k X k -X k X k G -GX k X k + G 2 = 1 n 2 E X k 2 X k X k -G 2 1 n 2 BG,
where H M means that M -H is positive semi-definite, we conclude that

σ 2 ≤ B n G ∞ .
Applying the matrix Bernstein inequality with d 1 = d 2 = d and R = 2B n , we obtain, with probability at least 1 -2 ,

Ḡ -G ∞ ≤ 1 3   2B n log(d/ ) + 4B 2 n 2 log(d/ ) 2 + 18 B G ∞ n log(d/ )   . (2.1)
We observe that the operator norm

Ḡ -G ∞ = sup θ∈S d |θ Ḡθ -θ Gθ|
rewrites, according to the notation of chapter 1, as

Ḡ -G ∞ = sup θ∈S d | N (θ) -N (θ)|,
where S d denotes the unit sphere of R d .

In this chapter we present, with the help of the robust estimator N introduced in chapter 1, a uniform bound on the approximation error |N (θ) -N (θ)| which does not depend explicitly on the dimension d ( Proposition 2.7).

Estimate of the Gram matrix via the empirical Gram matrix

In this section we aim at estimating the quadratic form N (θ) = θ Gθ by the empirical estimator

N (θ) = 1 n n i=1 θ, X i 2 , θ ∈ R d ,
from an i.i.d. sample X 1 , . . . , X n ∈ R d drawn according to P. We will use the robust estimator constructed in chapter 1 to provide a uniform bound on the approximation error. We briefly recall its definition.

For θ ∈ R d and λ > 0, we consider the empirical criterion

r λ (θ) = 1 n n i=1 ψ θ, X i 2 -λ ,
where the influence function ψ is introduced in equation (1.2) on page 16, and we define α(θ) = sup {α ∈ R + | r λ (αθ) ≤ 0} . We consider the family of estimators

N λ (θ) = λ α(θ) 2 , θ ∈ R d . Let Λ ⊂ R + \ {0}
2 be the finite set defined in equation (1.27) on page 34. We observe that, according to the definition of α(θ), for any threshold σ ∈ R + ,

1 n n i=1 ψ λ max{ N λ (θ), σ} -1 θ, X i 2 -1 ≤ r λ 1/2 N λ (θ) -1/2 θ = r λ α(θ) θ ≤ 0,
where we have used the fact that the function ψ is non-decreasing. Moreover r λ α(θ) θ = 0 as soon as α(θ) < +∞ and this holds true, according to Proposition 1.12 on page 28, with probability at least 1 -, for any θ ∈ S d and any (λ, β) ∈ Λ such that Φ + N (θ) > 0. Indeed, by Proposition 1.12, with probability at least 1 -, for any θ ∈ S d , any (λ, β) ∈ Λ,

N λ (θ) ≥ Φ + N (θ) . Let us put R = max i=1,...,n X i (2.2)
and

τ λ (t) = λ 2 R 4 3 max{t, σ} 2 , t ∈ R + . Proposition 2.4. Let σ ∈ R + be some threshold. With probability at least 1 -, for any θ ∈ S d , any (λ, β) ∈ Λ, such that Φ + N (θ) > 0, N (θ) max{ N λ (θ), σ} ≤ 1 -τ λ N λ (θ) -1 + ,
with the convention that 1 0 = +∞. Moreover, with probability at least 1 -, for any θ ∈ S d , any (λ, β) ∈ Λ, such that Φ + N (θ) > 0, N (θ)

N λ (θ) ≥ 1 - λ 2 3 .
Proof. According to Proposition 1.12, with probability at least 1 -, for any θ ∈ S d , any (λ, β) ∈ Λ, such that Φ + N (θ) > 0

1 n n i=1 ψ λ N λ (θ) -1 θ, X i 2 -1 = r λ α(θ) θ = 0. Defining g(z) = z -ψ(z), we get N (θ) max{ N λ (θ), σ} -1 ≤ 1 nλ n i=1 g λ θ, X i 2 max{ N λ (θ), σ} -1 -1 . (2.3)
In the same way, with probability at least 1 -, for any θ ∈ S d , any (λ, β) ∈ Λ such that Φ + N (θ) > 0, we obtain

1 - N (θ) N λ (θ) ≤ 1 nλ n i=1 g λ 1 -θ, X i 2 N λ (θ) -1 . (2.4)
We remark that the derivative of g is

g (z) = 1 -ψ (z) =                  1 if z / ∈ [-1, 1] z 2 2 1 + z + z 2 2 if z ∈ [-1, 0] z 2 2 1 -z + z 2 2 if z ∈ [0, 1],
showing that 0 ≤ g (z) ≤ z 2 , and therefore that g is a non-decreasing function satisfying

g(z) ≤ 1 3 z 3 + , ( 2.5)
where z + = max{z, 0}. Indeed, it is sufficient to observe that g(z) ≤ 0 if z ≤ 0, while, for z ≥ 0,

g(z) = z 0 g (s)ds ≤ z 0 s 2 ds = 1 3 z 3 .
Applying equation (2.5) to equation (2.3) we obtain

N (θ) max{ N λ (θ), σ} -1 ≤ λ 2 3n n i=1 θ, X i 2 max{ N λ (θ), σ} -1 -1 3 + ≤ λ 2 3n max{ N λ (θ), σ} 3 n i=1 θ, X i 6 ,
where we have used the fact (z 2 -1) + ≤ z 2 .

Since, by the Cauchy-Schwarz inequality, θ, X i 2 ≤ θ 2 R 2 , we get

N (θ) max{ N λ (θ), σ} -1 ≤ λ 2 3n max{ N λ (θ), σ} 3 θ 4 R 4 n i=1 θ, X i 2 = λ 2 3 × θ 4 R 4 max{ N λ (θ), σ} 2 × N (θ) max{ N λ (θ), σ} .
Similarly, since g in non-decreasing, we obtain that, with probability at least 1 -, for any

θ ∈ S d , any (λ, β) ∈ Λ such that Φ + N (θ) > 0, 1 - N (θ) N λ (θ) ≤ 1 nλ n i=1 g(λ) ≤ λ 2 3 ,
where the last inequality follows from equation (2.5).

The above result relates the behavior of empirical estimator N to that of N λ and the accuracy of the approximation depends on R, defined in equation (2.2). At the end of the section we mention some assumptions under which it is possible to give a non-random bound on R.

Before stating the bound on the approximation error we introduce the following result.

Proposition 2.5. With probability at least 1 -2 , for any θ ∈ S d , any

(λ, β) ∈ Λ, any σ > 0, max{ N (θ), σ} ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ) -1 + max N (θ), σ ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ) 1 -τ λ (σ) + -1 + N (θ) ≥ 1 - λ 2 3 + Φ + N (θ)
where Φ + and Φ -are defined in Proposition 1.12 on page 28.

Proof. We consider the threshold

σ = Φ -1 -max{N (θ), σ}) ≥ max{N (θ), σ},
where we have used the fact that, by definition, Φ -(t) -1 ≥ t, for any t ∈ R + . We assume that we are in the intersection of the two events of Proposition 1.12, which holds true with probability at least 1 -2 , so that σ ≥ max{N (θ), σ, N λ (θ)}.

(2.6) By Proposition 2.4, choosing as threshold max{σ, σ }, we get

N (θ) max{ N λ (θ), σ, σ } ≤ 1 -τ λ max{ N λ (θ), σ }) -1 + ,
(where τ λ is still defined with respect to σ), so that, according to equation (2.6),

N (θ) ≤ σ 1 -τ λ σ -1 + . (2.7)
As a consequence, recalling the definition of σ , we have

N (θ) ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ) -1 + .
Thus, observing that

σ ≤ Φ -1 -(σ) ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ) -1 + ,
we obtain the first inequality of the proposition. To prove the second inequality, we use equation (2.7) once to see that

σ ≥ N (θ) 1 -τ λ (σ ) + ≥ N (θ) 1 -τ λ (σ) + ,
and we use it again to get

N (θ) ≤ Φ -1 -max{N (θ), σ} 1 -τ λ (σ ) -1 + ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ)[1 -τ λ (σ)] + -1 + .
To complete the proof of the second inequality, it is enough to remark that

σ ≤ Φ -1 -(σ) ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ) 1 -τ λ (σ) + -1 + .
To prove the last inequality, it is sufficient to remark that N λ (θ) ≥ Φ + N (θ) by Proposition 1.12 and hence, when

Φ + N (θ) > 0, N (θ) ≥ 1 - λ 2 3 + Φ + N (θ) .
On the other hand, when Φ + N (θ) = 0, this inequality is also obviously satisfied.

According to the previous result we present a bound on the approximation error which uses the notation introduced in chapter 1.

Corollary 2.6. With probability at least 1-2 , for any θ ∈ S d , any (λ, β) ∈ Λ, any σ > 0, max{ N (θ), σ} max{N (θ), σ}

-1 ≤ B λ,β N (θ) + τ λ N (θ) 1 -τ λ N (θ) + 1 -B λ,β N (θ) + , 1 - max{ N (θ), σ} max{N (θ), σ} ≤ B λ,β N (θ) + λ 2 3 ,
where B λ,β is defined in equation (1.20) on page 31 and B λ,β on page 33. In particular

max{ N (θ), σ} max{N (θ), σ} -1 ≤ B λ,β N (θ) + τ λ N (θ) 1 -τ λ N (θ) + 1 -B λ,β N (θ) + . (2.8)
Proof. We observe that, by Proposition 2.5,

max{ N (θ), σ} ≤ Φ -1 -max{N (θ), σ} 1 -τ λ N (θ) -1 + ≤ max{N (θ), σ} 1 -τ λ N (θ) + 1 -B λ,β N (θ) + , which implies max{ N (θ), σ} max{N (θ), σ} -1 ≤ B λ,β N (θ) 1 -B λ,β N (θ) + + τ λ N (θ) 1 -τ λ N (θ) + 1 -B λ,β N (θ) + .
Applying Lemma 1.15 on page 33 we obtain the first inequality.

To prove the second inequality we observe that, by Proposition 2.5, max{ N (θ), σ} ≥ 1 -

λ 2 3 + Φ + max{N (θ), σ} = 1 - λ 2 3 + max{N (θ), σ} 1 + B λ,β N (θ) -1 ,
where we have used the fact that Φ + max{z, σ} = max{z, σ} 1 + B λ,β (z)

-1
as shown in equation (1.23) on page 31. Thus we conclude that

1 - max{ N (θ), σ} max{N (θ), σ} ≤ B λ,β N (θ) + λ 2 /3 1 + B λ,β N (θ) ≤ B λ,β N (θ) + λ 2 3 .
We observe that the fact of introducing the robust estimator N λ in the computation of the approximation error has led to provide bounds which do not depend explicitly on the dimension d.

According to the results introduced in chapter 1 we conclude the section by presenting a more explicit (dimension-free) bound on the approximation error.

We recall some notation. Let a > 0 and let

K = 1 + a -1 log n 72(2 + c)κ 1/2 where κ = sup θ∈R d E( θ,X 2 )>0 E θ,X 4 E θ,X 2 2 and c = 15 8 log(2)( √ 2-1) exp 1+2 √ 2 2
. We put

B * (t) =      n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ}) 6 + (κ -1) -1 ζ(max{t, σ}) ≤ √ n +∞ otherwise
where

ζ(t) = 2(κ -1) (2 + 3c) E X 4 1/2 4(2 + c)κ 1/2 t + log(K/ ) cosh(a/4) + 2(2 + c)κ 1/2 E X 4 1/2 t cosh(a/2).
The following proposition holds true. Let us assume that, for some exponent p ≥ 1 and some positive constants α and η,

Proposition 2.7. Consider any threshold

σ ∈ R + such that σ ≤ E X 4 1/2 . Define τ * (t) = λ * (t) 2 exp(a/2)R 4 3 max{t, σ} 2 , t ∈ R + ,
E exp α 2 X 2/p Tr(G) 1/p -1 -η 2/p ≤ 1.
In this case, with probability at least 1 -,

R ≤ Tr(G) 1/2 1 + η 2/p + 2α -1 log n/ p/2 , ( 2.9)
where we recall that

Tr(G) = E X 2 .
To give a point of comparison, in the centered Gaussian case where X ∼ N (0, G) is a Gaussian vector, we have, for any α

∈ [0, λ -1 1 Tr(G)[, E exp α 2 X 2 Tr(G) + 1 α d i=1 log 1 - αλ i Tr(G) = 1,
where λ 1 ≥ • • • ≥ λ d are the eigenvalues of G. Therefore, with probability at least 1 -,

R ≤ Tr(G) 1/2 - 1 α d i=1 log 1 - αλ i Tr(G) + 2 log(n/ ) α 1/2 .
Moreover we observe that

lim α→0 + - 1 α d i=1 log 1 - αλ i Tr(G) = 1.
In order to replace equation (2.9) with some polynomial assumptions we need to replace R by

R = 1 n n i=1 X i 6 1/6
. Indeed, by the Bienaymé Chebyshev inequality, we get that, with probability at least 1-,

R ≤ E X 6 + E X 12 n 1/2 1/6 ≤ 1 + (n ) -1/2 1/6 E X 12 1/12
and hence, with probability at least 1

-n -1 , R ≤ 2 1/6 E X 12 1/12 .
Also in the case where we consider this new quantity R we obtain a bound of the form of Proposition 2.7. Indeed we observe that another way to take advantage of equation ( 2.3) is to write

N (θ) max{ N λ (θ), σ} -1 ≤ λ 2 θ 6 3 max{ N λ (θ), σ} 3 1 n n i=1 X i 6 .
Thus, putting

ζ λ (t) = λ 2 R 6 3 max{t, σ} 3 , t ∈ R + , we get that, for any θ ∈ S d , N (θ) max{ N λ (θ), σ} ≤ 1 + ζ λ N λ (θ) .
The same reasoning used to prove Proposition 2.5 shows that, with probability at least 1 -, for any θ ∈ S d , any (λ,

β) ∈ Λ, any σ > 0, max{ N (θ), σ} ≤ Φ -1 -max{N (θ), σ} 1 + ζ λ N (θ) .
As a consequence, with probability at least 1 -2 , for any θ ∈ S d , any (λ,

β) ∈ Λ, max{ N (θ), σ} max{N (θ), σ} -1 ≤ B λ,β N (θ) + ζ λ N (θ) 1 -B λ,β N (θ) + .
We conclude by stating the analogous of Proposition 2.7.

Proposition 2.8. Let 0 < σ ≤ E X 4 1/2 and let us put

ζ * (t) = λ * (t) 2 exp(a/2) R 6 3 max{t, σ} 3 , t ∈ R + .
With probability at least 1 -2 , for any

θ ∈ S d , max{ N (θ), σ} max{N (θ), σ} -1 ≤ B * N (θ) + ζ * N (θ) 1 -B * N (θ) + .
Chapter 3

Principal Component Analysis

We use the quadratic estimator introduced in chapter 1 to construct robust estimators of the eigenvalues of the Gram matrix. Based on this result we propose a new approach that qualifies the stability of principal component analysis independently of the dimension of the ambient space.

Introduction

Principal Component Analysis (PCA) is a classical tool for dimensionality reduction. The basic idea of PCA is to reduce the dimensionality of a dataset by projecting it into the space spanned by the directions of maximal variance, that are called its principal components. Since this set of directions lies in the space generated by the eigenvectors associated with the largest eigenvalues of the covariance matrix of the sample, the dimensionality reduction is achieved by projecting the dataset into the space spanned by these eigenvectors, which in the following we call largest eigenvectors.

Given X ∈ R d a random vector distributed according to the unknown probability distribution P ∈ M 1 + (R d ), the goal is to estimate the eigenvalues and eigenvectors of the covariance matrix of X

Σ = E X -E(X) X -E(X)
from an i.i.d. sample X 1 , . . . , X n ∈ R d drawn according to P. We have already observed that in the case where the random vector X is centered the covariance matrix Σ is the Gram matrix G = E XX .

In this chapter we will consider the case of the Gram matrix but similar results can be deduced for the covariance matrix.

As in chapter 1 we introduce the quadratic form

N (θ) = θ Gθ, θ ∈ R d ,
and we observe that, denoting by p 1 , . . . , p d an orthonormal basis of eigenvectors of G, the i-th eigenvalue of the Gram matrix is λ i = N (p i ).

From now we denote by λ 1 ≥ • • • ≥ λ d the eigenvalues of G and we assume that the eigenvectors p 1 , . . . , p d are ranked according to the decreasing order of their eigenvalues.

Let G = Q + denote the positive part of the symmetric matrix Q defined in equation (1.37) on page 40. Let λ 1 ≥ • • • ≥ λ d be its eigenvalues and q 1 , . . . , q d the corresponding orthonormal basis of eigenvectors.

In this chapter we propose some robust versions of PCA. We first show in section 3.2 that each eigenvalue λ i of G is a robust estimator of the corresponding eigenvalue of the Gram matrix. As a consequence, the orthogonal projector on the largest eigenvectors of G can be estimated by the projector on the largest eigenvectors of G. The behavior of this estimator is related to the size of the gap in the spectrum of the Gram matrix and in order to have a good approximation we need to have a large eigengap as shown in section 3.3. To avoid the assumption of a large gap in the spectrum of G we propose in section 3.4 a robust version of PCA which consists in performing a smooth cut-off of the spectrum of the Gram matrix via a Lipschitz function. We provide bounds on the approximation error, in terms of the operator norm ( Proposition 3.8) and of the Frobenius norm ( Proposition 3.9), that replace the size of the eigengap by the inverse of the Lipschitz constant.

Estimate of the eigenvalues

In this section we prove that each eigenvalue of G is a robust estimator of the corresponding eigenvalue of the Gram matrix.

We first recall some notation. Let a > 0 and let

K = 1 + a -1 log n 72(2 + c)κ 1/2 where c = 15 8 log(2)( √ 2-1) exp 1+2 √ 2 2
and κ = sup

θ∈R d E( θ,X 2 )>0 E θ, X 4 E θ, X 2 2 .
Let s 2 4 = E X 4 1/2 and let σ ∈]0, s 2 4 ] be a threshold. We put

B * (t) =      n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ}) 6 + (κ -1) -1 ζ(max{t, σ}) ≤ √ n +∞ otherwise
where

ζ(t) = 2(κ -1) (2 + 3c) s 2 4 4(2 + c)κ 1/2 t + log(K/ ) cosh(a/4) + 2(2 + c)κ 1/2 s 2 4 t cosh(a/2). (3.1)
The following result holds.

Proposition 3.1. Let us assume that 8ζ(σ) ≤ √ n, σ ≤ s 2 4 and that κ ≥ 3/2. With probability at least 1 -2 , for any i = 1, . . . , d, the two following inequalities hold together

max{λ i , σ} -max{ λ i , σ} ≤ 2 max{λ i , σ}B * (λ i ) + 5δ G F , max{λ i , σ} -max{ λ i , σ} ≤ 2 max{ λ i , σ}B * min λ i , s 2 4 + 5δ G F .
Proof. We observe that, for any i ∈ {1, . . . , d}, the vector space span{q 1 , . . . , q i-1 } ⊥ ∩ span{p 1 , . . . , p i } ⊂ R d is of dimension at least 1, so that the set

V i = θ ∈ S d | θ ∈ span{q 1 , . . . , q i-1 } ⊥ ∩ span{p 1 , . . . , p i } ⊂ R d is non-empty. Indeed, putting A = span{q 1 , . . . , q i-1 } ⊥ and B = span{p 1 , . . . , p i }, we see that dim(A ∩ B) = dim(A) + dim(B) -dim(A + B) ≥ 1, since dim(A + B) ≤ dim(R d ) = d and dim(A) + dim(B) = d + 1.
Hence, there exists θ i ∈ V i and for such a θ i , we have

N (θ i ) ≥ λ i . It follows that max{λ i , σ} ≤ sup {max{N (θ), σ} | θ ∈ V i } ≤ sup max{N (θ), σ} | θ ∈ S d , θ ∈ span{q 1 , . . . , q i-1 } ⊥ .
Therefore, according to Lemma 1.25 on page 43,

max{λ i , σ} 1 -2B * (λ i ) ≤ sup max{θ Qθ, σ} | θ ∈ S d ∩ span{q 1 , . . . , q i-1 } ⊥ + 5δ G F ≤ max{ λ i , σ} + 5δ G F .
In the same way,

max{ λ i , σ} ≤ sup max{N (θ), σ} 1 + 2B * N (θ) | θ ∈ S d ∪ span{p 1 , . . . , p i-1 } ⊥ + 5δ G F ≤ max{λ i , σ} 1 + 2B * (λ i ) + 5δ G F , max{ λ i , σ} 1 -2B * min{ λ i , s 2 4 } ≤ sup max{N (θ), σ} | θ ∈ S d ∪ span{p 1 , . . . , p i-1 } ⊥ + 5δ G F ≤ max{λ i , σ} + 5δ G F , max{λ i , σ} ≤ sup max{θ Qθ, σ} 1 + 2B * min{θ Qθ, s 2 4 } | θ ∈ S d ∪ span{q 1 , . . . , q i-1 } ⊥ + 5δ G F , ≤ max{ λ i , σ} 1 + 2B * min{ λ i , s 2 4 } + 5δ G F .
In all these inequalities we have used the fact that

t → max{t, σ} 1 -2B * (min{t, s 2 4 }) t → max{t, σ} 1 + 2B * (min{t, s 2 4 })
are non-decreasing (the latter according to Lemma 1.22 on page 41) and that λ i ≤ s 2 4 . This proves the proposition for the eigenvalues of Q, and therefore also for their positive parts, that are the eigenvalues of Q + .

As a consequence, using the fact that

|λ i -λ i | ≤ max{λ i , σ} -max{ λ i , σ} + σ,
we obtain the following result.

Corollary 3.2. Under the same assumptions as in Proposition 3.1, with probability at least 1 -2 , for any

i = 1, . . . , d, |λ i -λ i | ≤ 2 max{λ i , σ}B * λ i + 5δ G F + σ, |λ i -λ i | ≤ 2 max{ λ i , σ}B * min λ i , s 2 4 + 5δ G F + σ.
We present two more bounds on the estimation error that depend on the largest eigenvalues λ 1 and λ 1 respectively. Corollary 3.3. Under the same assumptions as in Proposition 3.1, with probability at least 1 -2 , for any i = 1, . . . , d,

|λ i -λ i | ≤ 2 max{λ 1 , σ}B * λ 1 + 5δ G F + σ, |λ i -λ i | ≤ 2 max{ λ 1 , σ}B * min λ 1 , s 2 4 + 5δ G F + σ.
Proof. The proof follows from Lemma 1.22 on page 41.

In order to simplify notation, we define

B(t) = 2 max{t, σ}B * min{t, s 2 4 } + 7δ G F + σ. ( 3

Standard Principal Component Analysis

A method to determine the number of relevant components is based on the difference in magnitude between successive eigenvalues. In this section we study the projection on the r largest eigenvectors p 1 , . . . , p r of the Gram matrix, assuming that there is a gap in the spectrum of the Gram matrix, meaning that for some positive constant δ, λ r -λ r+1 ≥ δ.

(3.3)

We denote by Π r the orthogonal projector on the r largest eigenvectors p 1 , . . . , p r of G and similarly by Π r the orthogonal projector on the r largest eigenvectors q 1 , . . . , q r of its estimate G.

Our goal is to provide a bound on the approximation Π r -Π r ∞ that does not depend explicitly on the dimension d of the ambient space.

Since Π r and Π r have the same rank, we can write

Π r -Π r ∞ = sup θ∈S d θ∈Im( Πr) Π r θ -Π r θ as shown in Lemma B.5 in appendix B. Moreover, for any θ ∈ Im( Π r ) ∩ S d , we observe that Π r θ -Π r θ 2 = Π r θ -θ 2 = r i=1 θ, p i p i - d i=1 θ, p i p i 2 = d i=r+1 θ, p i 2 .
Since θ ∈ Im( Π r ), it can be written as θ = r k=1 θ, q k q k with r k=1 θ, q k 2 = 1, so that

Π r θ -Π r θ 2 = d i=r+1 r k=1 θ, q k q k , p i 2 .
Hence, by the Cauchy-Schwarz inequality, we get

Π r θ -Π r θ 2 ≤ d i=r+1 r k=1 θ, q k 2 r k=1 q k , p i 2 (3.4) = r k=1 d i=r+1 q k , p i 2 . ( 3
(λ i -λ k ) 2 q k , p i 2 ≤ 2B (λ 1 ) 2 (3.6) d i=1 λ i -λ k 2 q k , p i 2 ≤ B (λ 1 ) 2 , (3.7)
where B is defined in equation (3.2).

Proof. We observe that,

G -G ∞ = max sup θ∈S d θ G -G θ, sup θ∈S d θ G -G θ = sup θ∈S d N (θ) -θ Gθ .
Thus, by Corollary 1.28 on page 44, with probability at least 1 -2 ,

sup θ∈S d Gθ -Gθ ≤ B (λ 1 ) .
To prove equation (3.7) it is sufficient to observe that, since

Gθ -Gθ = d i,j=1 (λ i -λ j ) θ, q j p i , q j p i choosing θ = q k , with k ∈ {1, . . . , d}, Gq k -Gq k 2 = d i=1 (λ i -λ k ) 2 q k , p i 2 .
On the other hand, to prove equation (3.6), we observe that

Gθ -Gθ = d i=1 λ i θ, p i p i - d i=1 λ i θ, q i q i = d i=1 λ i ( θ, p i p i -θ, q i q i ) - d i=1 λ i -λ i θ, q i q i ≥ d i=1 λ i ( θ, p i p i -θ, q i q i ) - d i=1 λ i -λ i θ, q i q i
where, by Corollary 3.3,

d i=1 λ i -λ i θ, q i q i 2 = d i=1 λ i -λ i 2 θ, q i 2 ≤ B (λ 1 ) 2 .
Choosing again θ = q k , for k ∈ {1, . . . , d}, we get

d i=1 λ i ( q k , p i p i -q k , q i q i ) 2 = d i,j=1 (λ i -λ j ) q k , q j q j , p i p i 2 = d i=1 (λ i -λ k ) q k , p i p i 2 = d i=1 (λ i -λ k ) 2 q k , p i 2 ,
which concludes the proof.

We now apply Lemma 3.4 to our problem. The following proposition holds.

Proposition 3.5. With probability at least

1 -2 , Π r -Π r ∞ ≤ √ 2r λ r -λ r+1 B (λ 1 )
where B is defined in equation (3.2) and λ 1 is the largest eigenvalue of the Gram matrix.

Proof. We observe that, for any k ∈ {1, . . . , r},

d i=r+1 (λ k -λ i ) 2 q k , p i 2 ≥ d i=r+1 (λ r -λ i ) 2 q k , p i 2 ≥ (λ r -λ r+1 ) 2 d i=r+1 q k , p i 2 .
Then, by Lemma 3.4, with probability at least 1 -2 ,

(λ r -λ r+1 ) 2 d i=r+1 q k , p i 2 ≤ 2B (λ 1 ) 2 .
Applying the above inequality to equation (3.5) we conclude the proof.

We observe that the above proposition provides a bound on the approximation error Π r -Π r ∞ that does not depend on the dimension d of the ambient space, since B is dimension-free. However the result relates the quality of the approximation of the orthogonal projector Π r by the robust estimator Π r to the size of the spectral gap. In particular the larger the eigengap, the better the approximation is.

Robust Principal Component Analysis

In order to avoid the requirement of a large spectral gap, we replace the mere projection on the largest eigenvectors of the Gram matrix by a smooth cut-off of the spectrum of G via a Lipschitz function. In particular we consider a function which is one on the largest eigenvalues and is zero on the smallest ones. In the case there is a (sufficiently large) gap in the spectrum of G, the best Lipschitz constant is exactly the inverse size of the eigengap. In such a way we replace, in some sense, the size of the eigengap with the inverse of the Lipschitz constant.

Let f be a Lipschitz function with Lipschitz constant 1/L. We decompose the Gram matrix as

G = U DU where D = diag(λ 1 , . . . , λ d ) ∈ M d (R) is the diagonal matrix whose entries are the eigen- values of G and U ∈ M d (R) is the orthogonal matrix of eigenvectors of G. We define f (G) as f (G) = U f (D)U
where, for any i, j ∈ {1, . . . , d},

|f (λ i ) -f (λ j )| ≤ 1 L |λ i -λ j |.
We provide some results on the estimate of f (G), the image of the Gram matrix by the smooth cut-off f , in terms of the operator norm • ∞ and of the Frobenius norm • F . We recall that, given M ∈ M d (R) a symmetric matrix, the Frobenius norm of M is defined as

M 2 F = Tr(M M ) and that M ∞ ≤ M F .
Proposition 3.6. Let M, M ∈ M d (R) be two symmetric matrices. We denote by µ 1 , . . . , µ d the eigenvalues of M related to the orthonormal basis of eigenvectors p 1 , . . . , p d and by µ 1 , . . . , µ d the eigenvalues of M related to the orthonormal basis of eigenvectors q 1 , . . . , q d . We have

M -M 2 F = d i,k=1 (µ i -µ k ) 2 p i , q k 2 .
Proof. Since {p i } d i=1 is an orthonormal basis of eigenvectors of M and {µ i } d i=1 the corresponding eigenvalues, we can write M as

M = d i=1 µ i p i p i . Similarly, M = d i=1 µ i q i q i .
Our goal is to evaluate M -M F where, by definition,

M -M = d i=1 µ i p i p i - d k=1 µ k q k q k = d i,k=1 (µ i -µ k ) p i , q k q k p i .
We observe that M -M is a symmetric matrix and its Frobenius norm is

M -M 2 F = Tr((M -M ) (M -M )),
where

(M -M ) (M -M ) = d i,j,k=1 (µ i -µ k )(µ i -µ j ) p i , q k p i , q j q k q j .
Considering that Tr(q k q j ) = δ jk , we conclude that

M -M 2 F = d i,j,k=1 (µ i -µ k )(µ i -µ j ) p i , q k p i , q j δ jk = d i,k=1 (µ i -µ k ) 2 p i , q k 2 .
Corollary 3.7. Let f be a 1/L-Lipschitz function. We have

f (M ) -f (M ) F ≤ 1 L M -M F .
Proof. It is sufficient to observe that using twice Proposition 3.6 we obtain

f (M ) -f (M ) 2 F = d i,k=1 (f (µ i ) -f (µ k )) 2 p i , q k 2 ≤ 1 L 2 d i,k=1 (µ i -µ k ) 2 p i , q k 2 = 1 L 2 M -M 2 F .
We recall that we are more particularly interested in the Lipschitz function f that is one on the largest eigenvalues and zero on the smallest ones.

f (G) -f ( G) ∞ ≤ min r∈{1,...,d} L -1   B (λ 1 ) + 4rB (λ 1 ) 2 + 2 d i=r+1 λ 2 i   , where B is defined in equation (3.2) and λ 1 ≥ • • • ≥ λ d are the eigenvalues of G.
Proof. We observe that all the proof holds in the event of probability at least 1 -2 described in Proposition 3.1. Let H ∈ M d (R) be the matrix defined as

H = d k=1 λ k q k q k .
We observe that

f (G) -f ( G) ∞ ≤ f (G) -f (H) ∞ + f (H) -f ( G) ∞
and we look separately at the two terms. By definition of operator norm, we have

f (H) -f ( G) 2 ∞ = sup θ∈S d f (H)θ -f ( G)θ 2 = sup θ∈S d d k=1 (f (λ k ) -f ( λ k )) θ, q k q k 2 = sup θ∈S d d k=1 (f (λ k ) -f ( λ k )) 2 θ, q k 2 .
Since the function f is 1/L-Lipschitz, we get

f (H) -f ( G) 2 ∞ ≤ L -2 sup θ∈S d d k=1 (λ k -λ k ) 2 θ, q k 2
and then, applying Corollary 3.3, with probability at least 1 -2 , we obtain

f (H) -f ( G) 2 ∞ ≤ L -2 B (λ 1 ) 2 .
On the other hand, we have

f (G) -f (H) ∞ ≤ f (G) -f (H) F ≤ 1 L G -H F ,
as shown in Corollary 3.7. Hence, according to Proposition 3.6, we get

f (G) -f (H) 2 ∞ ≤ 1 L 2 d i,k=1 (λ i -λ k ) 2 p i , q k 2 where d i,k=1 (λ i -λ k ) 2 p i , q k 2 ≤   r i=1 d k=1 + d i=1 r k=1 + d i,k=r+1   (λ i -λ k ) 2 p i , q k 2 .
Since λ i ≥ 0, for any i ∈ {1, . . . , d}, we get

d i,k=r+1 (λ i -λ k ) 2 p i , q k 2 ≤ 2 d i=r+1 λ 2 i .
Moreover, by Lemma 3.4, we have

r k=1 d i=1 (λ i -λ k ) 2 p i , q k 2 ≤ 2rB (λ 1 ) 2
and since the same bound also holds for

r i=1 d k=1 (λ i -λ k ) 2 p i , q k 2 ,
we conclude the proof.

Slightly changing the definition of the estimator we present a bound for the approximation error in terms of the Frobenius norm. Instead of considering G = d i=1 λ i q i q i we consider the matrix G = d i=1 λ i q i q i with eigenvectors q 1 , . . . , q d and eigenvalues

λ i = λ i -B( λ i ) +
where we recall that B is defined in equation (3.2). We observe that, in the event of probability at least 1 -2 described in Corollary 3.3, for any i = 1, . . . , d,

λ i ≤ λ i .
We first present a result on the approximation error G -G F .

Proposition 3.9. (Frobenius norm) With probability at least

1 -2 , G -G F ≤ min r∈{1...,d} 13rB(λ 1 ) 2 + 2 d i=r+1 λ 2 i , where B is defined in equation (3.2) and λ 1 ≥ • • • ≥ λ d are the eigenvalues of G.
Proof. During the whole proof, we will assume that the event of probability at least 1 -2 described in Proposition 3.1 is realized. According to Proposition 3.6, we have

G -G 2 F = d i,k=1 (λ i -λ k ) 2 p i , q k 2 ,
where

d i,k=1 (λ i -λ k ) 2 p i , q k 2 ≤   r i=1 d k=1 + r k=1 d i=1 + d i,k=r+1   (λ i -λ k ) 2 p i , q k 2 .
Since, by definition,

λ i ≤ λ i , it follows that d i,k=r+1 (λ i -λ k ) 2 p i , q k 2 ≤ d i,k=r+1 (λ 2 i + λ 2 k ) p i , q k 2 ≤ 2 d i=r+1 λ 2 i .
Furthermore, we observe that

r k=1 d i=1 (λ i -λ k ) 2 p i , q k 2 ≤ 2 r k=1 d i=1 (λ i -λ k ) 2 p i , q k 2 + 2 r k=1 d i=1 B( λ k ) 2 p i , q k 2 ,
where, by Lemma 3.4,

r k=1 d i=1 (λ i -λ k ) 2 p i , q k 2 ≤ rB(λ 1 ) 2 .
and B( λ k ) ≤ B( λ 1 ). We have then proved that

r k=1 d i=1 (λ i -λ k ) 2 p i , q k 2 ≤ 2rB(λ 1 ) 2 + 2rB( λ 1 ) 2 .
Applying Corollary 3.3 and using the fact that B(t + a) ≤ B(t) + a/2, as explained after equation (3.2) on page 80, we deduce that

B( λ 1 ) ≤ B λ 1 + B(λ 1 ) ≤ 3B(λ 1 )/2.
This proves that

r k=1 d i=1 (λ i -λ k ) 2 p i , q k 2 ≤ 2r B(λ 1 ) 2 + 9B(λ 1 ) 2 /4 = 13rB(λ 1 ) 2 /2.
Considering that the same bound holds for

r i=1 d k=1 (λ i -λ k ) 2 p i , q k 2 ,
we conclude the proof.

To obtain a bound on f (G) -f ( G) F it is sufficient to combine the above proposition with Corollary 3.7.

Corollary 3.10. With the same notation as in Proposition 3.9, with probability at least 1 -2 , for any

1/L-Lipschitz function f , f (G) -f ( G) F ≤ min r∈{1,...,d} L -1 13 rB(λ 1 ) 2 + 2 d i=r+1 λ 2 i .

Notes: How to choose r in the previous bounds

In the previous bounds, the optimal choice of the dimension parameter r depends on the distribution of the eigenvalues of the Gram matrix G. Nevertheless, it is possible to upper bound what happens when this distribution of eigenvalues is the worst possible.

Observe that

d i=r+1 λ 2 i ≤ λ r+1 Tr(G)
and also rλ r+1 ≤ Tr(G), so that

d i=r+1 λ 2 i ≤ r -1 Tr(G) 2 .
Hence, if we consider for example the case where the approximation error is evaluated in terms of the Frobenius norm, the worst case formulation of Corollary 3.10 is obtained choosing r = 2/13 Tr(G)B(λ 1 ) -1 and, in this case, it can be restate as follows.

Corollary 3.11. With probability at least

1 -2 , f (G) -f ( G) F ≤ L -1 11 Tr(G)B(λ 1 ) + 13B(λ 1 ) 2 .
This proposition shows that the worst case speed is not slower than n -1/4 . We do not know whether this rate is optimal in the worst case. We could in the same way obtain a worst case corollary for Proposition 3.8.

Chapter 4

Spectral Clustering

Based on the results proved in chapter 1, we propose a new algorithm for spectral clustering. This new algorithm can be viewed as a change of representation in a reproducing kernel Hilbert space, followed by a (greedy) classification.

Introduction

Clustering is the task of grouping a set of objects into classes, called clusters, in such a way that objects in the same group are more similar to each other than to those in other groups. Spectral clustering algorithms use the spectrum of some data-dependent matrices to perform clustering. These matrices can be either the affinity (or similarity) matrix [START_REF] Donath | Lower bounds for the partitioning of graphs[END_REF] or the Laplacian matrix [START_REF] Fiedler | A property of eigenvectors of nonnegative symmetric matrices and its application to graph theory[END_REF].

Given X 1 , . . . , X n , a set of points to cluster, the affinity matrix A ∈ M n (R) measures the similarity, in terms of distance, between each pair of points and a common definition is

A i,j = exp - X i -X j 2 2σ 2 (4.1)
where σ is a free scale parameter. The Laplacian matrix is obtained by rescaling the affinity matrix by its row sums. Among all the definitions of the Laplacian matrix, we choose the normalized Laplacian matrix

L = I -D -1/2 AD -1/2 ,
where D is the diagonal matrix whose i-th entry is D ii = j A ij . The Laplacian matrix L is positive semi-definite, its smallest eigenvalue is zero and the corresponding eigenvector is the constant vector 1 = (1, . . . , 1) .

We observe that the Laplacian matrix has the same eigenvectors as D -1/2 AD -1/2 (and related eigenvalues) and in the following we consider for simplicity

L = D -1/2 AD -1/2
and we refer to it as the Laplacian matrix.

Many spectral clustering algorithms have been proposed, e.g. [START_REF] Shi | Normalized cuts and image segmentation[END_REF], [START_REF] Meila | Learning segmentation by random walks[END_REF], [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF]. We briefly describe one of the most successful ones, introduced by Ng, Jordan and Weiss [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF].

The goal is to cluster a set of points into c classes according to their similarity. Assuming that the number c of clusters is known, the authors use simultaneously the c largest eigenvectors of the Laplacian matrix to group the data points into classes. More precisely, denoting by X 1 , . . . , X n the dataset, the algorithm goes as follows 1. Form the affinity matrix A with entries A i,j = exp -X i -X j 2 /2σ 2 if i = j and A i,i = 0.

Construct the Laplacian matrix L.

3. Compute the c largest eigenvectors of L and form the n × c matrix X whose columns are these largest eigenvectors.

4. Renormalize each row of X to have unit length and treat each row as a point in R c .

5. Cluster points according to the new representation.

Spectral clustering can also be related to Markov chains as shown in [START_REF] Meila | Learning segmentation by random walks[END_REF] and in particular the affinity matrix A can be converted into a transition probability matrix

R = D -1 A.
The Laplacian matrix L = D 1/2 RD -1/2 has the same eigenvalues as R and related eigenvectors.

Our approach relies on viewing those matrices as empirical versions of underlying integral operators.

Let X ⊂ R d be a set endowed with the probability distribution P ∈ M 1 + (X ) and let L 2 P be the space of square integrable functions on X with respect to P. The operator related to the matrix L can be interpreted as the discrete version of the integral operator L K : L 2 P → L 2 P defined by

L K f (x) = f (y) K(x, y) ( K(x, z)dP(z)) 1 2 ( K(y, z)dP(z)) 1 2 dP(y), f ∈ L 2 P ,
where K ∈ L 2 P × L 2 P is a symmetric positive semi-definite kernel on X and in particular

K(x, y) = exp - x -y 2 2σ 2
in the case where the affinity matrix A has the form described in equation (4.1).

Connections between empirical operators and their continuous counterpart have been studied in many works, e.g. [START_REF] Rosasco | On learning with integral operators[END_REF], [START_REF] Ulrike Von Luxburg | Consistency of spectral clustering[END_REF].

Starting from the kernel

K(x, y) = K(x, y) ( K(x, z)dP(z)) 1 2 ( K(y, z)dP(z)) 1 2
we introduce, in section 4.2, an ideal spectral clustering algorithm that would require the exact knowledge of the law P. It consists in replacing the projection on the c largest eigenvectors of L done in step 3 of the algorithm introduced in [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF] by computing some suitable power of the operator L K . Indeed, K is related to a Markov chain, and iterating K (or rather the operator L K ) performs some kind of soft truncation of its eigenvalues and is related to the computation of the marginal distribution at time t of the Markov chain. This leads to a natural dimensionality reduction and allows us to propose an algorithm that automatically estimates the number c of clusters when it is not known in advance.

In order to obtain a practical algorithm, all we have to do is to compute some estimate of the operator L K , using the estimators and the generalization bounds introduced in the previous chapters, as shown in section 4.3.

We conclude providing in section 4.4 some experiments in the setting of image analysis. We will show that a very simple greedy classification algorithm can be used to perform the final classification, once the change of representation has been operated.

For more details on reproducing kernel Hilbert spaces we refer to appendix C.

Description of an ideal algorithm

Let X be a compact subset of R d , or more generally of some separable Hilbert space of possibly infinite dimension, endowed with the (unknown) probability distribution P ∈ M 1 + (X ) and let K : X × X → R be a symmetric positive semi-definite kernel. We define

K(x, y) = K(x, y) 2 , x, y ∈ X ,
and we observe that K is itself a symmetric positive semi-definite kernel.

Our approach relies on interpreting clustering as a change of representation in a reproducing kernel Hilbert space, described by a change of kernel.

We start with an ideal version of this algorithm which uses the unknown distribution P. For any x, y ∈ X , 1. Form the kernel (compute the Laplace operator)

K(x, y) = K(x, y) ( K(x, z) dP(z)) 1 2 ( K(y, z) dP(z)) 1 2
2. Construct the new kernel (iterate the Markov chain)

Km (x, y) = K(y, z 1 ) K(z 1 , z 2 ) . . . K(z m-1 , x) dP ⊗(m-1) (z 1 , . . . , z m-1 ),
where m > 0 is a free parameter 3. Make a last change of kernel (normalize the kernel, or equivalently project the representation on the unit sphere)

K m (x, y) = K2m (x, y) K2m (x, x) 1 2 K2m (y, y) 1 2
4. Cluster points according to the new representation defined by the kernel K m .

We will see that this new representation sends clusters to the neighborhood of an orthonormal basis, and therefore to a neighborhood of the vertices of a simplex, making subsequent classification an easy task.

Analysis of the algorithm

We first observe that the kernel K is symmetric and positive semi-definite. According to the Moore-Aronszajn theorem ( Proposition C.14), it defines a reproducing kernel Hilbert space H and a feature map φ : X → H, such that

K(x, y) = φ(y), φ(x) H . ( 4.2)
We introduce the Gram operator G : H → H defined by

Gv = v, φ(z) H φ(z) dP(z) (4.3)
and we observe that the kernel Km can be written in terms of some power of G.

Proposition 4.1. We have

Km (x, y) = G m-1 2 φ(x), G m-1 2 φ(y) H .
Proof. Since G is a positive operator, it is sufficient to prove, by induction, that

Km (x, y) = G m-1 φ(x), φ(y) H . ( 4.4)
Case m = 1. We observe that by definition

K1 (x, y) = K(x, y) = φ(x), φ(y) H .
We now assume that equation (4.4) holds for any < m and we prove the identity for Km . We observe that, by definition,

Km (x, y) = K(y, z 1 ) K(z 1 , z 2 ) . . . K(x, z m-1 ) dP ⊗(m-2) (z 2 , . . . , z m-1 ) dP(z 1 ) = K(y, z 1 ) Km-1 (x, z 1 ) dP(z 1 ).
Applying the induction hypothesis we obtain that

Km (x, y) = φ(y), φ(z 1 ) H G m-2 φ(x), φ(z 1 ) H dP(z 1 ) = G m-2 φ(x), φ(y), φ(z 1 ) H φ(z 1 ) dP(z 1 ) H = G m-2 φ(x), Gφ(y) H , which concludes the proof, since G is self-adjoint.
As a consequence, the kernel K m rewrites as

K m (x, y) = G m-1/2 φ(x) G m-1/2 φ(x) H , G m-1/2 φ(y) G m-1/2 φ(y) H H .
Accordingly, the representation of x ∈ X defined by the kernel

K m is isometric to the representation R(x) -1 H R(x) ∈ H where R(x) = G m-1/2 φ(x) ∈ Im(G) ⊂ H. (4.5)
In the sequel of this section we describe other isometric representations of x that can be seen as different realizations of the kernel space defined by K m .

The representation of x in the kernel space defined by K m is isometric to the representation R(x) -1 R(x) ∈ L 2 P , where R(x) ∈ L 2 P is defined as

R(x)(y) = Km (x, y) = G m-1 φ(x), φ(y) H . (4.6)
This is a consequence of the fact that

R(x), R(z) L 2 P = Km (x, y) Km (y, z) dP(y) = K2m (x, z).
We now describe a third isometric representation, related to a Markov chain, that is helpful to understand the clustering effect observed in favorable situations.

Let us assume that the support of P is made of several compact connected components and let

K(x, y) = exp -β x -y 2 ,
where β > 0 is large enough. The Markov chain related to the operator

M (f ) : x → K(x, y)f (y) dP(y) K(x, z) dP(z)
has a small probability to jump from one component to another one. Thus, for reasons that we will not try to prove here, for suitable values of m, the measure M m x : f → M m (f )(x) is close to be supported by the connected component which x belongs to and more precisely it is close to the restriction of the invariant measure of the operator M restricted to this component. As a result, the function x → M m (f )(x) is approximately constant on each connected component.

More precisely, let us put

µ(x) = K(x, z) dP(z). (4.7)
Lemma 4.2. We have

M m (f )(x) = µ(x) -1/2
Km (x, y)µ(y) 1/2 f (y) dP(y).

Proof. We prove the result by induction.

Case m = 1. Recalling the definition of the kernel K we get

M f (x) = µ(x) -1/2 K(x, y)µ(y) 1/2 f (y) dP(y) = µ(x) -1/2 K1 (x, y)µ(y) 1/2 f (y) dP(y).
We now assume that the identity holds for any < m. It follows that

M m f (x) = µ(x) -1 K(x, z)µ(z) -1/2 Km-1 (z, y)µ(y) 1/2 f (y) dP(y)dP(z) = µ(x) -1/2 K(x, z) Km-1 (z, y)µ(y) 1/2 f (y) dP(y)dP(z) = µ(x) -1/2 Km (x, y)µ(y) 1/2 f (y) dP(y)
which concludes the proof.

As a consequence, the probability measure

M m x : f → M m (f )(x) has density dM m x dP (y) = µ(x) -1/2 Km (x, y)µ(y) 1/2 .
Moreover, if we denote by M (x, y) the kernel corresponding to the operator M , we observe that µ(x)M (x, y)f (y) dP(x)dP(y) = K(x, y)f (y) dP(x)dP(y) = µ(y)f (y) dP(y).

This shows that the invariant measure Q of the operator M has a density with respect to P equal to dQ dP = µ.

Lemma 4.3. We have

dM m x dQ , dM m z dQ L 2 Q = µ(x) -1/2 K2m (x, z)µ(z) -1/2 .
Proof. The proof follows from the identity

dM m x dQ (y) = µ(x) -1/2 Km (x, y)µ(y) -1/2 .
Consider the representation of x ∈ X given by R(x

) -1 R(x) ∈ L 2 Q , where R(x) = µ(x) 1/2 dM m x dQ ∈ L 2 Q .
Lemma 4.3 shows that the representation of x ∈ X in the kernel space defined by K m is also isometric to R(x).

We conclude the section observing that

K m (x, z) = K2m (x, x) -1/2 K2m (x, z) K2m (z, z) -1/2
is the cosine of the angle formed by the two vectors representing x and z. Moreover, as explained above, when β and m are chosen in a suitable way, the supports of the probability measures M m x and M m z are almost disjoint if x and z belong to two different clusters, whereas M m x and M m z are almost the same when x and z belong to the same cluster. Therefore, in this setting, the cosine K m is either close to 0 or close to 1. This means that in the Hilbert space defined by the normalized kernel K m , the c clusters are concentrated around c orthogonal unit vectors, forming the vertices of a regular simplex. Equation (4.5) shows that this simplex is necessarily (in first approximation) contained in the linear span of the c largest eigenvectors of the Gram operator.

Choice of the scale parameter

In the case when the influence kernel K is of the form

K β (x, y) = exp -β x -y 2 ,
we propose to choose β as the solution of the equation

F (β) := K β (x, y) 2 dP(x)dP(y) = h,
where h is a suitable parameter which measures the probability that two independent points drawn according to the probability P are close to each other.

Introducing the Gram operator G β : L 2 P → L 2 P defined by K β such that

G β (f )(x) = K β (x, y)f (y) dP(y), x ∈ X ,
the parameter h is equal to the square of the Hilbert Schmidt norm of G β . As reminded in appendix C ( Proposition C.2), G β has a discrete spectrum

λ 1 ≥ λ 2 ≥ • • • . Since K β (x, x) = 1, it satisfies +∞ i=1 λ i (β) = K β (x, x) dP(x) = 1, +∞ i=1 λ 2 i (β) = F (β) ≤ 1.
Therefore, F (β) governs the spread of the eigenvalues of G β . We observe that lim β→0

F (β) = 1 which implies that λ 1 (β) -→ β→0 1 whereas λ i (β) -→ β→0 0 for i ≥ 2.
On the other hand we have lim β→+∞ F (β) = 0, so that when β grows, the eigenvalues are spread more and more widely. For these reasons, the value of the parameter h = F (β) controls the effective dimension of the representation of the distribution of points P in the reproducing kernel Hilbert space defined by K β .

Experiments show that we want this effective dimension to be pretty large, meaning that we will impose a small value of the parameter h.

Estimation of the ideal algorithm by an empirical one

By empirical algorithm we mean an algorithm based on the empirical distribution

1 n n i=1 δ X i ,
where X 1 , . . . , X n is an i.i.d. sample drawn from the probability measure P. The algorithm described in section 4.2 is ideal since it depends on the unknown probability distribution P.

In this section we use the results obtained in the previous chapters to present an estimated version of this ideal algorithm and to prove its convergence by non-asymptotic bounds.

Let X be a compact subset of a separable Hilbert space. In this section we consider the Gaussian kernel

K(x, y) = exp -β x -y 2 , x, y ∈ X ,
where β > 0 is a free parameter. However the construction we are going to present holds for any kernel K = K 2 where K is a symmetric positive semi-definite kernel, normalized in such a way that K(x, x) = 1, for any x in the state space X .

In this setting the kernel K is defined by

K(x, y) = exp - β 2 x -y 2 , x, y ∈ X .
To fix β we observe that, denoting by A the matrix of entries

A i,j = exp -β X i -X j 2 , i, j = 1, . . . , n,
we have, with the notation of section 4.2.2,

F (β) = 1 n 2 Tr(A 2 ) = 1 n + f (β)
where

f (β) = 1 n 2 i,j i =j exp -2β X i -X j 2 .
We choose β as the solution of f (β) = 0.005.

Let H be the reproducing kernel Hilbert space defined by K and let φ : X → H be the corresponding feature map, so that

K(x, z) = φ(x), φ(z) 2 H .
Lemma 4.4. The density µ(x) of the invariant measure, defined in equation (4.7), can be written as

µ(x) = φ(x), v 2 H d P(v) =: N H (x), x ∈ X ,
where

P = P • φ -1 ∈ M 1 + ( H).
Proof. It is sufficient to observe that

µ(x) = K(x, z) dP(z) = φ(x), φ(z) 2 H dP(z) = φ(x), v 2 H d P(v) = N H (x).
Therefore, according to its definition on page 91, the kernel K can be expressed as

K(x, y) = N H (x) -1/2 K(x, y)N H (y) -1/2 .
Rather than an estimate of K itself, we provide an estimate of its thresholded approximation Kσ (x, y) = max{N H (x), σ} -1/2 K(x, y) max{N H (y), σ} -1/2 , σ > 0.

When the threshold σ is suitably small this approximation retains the clustering properties of K. To support this claim, at least on a heuristic basis, we remark that Kσ is related to the sub-Markov operator

M σ (f ) : x → K(x, y)f (y) dP(y) max K(x, z) dP(z), σ
and when σ is small, M σ modifies M only in the regions where the density of P is the smallest. The fact that M σ kills the diffusion of M in regions of low probability density may be expected to preserve an interesting clustering effect.

It should also be remarked that since x → K(x, y) is continuous, also x → K(x, y) dP(z) is continuous, by the Lebesgue dominated convergence theorem, and hence, as we assume that X is compact,

σ * = inf x∈supp(P) K(x, z) dP(z) > 0.
This means that M σ (f ) = M (f ) in L 2 P when 0 < σ ≤ σ * . As a consequence the Markov chains defined by the two operators (M and M σ ) have the same distribution for any starting point in the support of P and have the same eigenvalues and eigenvectors (since L 2 P identifies functions that coincide on the support of P). In other words, Kσ (x, y) = K(x, y) for any x, y ∈ supp(P).

From now on we assume that 0 < σ ≤ σ * .

As a consequence of Lemma 4.4 we can estimate µ(x) = N H (x) by the empirical estimator

µ(x) = N H (x) = 1 n n i=1 φ(x), φ(X i ) 2 = 1 n n i=1 K(x, X i ). (4.8)
We introduce the estimated kernel

K(x, y) = max{ N H (x), σ} -1/2 K(x, y) max{ N H (y), σ} -1/2
which provides an estimate of the thresholded kernel Kσ . We remark that an explicit change of feature maps allows us to identify the Hilbert space defined by the kernel K (which we have denoted by H) and the one defined by K. Indeed, we have

K(x, y) = φ(x), φ(y) H , x, y ∈ X ,
where the estimated feature map φ is defined as

φ(x) = χ(x)φ(x), (4.9) with χ(x) = N H (x) max N H (x), σ 1/2
. Further, when x ∈ supp(P) and σ ≤ σ * , we have

χ(x) = max{N H (x), σ} max N H (x), σ 1/2 . ( 4.10)
Next proposition provides a bound on the accuracy of the estimation of Kσ by K.

Let us first recall some notation. Let a > 0 and let

K = 1 + a -1 log n 72(2 + c)κ 1/2 where c = 15 8 log(2)( √ 2-1) exp 1+2 √ 2 2 and κ = sup θ∈ H E[ θ, φ(X) 2 H ]>0 E θ, φ(X) 4 H E θ, φ(X) 2 H 2
with X ∈ X a random vector of law P. We put In this setting we take R = 1 since we are working with a normalized kernel K which sends the observations on the unit sphere of the corresponding Hilbert space H.

B * (t) =      n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(
Proposition 4.5. Assume 0 < σ ≤ σ * and let

ξ(t) = B * (t) + τ * (t) 1 -τ * (t) + 1 -B * (t) + .
With probability at least 1 -2 , for any x ∈ supp(P), Kσ (x, y)

K(x, y) -1 ≤ ξ µ(x) + ξ µ(y) 2 + ξ µ(x) 2 + ξ µ(y) 2 2 χ(x) -1 ≤ ξ µ(x) + ξ µ(x) 2 2 1 - ξ µ(x) + ξ µ(x) 2 2 -1 + .
Proof. By Proposition 2.7, with probability at least 1 -2 , for any x ∈ X , we have

1 - ξ µ(x) + ξ µ(x) 2 2 ≤ 1 -ξ µ(x) 1/2 ≤ φ(x) φ(x) ≤ 1 + ξ µ(x) 1/2 ≤ 1 + 1 2 ξ µ(x) ,
which proves the second inequality. To obtain the first bound we observe that, for any x, y ∈ X , Kσ (x, y)

K(x, y) ≤ 1 + ξ µ(x) 1/2 1 + ξ µ(y) 1/2 ≤ 1 + 2 -1 ξ µ(x) 1 + 2 -1 ξ µ(y) ≤ 1 + 2 -1 ξ µ(x) + ξ µ(x) + 4 -1 ξ µ(x) ξ µ(y) ≤ 1 + 2 -1 ξ µ(x) + ξ µ(x) + 8 -1 ξ µ(x) 2 + ξ µ(y)
and moreover Kσ (x, y)

K(x, y) ≥ 1 -ξ µ(x) 1/2 1 -ξ µ(y) 1/2 ≥ 1 - ξ µ(x) + ξ µ(x) 2 2 1 - ξ µ(y) + ξ µ(y) 2 2 ≥ 1 - ξ µ(x) + ξ µ(y) 2 - ξ µ(x) 2 + ξ µ(y) 2 2 .
Replacing K by the estimated kernel K in the definition of Km on page 91, we obtain the new kernel Hm (x, y) = K(y, z 1 ) K(z 1 , z 2 ) . . . K(x, z m-1 ) dP ⊗(m-1) (z 1 , . . . , z m-1 ).

This new kernel is still not observable since it makes use again of the sample distribution P. To make the mathematical discussion simpler, we study the estimation of Hm from a second i.i.d. sample X n+1 , . . . , X 2n ∈ R d drawn according to P and independent from the first sample X 1 , . . . , X n .

In this simplified split sample setting, K can be treated as a non-random kernel since the distribution of the second sample, conditioned on the value of the first sample, P X n+1:2n |X 1:n , is the same product distribution P ⊗n as the distribution of the second sample P X n+1:2n .

Similarly to Proposition 4.1, we observe that the kernel Hm can be written as

Hm (x, y) = G m-1 φ(x), φ(y) H , (4.12)
where G : H → H is the Gram operator defined by

Gv = v, φ(z) H φ(z) dP(z). (4.13)
We remark that, since by definition χ(x) > 0, To provide an estimator of the Gram operator G we use the construction done in section 1.3. We consider a confidence region for Gθ, θ H that we denote by

ker( G) = u ∈ H, u, φ(x) 2 dP(x) = 0 = u ∈ H, u, φ(x) 2 dP(x) = 0 = ker(G).
B -(θ) ≤ G(θ), θ H ≤ B + (θ), θ ∈ H,
and we define the estimator Q of G as in equation (1.45) on page 48. Let Q = Q + be the positive part of Q so that, according to Proposition 1.33 on page 49, with probability at least 1 -2 , for any θ ∈ S H (the unit sphere of H), max{ θ, Qθ , σ} -max{ θ, Gθ , σ} ≤ 2 max{ θ, Gθ , σ}B * θ, Gθ + η, (4.15) for a small parameter η. Hence we get the estimated kernel

H m (x, y) = Q m-1 φ(x), φ(y) H ,
and we consider as an estimator of K m (x, y) = K2m (x, x) -1/2 K2m (x, y) K2m (y, y) -1/2 the kernel

H m (x, y) = H 2m (x, x) -1/2 H 2m (x, y) H 2m (y, y) -1/2 .
The accuracy of this estimation is proved in the following proposition. We introduce

E k (x, y) = H 2m (x, y) -K2m (x, y) , x, y ∈ supp(P),
and the non-normalized quadratic error on the non-normalized kernel produced by the clustering algorithm

E k (x) = E k (x, y) 2 dP(y) 1/2
, x ∈ supp(P).

We recall that K2m (x, y) = G 2m-1 φ(x), φ(y)

H
where G is the Gram operator defined in equation (4.3).

Proposition 4.6. For any x, y ∈ supp(P),

E k (x, y) ≤ max{1, χ ∞ } 2 µ(x) 1/2 µ(y) 1/2 Q 2m-1 -G 2m-1 ∞ + 2 χ -1 ∞ E k (x) ≤ max{1, χ ∞ } 2 µ(x) 1/2 Q 2m-1 -G 2m-1 ∞ + 2 χ -1 ∞ ,
where χ and µ are defined in equation (4.10) and in equation (4.7) respectively.

Proof. We observe that, by definition,

E k (x, y) = Q 2m-1 φ(x), φ(y) H -G 2m-1 φ(x), φ(y) H ≤ ( Q 2m-1 -G 2m-1 ) φ(x), φ(y) H + G 2m-1 φ(x) -φ(x) , φ(y) H + G 2m-1 φ(x), φ(y) -φ(y) H .
Recalling the definition of φ we get

E k (x, y) ≤ Q 2m-1 -G 2m-1 ∞ χ 2 ∞ φ(x) H φ(y) H + χ -1 ∞ χ ∞ + 1 φ(x) H φ(y) H .
Since, by definition, K(x, x) = 1, we obtain

φ(x) 2 H = K(x, x) = K(x, x) µ(x) 1/2 µ(x) 1/2 = 1 µ(x)
which proves the first bound. We now consider E k . We have

E k (x) = Q 2m-1 φ(x), φ(y) H -G 2m-1 φ(x), φ(y) H 2 dP(y) 1/2 ≤ Q 2m-1 -G 2m-1 φ(x), φ(y) 2 H dP(y) 1/2 + G 2m-1 φ(x) -φ(x) , φ(y) 2 H dP(y) 1/2 + G 2m-1 φ(x), φ(y) -φ(y) 2 H dP(y) 1/2 .
Recalling that

G 1/2 u 2 H = Gu, u H = u, φ(y) 2 H dP(y) ≤ χ 2 ∞ Gu, u , u ∈ H,
we conclude

E k (x) ≤ G 1/2 Q 2m-1 -G 2m-1 φ(x) H + G 1/2 G 2m-1 φ(x) -φ(x) H + χ -1 ∞ G 2m-1/2 φ(x) H .
Before providing a bound on the quantities that appear at the right-hand side of the above inequalities, we present another way to construct the estimator H m that will lead to obtain tighter bounds on the approximation error.

We recall that the kernel K m can be written as

K m (x, y) = R(x) -1 R(x), R(y) -1 R(y) L 2 P ,
where the representation R(x) is defined in equation (4.6), and we separately estimate the representation R(x) -1 R(x) and the scalar product in L 2 P . Our estimated representation of x in L 2 P is N(x) -1 R(x), where

R(x)(y) = Q m-1 φ(x), φ(y) H and N(x) = Q 2m-1 φ(x), φ(x) 1/2 H .
This representation is not fully observable because of the presence of the ideal feature map φ in the definition of R. Nevertheless, it can be used in practice since, as we will explain in the following, it is possible to derive an observable estimate of the norm R(x) -R(x) L 2 P .

More precisely we show that we can estimate the norm in Im (G) ⊂ L 2 P , where G : L 2 P → L 2 P is the Gram operator defined by the formula

G(f )(x) = K(x, y)f (y) dP(y)
and such that

G ∞ = G ∞ = 1.
Proof. We observe that, by definition, As a consequence, combining the two previous results, the kernel

R(x) = S Q m-1 φ(x) and R(x) = SG m-1 φ(x).
R(x), R(y) L 2 P = S Q m-1 φ(x), S Q m-1 φ(y) L 2 P is estimated by H 2m (x, y) = Q 2m-1 φ(x), φ(y) H
and therefore we obtain, as before, as an estimator of K m the kernel

H m (x, y) = H 2m (x, x) -1/2 H 2m (x, y) H 2m (y, y) -1/2 .
Let us mention that the representation R(x) ∈ L 2 P is isometric to the fully observed representation

R(x) = Q m-1 φ(x), x ∈ X , in the Hilbert space Im(G), • G = Im( G),
• G with the non-observed Hilbert norm u G = Gu, u 1/2 (that could be estimated by Qu, u 1/2 according to Lemma 4.8). However we do not consider the representation

Q m-1/2 φ(x) ∈ H, • H
because its study would require a bound for the operator norm Q 3/2 -G 3/2 ∞ and the one we obtained in previous chapters is less precise than what we got for Q -G ∞ itself (that is instead what we get choosing the representation R(x)).

In the following we provide non-asymptotic bounds showing that the estimated nonnormalized representation R converges towards the ideal non-normalized representation R, defined in equation (4.5), when the sample size goes to infinity.

More specifically we give non-asymptotic bounds for the non-normalized representation error in

L 2 P norm E r (x) = R(x) -R(x) L 2 P , x ∈ supp(P),
and

E c (f ) = R(x) -R(x), f 2 L 2 P dP(x) 1/2 , f ∈ L 2 P ,
that qualifies the convergence of the coordinates of the representation. We introduce these two errors because we will be able to produce a tighter bound for E c . In favorable situations the ideal representation R lives in a neighborhood of a low-dimensional space so that the speed of convergence is well reflected by the speed of convergence along those few dimensions.

Proposition 4.10. The two errors defined above are such that, for any x ∈ supp(P), any

f ∈ L 2 P , E r (x) ≤ µ(x) -1/2 χ ∞ Q m-1 -G m-1 ∞ + χ -1 ∞ E c (f ) ≤ f L 2 P χ ∞ Q m-1 -G m-1 ∞ + χ -1 ∞ ,
where χ and µ are defined in equation (4.10) and in equation (4.7) respectively.

Proof. In order to prove the first inequality, we observe that, by definition,

E r (x) = S Q m-1 φ(x) -SG m-1 φ(x) L 2 P ≤ S Q m-1 -G m-1 φ(x) L 2 P + SG m-1 φ(x) -φ(x) L 2 P .
Moreover, since Su 2

L 2 P = S * Su, u H = Gu, u H , we get S ∞ = G 1/2 ∞ = 1.
As a consequence, recalling the definition of φ, we get

E r (x) ≤ Q m-1 -G m-1 ∞ χ ∞ φ(x) H + χ -1 ∞ φ(x) H .
We now prove the second bound. Let Φ : L 2 P → L 2 P be the orthogonal projector on Im(G). Since, according to Lemma 4.9, almost surely, R(x) -R(x) ∈ Im(G), for any x ∈ X , then R(x) -R(x), f

L 2 P = R(x) -R(x), Φ(f ) L 2 P almost surely. Moreover, since Im(G) = S Im(G) , there is u ∈ Im(G) such that Φ(f ) = Su. We can then write R(x) -R(x), f L 2 P = R(x) -R(x), Su L 2 P = S Q m-1 φ(x) -G m-1 φ(x) , Su L 2 P = Q m-1 φ(x) -G m-1 φ(x), Gu H = Q m-1 -G m-1 φ(x), Gu H + G m-1 φ(x) -φ(x) , Gu H .
Therefore, similarly as before, we get

E c (f ) ≤ G 1/2 Q m-1 -G m-1 Gu H + χ -1 ∞ G m-1 φ(x), Gu 2 H dP(x) 1/2 = G 1/2 Q m-1 -G m-1 Gu H + χ -1 ∞ G m+1/2 u H ≤ χ ∞ Q m-1 -G m-1 ∞ G 1/2 u H + χ -1 ∞ G 1/2 u H .
We conclude the proof observing that

G 1/2 u H = Φ(f ) L 2 P ≤ f L 2 P ,
where we recall that G = S * S and Φ(f ) = Su.

In order to provide observable bounds for the four errors defined previously on page 100 and on page 104 we need to introduce some technical results.

Lemma 4.11. Let m > 0. We have

Q m -G m ∞ ≤ 1 + Q -G ∞ m -1 ≤ m Q -G ∞ 1 + Q -G ∞ m-1 .
Proof. We recall that G ∞ = 1. Thus, since

Q m -G m = m-1 k=0 Q k ( Q -G)G m-k-1 , we get Q m -G m ∞ ≤ Q -G ∞ m-1 k=0 Q k ∞ . Hence, since Q ∞ ≤ 1 + Q -G ∞ , we have Q m -G m ∞ ≤ 1 + Q -G ∞ m -1 = Q -G ∞ m-1 k=0 1 + Q -G ∞ k ≤ m Q -G ∞ 1 + Q -G ∞ m-1 .
We introduce

δ 1 = ξ essinf P µ + ξ essinf P µ 2 ,
where we recall that ξ is defined in Proposition 4.5 and µ(x) = K(x, z) dP(z). We also put

δ 2 = 2 max G ∞ , σ B * G ∞ + σ + η,
where B * is defined in equation (4.11), G is the Gram operator defined in equation (4.13) and η > 0 is introduced in equation (4.15) on page 99.

Proposition 4.12. Assume 0 < σ ≤ σ * . With probability at least 1 -2 , the following inequalities hold together

χ -1 ∞ ≤ δ 1 (1 -δ 1 ) + , G -G ∞ ≤ δ 1 (2 -δ 1 ) (1 -δ 1 ) 2 + , so that χ ∞ ≤ 1 -δ 1 -1 + and G ∞ ≤ 1 -δ 1 - 2
+
where χ is defined in equation (4.10). Moreover, with probability at least

1 -2 , Q -G ∞ ≤ δ 2 .
Proof. According to Proposition 4.5, we have

χ -1 ∞ ≤ δ 1 (1 -δ 1 ) + (4.17) so that χ ∞ ≤ 1 -δ 1 -1 + . Moreover, we observe that G -G u, u H = u, φ(x) 2 -u, φ(x) 2 dP(x) = (χ(x) 2 -1) u, φ(x) 2 dP(x) ≤ χ 2 -1 ∞ u, φ(x) 2 dP(x) = χ 2 -1 ∞ Gu, u . Therefore, since G ∞ = 1, G -G ∞ ≤ χ 2 -1 ∞ G ∞ ≤ 1 -χ ∞ 2 + 1 -χ ∞ .
The last inequality follows immediately from Proposition 1.33 on page 49.

Corollary 4.13. With probability at least 1 -4 , in addition to the previous inequalities,

δ 2 ≤ 2(1 -δ 1 ) -2 + B * (1 -δ 1 ) -2 + + σ + η and Q -G ∞ ≤ δ 1 (2 -δ 1 ) (1 -δ 1 ) 2 + + δ 2 . Therefore E k (x, y) ≤ µ(x) -1/2 µ(y) -1/2 2m(2δ 1 + δ 2 ) 1 + 2δ 1 + δ 2 2m-2 (1 -δ 1 ) 4m + ,
x, y ∈ supp(P),

E k (x) ≤ µ(x) -1/2 2m(2δ 1 + δ 2 ) 1 + 2δ 1 + δ 2 2m-2 (1 -δ 1 ) 4m + , x ∈ supp(P), E r (x) ≤ µ(x) -1/2 m -1/2 2δ 1 + δ 2 (1 + 2δ 1 + δ 2 ) m-2 (1 -δ 1 ) 2m-1 + , x ∈ supp(P), E c (f ) ≤ f L 2 P m -1/2 2δ 1 + δ 2 (1 + 2δ 1 + δ 2 ) m-2 (1 -δ 1 ) 2m-1 + , f ∈ L 2 P ,
where the four errors above are defined on page 100 and on page 104.

Proof. The bound for δ 2 is a consequence of Lemma 1.22 on page 41. The second bound comes from the triangle inequality applied to the operator norm

Q -G ∞ ≤ Q -G ∞ + G -G ∞ .
Combining previous bounds, we get

E k (x, y) ≤ µ(x) -1/2 µ(y) -1/2 (1 -δ 1 ) -2 (2m -1) 2δ 1 (1 -δ 1 ) 2 + δ 2 × 1 + 2δ 1 (1 -δ 1 ) 2 + δ 2 2m-2 + 2δ 1 1 -δ 1 ,
which is lower than the simplified bound stated in the corollary. A similar computation holds for E k (x), replacing µ(y) -1/2 by the constant one. Concerning E r (x) we get

E r (x) ≤ µ(x) -1/2 (m -1) (1 -δ 1 ) + 2δ 1 (1 -δ 1 ) 2 + + δ 2 × 1 + 2δ 1 (1 -δ 1 ) 2 + + δ 2 m-2 + δ 1 (1 -δ 1 ) + ,
leading to the simplified bound presented in the corollary. The computation for E c (f ) is the same replacing µ(x) -1/2 with f L 2 P .

To conclude we observe that the last renormalization to obtain the kernel H m can also be seen as a projection on the sphere. Indeed the feature map of H m can be taken to be the projection on the sphere of the feature map of H 2m .

From this kernel we derive the new classification by thresholding the distances between the points. We have already said that taking a suitable power of the kernel K, which is here approximated by K, will threshold to zero its smallest eigenvalues. This leads to a natural dimensionality reduction which allows us to automatically estimate the number c of classes. We produce a greedy classification algorithm which consists in taking an index i ∈ {1, . . . , n} at random, constructing the corresponding class

C i = {j | H m (X i , X j ) ≥ s},
where s is a fixed threshold, and then starting again with the remaining indices. In such a way we construct a sequence of indices i 1 , . . . , i c , where the number c of classes is then estimated automatically. The k-th class is hence defined by

C k = C i k \ <k C i , k = 1, . . . , c.
In some configurations, this obviously yields an unstable classification, but in practice, when the classes are clearly separated from each other, this simple scheme is successful.

Implementation of the algorithm

Here we describe a simplified algorithm freely inspired by the algorithms for which we proved theoretical learning bounds in the previous sections.

First recall that we estimate µ(x) by µ(x) as described in equation (4.8) on page 97 to form the estimated kernel

K(x, y) = max{ µ(x), σ} -1/2 K(x, y) max{ µ(y), σ} -1/2
of K. Then we construct an estimator of the kernel Hm (x, y) = K(y, z 1 ) K(z 1 , z 2 ) . . . K(x, z m-1 ) dP ⊗(m-1) (z 1 , . . . , z m-1 ), based on a second independent sample X n+1 , . . . , X 2n ∈ R d drawn according to P.

According to equation (4.12), we have

Hm (x, y) = G m-1 φ(x), φ(y) H
where H is the reproducing kernel Hilbert space defined by K, the Gram operator G is defined in equation (4.13) and the feature map φ in equation (4.9).

Instead of estimating G by Q constructed as explained on page 99 (according to equation (1.45) on page 48) we use here a simpler iterative scheme to compute Q. More precisely, we use the polarization formula to estimate the coefficients of Q in an orthonormal basis of span{ φ(X n+1 ), . . . , φ(X 2n )} and then update this basis by replacing it iteratively by a diagonalization basis of the previous estimate.

We introduce the n × n matrix M defined as

M i,j = K X n+i , X n+j , i, j = 1, . . . , n.
We denote by λ 1 ≥ λ 2 ≥ • • • ≥ λ n ≥ 0 its eigenvalues and we decompose M as

M = U DU with U U = U U = I and D = diag(λ 1 , λ 2 , . . . , λ n ). Let r be the rank of M so that λ r > 0 and λ r+1 = λ r+2 = • • • = λ n = 0. We define V i,j = λ -1/2 i U j,i , 1 ≤ i ≤ r, 1 ≤ j ≤ n.
Lemma 4.14. An orthonormal basis of span φ(X n+1 ), . . . , φ(X 2n ) is given by

q i = n j=1 V i,j φ(X n+j ), i = 1, . . . , r,
where the sum is performed in H. Moreover, we can reconstruct φ(X n+i ), for i = 1, . . . , n, as

φ(X n+i ) = r k=1 U i,k λ 1/2 k q k . ( 4

.18)

Proof. To prove the first part of the lemma, it is sufficient to observe that the number of vectors is r and that, for any i, j ∈ {1, . . . , r},

q i , q j H = V M V i,j = λ -1/2 i U U DU U i,j λ -1/2 j = λ -1/2 i D i,j λ -1/2 j = I i,j .
To prove equation (4.18), we observe that

φ(X n+i ) = r k=1 φ(X n+i ), q k q k = r k=1 n j=1 M i,j U j,k λ -1/2 k q k = r k=1 U DU U i,k λ -1/2 k q k = r k=1 U i,k λ k λ -1/2 k q k .
We start by estimating in the basis {q i } r i=1 the Gram matrix G i,j = G(q i ), q j H , i, j = 1, . . . , r.

Introducing the Gram kernel

H i,j = G φ(X n+i ), φ(X n+j ) H , i, j = 1, . . . , n,
we see, from the definition of q i , that

G = V HV . (4.19)
We then use the polarization formula to get G i,j = 1 4 G(q i + q j ), q i + q j H -G(q i -q j ), q i -q j H , i, j = 1, . . . , r.

Given a vector p = (p 1 , . . . , p n ) of real values (that will be of the form θ, φ(X n+i ) H ) and a scale parameter λ > 0 (that we choose uniform here for simplicity), we consider the solution S(p, λ) of the equation

n i=1 ψ λ S(p, λ) -1 p 2 i -1 = 0, (4.20)
where the function ψ is defined in equation (1.2) on page 16. In practice we use a few iterations of the Newton algorithm to solve this equation. We estimate

G(q i + σq j ), q i + σq j H = q i + σq j , φ(x) 2 H dP(x) = n k=1 V i,k + σV j,k φ(X n+k ), φ(x) H 2 dP(x),
where σ ∈ {-1, +1}, by N (q i + σq j ), which we approximate by

S n k=1 V i,k + σV j,k M k, , 1 ≤ ≤ n , λ . (4.21)
We then remark that for any couple of indices (i, ), such that 1 ≤ i ≤ r and 1 ≤ ≤ n,

n k=1 V i,k M k, = D -1/2 U M i = D -1/2 U U DU i, = D 1/2 U i, , so that equation (4.21) becomes S D 1/2 U i, + σ D 1/2 U j, , 1 ≤ ≤ n , λ .
Taking notations, for any r × n matrix W , we denote by C(W ) the r × r matrix of entries

C(W ) i,j = 1 4 S W i, + W j, , 1 ≤ ≤ n , λ -S W i, -W j, , 1 ≤ ≤ n , λ .
We define a first estimate Q 0 of the r × r matrix G of entries G i,j = G(q i ), q j by

Q 0 = C Π r,n D 1/2 U ,
where Π r,n is the r×n matrix of the projection on the r first coordinates, that is (Π r,n ) i,j = δ i,j , for 1 ≤ i ≤ r and 1 ≤ j ≤ n. We then start the iterative scheme. We decompose

Q 0 = O 0 D 0 O 0 , where O 0 O 0 = O 0 O 0 = I
, V = Π r,n D 1/2 U . This yields a new estimator of G = O 0 O 0 GO 0 O 0 equal to Q 1 = O 0 C O 0 Π r,n D 1/2 U O 0 .
The inductive update step is more generally the following. We decompose

Q k = O k D k O k ,
where

O k O k = O k O k = I
and D k is a diagonal matrix, and we define the new updated estimator of G as

Q k+1 = O k C O k Π r,n D 1/2 U O k .
We stop this iterative estimation scheme when Q k -Q k-1 F falls under a suitable threshold and we take as our robust estimator of G m (q i ), q j H the last update (Q m k ) ij . The choice of the number of iterations is done automatically. We denote by λ 1 ≥ λ 2 ≥ • • • ≥ λ r the eigenvalues of Q k . Let p be the maximum number of classes we consider. As it can be seen in the simulations, the choice of p is robust, meaning that p can be harmlessly overestimated. We choose the number m of iterations by solving

λ p λ 1 m ζ
where ζ > 0 is a given small parameter. The choice of ζ is also robust and 1/100 is a reasonable value for it.

Once we have estimated G m (q i ), q j H by (Q m k ) i,j , the matrix H2m

(X n+i , X n+j ) = G 2m-1 φ(X n+i ), φ(X n+j ) H , i, j = 1, . . . , n,
is approximated by

H = U D 1/2 Π r,n Q 2m-1 k Π r,n D 1/2 U = U D 1/2 Π r,n O k D 2m-1 k O k Π r,n D 1/2 U ,
according to equation (4.18).

As a last step we renormalize H to form the new kernel

C i,j = H i,i -1/2 H i,j H j,j -1/2 , i, j = 1, . . . , n. (4.22)
This renormalization can also be interpreted as the projection on the sphere of the embedded sample, or equivalently the kernel C can be seen as the correlation matrix of the embedded dataset.

A further (optional) step is to center the representation by computing

C = C -n -1 1 n,n C -n -1 C1 n,n + n -2 1 n,n C1 n,n
where 1 n,n is the n × n matrix whose entries are all equal to one, and by normalizing once again to form

C i,j = C i,i -1/2 C i,j C j,j -1/2 , i, j = 1, . . . , n. (4.23)
From the kernel C (or as an alternative from C), we derive the new classification

C i = j | C i,j ≥ s ,
where s is a threshold that we took equal to 0.1 in all our experiments. This may not be symmetric or transitive, but this is not a problem in practice. We recall that our (greedy) classification algorithm consists in taking an index i 1 ∈ {1, . . . , n} at random, form the corresponding class C i 1 and then start again with the remaining indices.

In such a way, denoting by c the number of classes, we construct a sequence of indices i 1 , . . . , i c , and we define the k-th class as

C k = C i k \ <k C i .
It is worth underlying that in our method the number of clusters is estimated automatically, while in the spectral clustering algorithms mentioned above, [START_REF] Shi | Normalized cuts and image segmentation[END_REF], [START_REF] Meila | Learning segmentation by random walks[END_REF], [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], it is assumed to be known in advance.

Empirical results

In this section we give some examples. We first show how the algorithm described in the previous section groups the points to cluster at the vertices of a simplex, simplifying the geometry of the classes. Then, we test the algorithm in the setting of image analysis.

A first example

We briefly describe the setting in which we test the algorithm introduced in the previous section. We consider an i.i.d. set of n = 900 points to cluster, whose configuration is shown in figure 4.1 and we fix the maximum number of classes p = 7.

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q qq q q q qq q q qq q q q q q qq qq q q q q q q q q q q qq q q q q q q q q q q qq qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q qq q q q q q q q qq q q q qq q qq qq -1.0 -0. final representation q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1 0 1 2 3 In figure 4.3 we plotted the first eigenvalues of M in black (joined by a solid line), the eigenvalues of its iteration G 2m-1 φ(X n+i ), φ(X n+j ) H estimated by H i,j in blue (joined by a dashed line) and the eigenvalues of the covariance matrix of the final representation (defined by C) in red (joined by a dash dotted line). We observe that the first eigenvalues of M are close to one, while there is a remarkable gap between the eigenvalues of its iteration. In particular the size of the gap is larger once we have renormalized using the matrix C.

The number of iterations is automatically estimated and it is equal to 1441.

We obtain similar results by considering the matrix C instead of C. q q q q q q q q q q 2 4 6 8 10 0.0 0.2 0.4 0.6 0.8 1.0 Eigenvalues ratios q q q q q q q q q q q q q q q q q q q q Figure 4.3: In black (solid line) the eigenvalues of M , in blue (dashed line) those of H and in red (dash-dot line) the eigenvalues of C.

We can also test the simpler clustering algorithm obtained by estimating the Gram operator by its empirical counterpart. More precisely, instead of computing the robust estimator with the iterative scheme described in the previous section, we simply construct the matrix M as in equation (4.24) and we raise it to the power m to obtain the new matrix

M i,j = (M m ) -1/2 i,i (M m ) i,j (M m ) -1/2 j,j .
We then apply the classification algorithm to the matrix M .

.4 shows that in this special case the use of a robust estimator is not needed. Indeed, in the representation defined by M , the data points are also grouped at the vertices of the simplex generated by the first eigenvectors of M and the same classification is computed.

The number of iterations is automatically estimated and it is equal to 2090. q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q qq q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qqq q q q qq q q q q q q q q q q q qq q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q qq q q q q qq q q q q q qq q q qq q q q q q q q q q -1.0 -0. Final representation q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1 0 1 2 3

Some perspectives on invariant shape analysis

In this section we present a small example about image classification. The goal of image analysis is to recognize and classify patterns. These patterns may have been transformed by some set of transformations, such as translations, rotations or scaling and are more generally affected by the conditions in which the images to be analysed have been taken. The challenge is then to find a representation that is not hampered by the variations of the same pattern from image to image. Making a representation invariant to a set of transformations is a challenging task, even in the simpler case of translations. Indeed, it is sufficient to observe that the set of functions obtained by translating a single pattern in various directions typically spans a vector space of high dimension, meaning that the shapes (here the functions) that we would like to put in the same category do not even live in a common low-dimensional subspace.

A possible approach is to study representations that leave invariant some group of transformations. For instance, the Fourier transform has translation invariant properties, since its modulus is translation invariant. However it is unstable to small deformations at high frequencies. Wavelet transforms provide a workaround. Scattering representations proposed by Mallat [START_REF] Mallat | Group invariant scattering[END_REF] compute translation invariant representations by cascading wavelet transforms and modulus pooling operators. They bring improvements for audio classification [START_REF] Andén | Multiscale scattering for audio classification[END_REF] and for image classification [START_REF] Bruna | Invariant scattering convolution networks[END_REF].

This kind of careful mathematical study has to be repeated for any kind of transformations and in image analysis the pattern transformations we would like to take into account are numerous and not easy to formalize, since they may come from changes in the perspective or in illumination, from partial occlusions, from object deformations, etc.

Coupling spectral clustering with a change of representation in a reproducing kernel Hilbert space may lead to a generic and rather radical alternative. Instead of deriving the representation of a pattern from a mathematical study, the idea is to learn the representation itself from examples of patterns that sample in a sufficiently dense way the orbits of the set of transformations at stake. Developing a convincing toolbox for unsupervised invariant shape analysis is beyond the scope of this study and it will be carried on elsewhere. Our purpose in this section is just to hint that spectral clustering, coupled with some preliminary change of representation in a reproducing kernel Hilbert space, can be a winning tool to bring down the representation of classes to a low-dimensional space. We suggest with a small example how this approach may lead to learn transformation invariant representations from datasets containing small successive transformations of a same pattern.

We briefly describe this approach to get a hint of its potential. We consider two images (figure 4.5) and we create our patterns by translating a subwindow of a given size in each image repeatedly, using a translation vector smaller than the subwindow size. In such a way we create a sample consisting of two classes of connected images, as shown in figure 4.6. This notion of translation cannot be grasped easily by a mathematical definition, because we do not translate a function (here the image defined as a function sampled on a rectangular grid of pixels), but a window of observation. Hence in this case the translation depends on the image content and it may not be easy to model in any realistic situation. We present on an example a successful scenario using first a change of representation in a reproducing kernel Hilbert space to better separate the two classes, and then spectral clustering to shrink each class to a tight blob. We show that the so called kernel trick, introduced to better separate classes in the supervised learning framework addressed by support vector machines (SVMs), also works in an unsupervised context. In this setting, we do not separate classes using hyperplanes, since we do not know class labels which would be necessary to run a SVM, but, instead, we use spectral clustering to finishing the work. Recalling the definition of k 1 , the kernel k 2 rewrites as

k 2 (x, y) = exp -2β 2 1 -exp -β 1 x -y 2
where β 1 , β 2 > 0 are obtained as described in section 4.2.2, based on the trace estimation. We define the new kernel

K(x, y) = exp -β x -y 2 H 2 ,
where the parameter β > 0 is chosen again as in section 4.2.2, and we introduce K(x, y) = K(x, y) ( K(x, z) dP(z)) .

Proceeding as already done in the previous sections, we estimate K with the kernel K and we apply the classification algorithm to K m . As in the framework of SVMs, we use the kernel trick to embed the sample in a higher-dimensional space in which the geometry of the classes is simpler. However, as already said, we do not use hyperplanes but spectral clustering to separate the two classes.

In figure 4.7 we compare the representation of the images in the initial space and in the space H 2 .

On the left we present the projection of the sample onto the space spanned by the first two largest eigenvectors of the matrix of inner products between images X i , X j . On the right we plot the projection onto the space spanned by the two largest eigenvectors of the matrix of inner products k 2 (X i , X j ) in H 2 .

We observe that in the first representation the two classes intersect each other while in the second one, after the change of representation, they are already separated.

To conclude, figure 4.8 shows the final representation. Here the data points are projected onto the space spanned by the two largest eigenvectors of the matrix K(X i , X j ) m .

A.1 Introduction

Let P ∈ M 1 + (R d ) be an unknown probability measure and let us assume it has a density with respect to the Lebesgue measure. The goal of density estimation is to construct an estimator of the density function f from an i.i.d. sample X 1 , . . . , X n ∈ R d distributed according to P. To achieve this aim, we need to introduce the concept of approximate identities (also called approximations to the identity). Let ϕ be a function defined on R d and let > 0. We assume that ϕ is a non-negative measurable function such that ϕ(x)dx = 1.

Let us denote by L p (R d ) the space of p-integrable functions with respect to the Lebesgue measure. We now consider as approximate identities the family of Gaussian distributions π α ∼ N (0, α -1 I) of means 0 and covariance matrix α -1 I. By the properties of approximate identities, the convolution

f * π α (x) = α 2π d 2 exp - α 2 z -x 2 dP(z)
converges to f as α grows to infinity. However, since P is unknown, also f * π α is unknown.

A possible approach to estimate f * π α is to replace the distribution P by the empirical measure 1 n n i=1 δ X i . This leads to the kernel density estimator fα (x) = 1 n α 2π

d 2 n i=1 exp - α 2 X i -x 2 ,
where 1/ √ α is usually called bandwidth. We discuss the use of more general kernels in appendix A.3.

In the next section we construct, with the help of the results presented in chapter 1, a robust estimator f α of the convolution f * π α .

A.2 Estimate of the density function

In this section, our goal is construct a robust estimator of In analogy to the notation of chapter 1, we introduce

N α φ α (x) = φ α (z), φ α (x) 2 H dP(z), so that f * π α (x) = α 2π d 2
N α φ α (x) .

Let N α be the estimator of N α introduced in section 1.3 on page 46. We define

f α (x) = α 2π d 2 N α φ α (x) .
Next proposition provides a uniform bound on the approximation error f * π α (x) -f α (x) .

Proposition A.2. Let X ∈ R d be a random vector of law P. Let X be a subset of R d . Let a > 0 and let Proof. Remark that s 2 4 = E φ α (X) 4 1/2 = 1, since φ α (x) H = 1, for any x ∈ R d . To prove the result, it is sufficient to bound, uniformly on x ∈ X , the quantity max{N α φ α (x) , σ} max{ N α φ α (x) , σ} -1 .

K = 1 + a -
We can then remark that if in Proposition 1.17 on page 35 and consequently in Proposition 1.31 on page 47 we restrict the statement to any θ ∈ Θ, where Θ ⊂ R d is some subset of R d , then, in the definition of κ (equation (1.10) on page 26), we can also restrict the supremum to the same subset. Using this variant of Proposition 1.31 on page 47 where we restrict to θ ∈ φ α (X ), we obtain that with probability at least 1 -2 , for any x ∈ X , 2. S is injective 3. The operator SS * : L 2 (X, µ) → L 2 (X, µ) is such that

SS * f (x) = K(x, z)f (z) dµ(z)
and it is a positive trace class operator.

4. S * S is the expectation with respect to µ of the operator

K x ⊗ K x , that is, it is such that S * Sf, g K = f, K x K g, K x K dµ(x)
and it is a positive trace class operator.

In particular, since SS * is a self-adjoint positive trace class operator, it is compact. Hence, by the Hilbert-Schmidt therorm ( Proposition C.3), there exists an orthonormal basis {ϕ i } i of L 2 (X, µ) consisting of eigenvectors of SS * with positive eigenvalues, that is

SS * ϕ i = σ 2 i ϕ i , σ i ≥ 0.
By the Riesz-Schauder theorem ( Proposition C.2) every non-zero eigenvalue has finite multiplicity and 0 is (eventually) the only limit point. We introduce I = {i | σ i > 0} the set of indices corresponding to the non-zero eigenvalues and for any i ∈ I, we define

f i = 1 σ i S * ϕ i .
Proposition C.16. (Mercer theorem) Let K be a continuous and of positive type kernel and let K be a reproducing kernel Hilbert space with reproducing kernel K. There exists a unique choice of the orthonormal basis {ϕ i } i such that the eigenfunctions corresponding to non-zero eigenvalues are continuous. With this choice,

f i (x) = σ i ϕ i (x), i ∈ I,
and {f i } i∈I is an orthonormal basis of K. Moreover, the kernel K has the representation

K(x, x ) = i σ 2 i ϕ i (x)ϕ i (x ), x, x ∈ X,
where the convergence is absolute and uniform.

1 x

 1 + min{b, y} = log a exp min{b, y} ≤ log a + a min{b, y} exp(b) -1 b , since the function x → e x -is non-decreasing for x ≥ 0. It follows that log(a) + min{b, y} ≤ log a + ya exp(b) -1 /b .

Finally, when min{h ρ ,

 ρ 0} dρ > -∞ and max{h ρ , 0} dρ = +∞, from the Lebesgue monotone convergence theorem and the Jensen inequality, +∞ = exp h ρ dρ = sup k∈N exp min h ρ , k dρ ≤ sup k∈N exp min{h ρ , k dρ = exp h ρ dρ.

Lemma 1 . 21 .

 121 Tr(G 2 ). (1.38) Let us put η = 2δ Tr(G 2 ). With probability at least 1 -2 , for any θ ∈ S d and any (λ, β) ∈ Λ,

Lemma 1 . 25 . 4 +

 1254 Let us assume that 8ζ(σ) ≤ √ n, σ ≤ s 2 4 and that κ ≥ 3/2. With probability at least 1 -2 , for any θ ∈ S d , max{θ Qθ, σ} -max{N (θ), σ} ≤ 2 max N (θ), σ B * N (θ) + 5δ Tr(G 2 ), max{θ Qθ, σ} -max{N (θ), σ} ≤ 2 max θ Qθ, σ B * min{θ Qθ, s 2 5δ Tr(G 2 ).

k→+∞N

  Π k (θ) = N (θ) and of the continuity of B * .

Proposition 1 .

 1 40. Let N (θ) be the estimator defined in Proposition 1.38, using the entropy bound defined in terms of E Tr(A) = Tr Σ by equation (1.48) in the definition of the grid of parameters Λ. For any energy level σ ∈ R + such that σ ≤ Tr(Σ), with probability at least 1 -2 , for any θ ∈ R d ,

Figure 1 . 1 :Figure 1 . 2 :

 1112 Figure 1.1: Two data samples projected onto the two first coordinates (above) and the second and third coordinates (below). Blue circles are drawn from the most frequent distribution and red triangles from the less frequent one.

Proposition 3 . 8 .

 38 (Operator norm) With probability at least 1 -2 , for any 1/L-Lipschitz function f ,

Consequently, since

 since Im(G) = ker(G) ⊥ , we get span φ supp(P ) = ker(G) ⊥ = Im(G) = Im( G) = span φ supp(P ) . (4.14)

  = S Im(G) = S span φ supp(P) , we conclude that R(x) ∈ Im(G). Moreover, since span φ(X n ), . . . , φ(X 2n ) = span φ(X n ), . . . , φ(X 2n ) ⊂ span φ supp(P) almost surely, also R(x) ∈ Im(G).

Figure 4 . 2 :

 42 Figure 4.2: On the left the simplex generated by the eigenvectors of C, on the right classification performed on C.

Figure 4 . 4 :

 44 Figure 4.4: On the left the simplex generated by the eigenvectors of M , on the right classification performed on M .

Figure 4 . 5 :

 45 Figure 4.5: The two original images.

Figure 4 . 6 :

 46 Figure 4.6: Our sample consisting of two classes of connected images, the first sequence is obtained with a horizontal translation, the second one with a diagonal translation.

1 2 (

 2 K(y, z) dP(z)) 1 2

Proposition A. 1 .

 1 (Theorem 1.18.[START_REF] Stein | Introduction to Fourier Analysis on Euclidean Spaces[END_REF])For any f ∈ L p (R d ), p ∈ [1, +∞[, we have lim →0 f * ϕ -f L p (R d ) = 0, where ϕ (x) = -d ϕ(x/ ) are called approximate identities. Moreover, for any f ∈ L p (R d ), p ∈ [1, +∞], lim →0 f * ϕ (x) = f (x) almost surely.As a consequence, assuming f ∈ L p (R d ), the convolutionf * ϕ (x) = ϕ (z -x)f (z) dz = ϕ (z -x) dP(z)converges to f almost surely as goes to zero.

f * π α (x) = α 2π d 2 exp - α 2 z2π d 2 φ

 222 -x 2 dP(z), x ∈ R d , from the i.i.d. sample X 1 , . . . , X n ∈ R d drawn according to P. We introduce the kernelk α (x, x ) = expα 4 x -x 2 , x, x ∈ R d ,and, according to the Moore-Aronszajn theorem ( Proposition C.14), it defines a reproducing kernel Hilbert space H and a feature map φ α : R d → H such thatk α (x, x ) = φ α (x ), φ α (x) H .It follows that the convolution rewrites as f * π α (x) = α α (z), φ α (x) 2 H dP(z).

Proposition C. 15 .

 15 max{N α φ α (x) , σ} max{ N α φ α (x) , σ} -1 ≤ B * N α φ α (x) ,which concludes the proof. The following properties hold 1. S and S * are Hilbert-Schmidt operators.

  With probability at least 1 -2 , for any θ ∈ S d ,

	Proposition 1.26. Let us assume that 8ζ(σ) ≤	√	n, σ ≤ s 2 4 and that κ ≥ 3/2.

  to Corollary 1.27 we obtain the following result.

	Corollary 1.28. Assume that 8ζ(σ) ≤	√	n, σ ≤ s 2

[START_REF] Burges | A tutorial on support vector machines for pattern recognition[END_REF]

and that κ ≥ 3/2. With probability at least 1 -2 , for any θ ∈ S d ,

  Using Lemma 1.24, we deduce the analogous of Proposition 1.26.

Proposition 1.32. For any threshold σ ∈ R + such that σ ≤ s

2

4 and 8ζ(σ) ≤ √ n, in the case when κ ≥ 3/2, with probability at least 1 -2 , for any θ ∈ S H , for any k,

  The proof follows directly by applying Corollary 2.6 to (λ  , β  ) ∈ Λ defined in equation (1.27) on page 34. Indeed, by Proposition 1.17 on page 35, for any t ∈ R + ,We conclude this section by mentioning assumptions under which it is possible to give a non-random bound for R, defined in equation (2.2).

	B λ
	Notes: non-random bounds for R

where R is defined in equation (2.2) and λ * in equation

(1.31)

on page 36. With probability at least 1 -2 , for any θ ∈ S d , max{ N (θ), σ} max{N (θ), σ} -1 ≤ B * N (θ) + τ * N (θ) 1 -τ * N (θ) + 1 -B * N (θ) + . Proof.  ,β  (t) ≤ B * (t)

and, by equation (1.34) on page 36, we have λ  ≤ λ * (θ) exp(a/4).

  .5) Before stating the main result, we introduce a technical lemma.

	Lemma 3.4. With probability at least 1 -2 , for any k ∈ {1, . . . , d}, the two inequalities
	hold together
	d
	i=1

  σ} 2 , t ∈ R + ,where λ * in equation (1.31) on page 36 and R = max i=1,...,n φ(X i ) H .

		max{t, σ})	6 + (κ -1) -1 ζ(max{t, σ}) ≤	√	n	(4.11)
	+∞		otherwise
	with
	ζ(t) = 2(κ -1)	(2 + 3c) Tr( G) 4(2 + c)t	+ log(K/ ) cosh(a/4) +	2(2 + c)κ Tr( G) t	cosh(a/2)
	where Gv = v, φ(z) H φ(z) dP(z) is the Gram operator on H. According to Proposi-tion 2.7 on page 74 we define
		τ * (t) =	λ * (t) 2 exp(a/2)R 4 3 max{t,

  and D 0 is a diagonal matrix, and we observe that D 0 is an estimator ofO 0 GO 0 = O 0 V HV O 0 ,according to equation (4.19), which can be re-estimated by C O 0 Π r,n D 1/2 U , since, by definition

  1 log n 72(2 + c)κ1/2

	with c =	15 8 log(2)( √	2-1) exp 1+2 √ 2	2	and
							E exp -α x -X 2
					κ = sup x∈X	E exp -	α 2	x -X 2 2 .
	Define
	ζ(t) =	2(κ -1)	(2 + 3c) 4(2 + c)κ 1/2 t	+ log(K/ ) cosh(a/4) +	2(2 + c)κ 1/2 t	cosh(a/2)
	and the bound
	B * (t) =	  	n -1/2 ζ(max{t, σ}) 1 -4 n -1/2 ζ(max{t, σ})	6 + (κ -1) -1 ζ(max{t, σ}) ≤	√	n
			 	+∞		otherwise,
							α 2	x -X 2	,
	where				B(t) =	α 2π	d/2	max{t, σ}	B * (t) 1 -B

depending on the threshold σ. For any choice of σ ∈]0, 1], with probability at least 1 -2 , for any x ∈ X ,

max{f * π α (x), σ α/(2π) d/2 } max{ f α (x), σ α/(2π) d/2 } -1 ≤ B * E expα 2 x -X 2 , (A.1)

and consequently

f * π α (x) -f α (x) ≤ B E exp -* (t) + σ .

Chapter 2

The Empirical Gram Matrix

We consider the problem of estimating the Gram matrix via the usual empirical estimator and we provide dimension-free bounds on the approximation error, using the robust estimator introduced in chapter 1 as a tool.

Introduction

Let X ∈ R d be a random vector distributed according to the unknown probability measure P ∈ M 1 + (R d ). We consider the problem of estimating the Gram matrix

from an i.i.d. sample X 1 , . . . , X n ∈ R d distributed according to P, by using the usual empirical estimator

It is known that, by the law of large numbers, the empirical Gram matrix Ḡ converges to G as n goes to infinity.

Many results concerning the Gram matrix estimate, or more generally the covariance matrix estimate, can be found in the literature, e.g. [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF], [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF], [START_REF] Rudelson | Random vectors in the isotropic position[END_REF]. We briefly describe the approach presented in [START_REF] Tropp | User-friendly tail bounds for sums of random matrices[END_REF] and [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF].

In Vershynin [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF] the problem of estimating the Gram matrix follows from the study of spectral properties of random matrices with independent rows. The author observes that, denoting by A the matrix whose i-th row is the vector X i , the empirical Gram matrix can be expressed as

where A has, by construction, independent rows. We denote by • ∞ the operator norm.

The following result holds.

Proposition 2.1. (Corollary 5.52, [START_REF] Vershynin | Introduction to the non-asymptotic analysis of random matrices[END_REF]) Assume that the probability distribution P is supported in some Euclidean ball of radius √ m. Let δ ∈]0, 1[. With probability at least

A tool for representing functions in Im(G) is the operator S : H → L 2 P defined by S(u)(x) = u, φ(x) H , x ∈ X .

We observe that with this notation G = SS * and G = S * S (4. [START_REF] Koltchinskii | Concentration inequalities and moment bounds for sample covariance operators[END_REF] where we recall that H is the reproducing kernel Hilbert space defined by K and ψ the corresponding feature map.

Lemma 4.7. We have

Proof. We observe that, since G is symmetric,

Therefore any f ∈ Im(G) is of the form f = Su, with u ∈ Im(G), so that we can estimate

by Qu, u H and the estimation error is bounded as described in the following lemma.

Lemma 4.8. For any u ∈ H, we have

More generally, for any u, v ∈ H,

Let us now come back to R(x) and R(x). Lemma 4.9. We have

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1 0 1 In these experiments, we have not used a split sample scheme and we have recycled the same sample twice to compute the normalized kernel as well as to estimate the Gram operator. This means that in the implementation of the algorithm described in section 4.3.1, we took X n+i = X i , for any i = 1, . . . , n. With the same notation, the matrix

where A is the kernel matrix A i,j = n -1 exp -β X i -X j 2 and D is the diagonal matrix with entries

with σ = 0.001. Chapter 4. Spectral Clustering q q q q q q q q q q q q q q q q q q q q q q qqq qq q q q q q q q q q q q q q q q q q q q q q q q -0. q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q -0. On the left the projection onto the space spanned by the two largest eigenvectors of X i , X j , on the right the projection onto the space spanned by the two largest eigenvectors of k 2 (X i , X j ).

q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q -1.0 -0. Appendix A

Density Estimation

In this appendix we consider the problem of estimating the density function of an unknown probability distribution and we present an estimator, based on the results obtained in chapter 1, that is more robust than the classical kernel density estimator.

Let us remark that we can see the influence of the dimension in this result by looking at equivalents when α goes to +∞. In this case (under suitable regularity assumptions)

A.3 Kernel density estimation

Kernel density estimation is a non-parametric method of estimating the probability density function of a continuous random variable which consists in associating to each point of a given dataset a kernel function centered on that point. The kernel density estimator is defined as the sum of the centered kernel functions, scaled by a parameter h, called bandwidth, to have unit area. To be more precise, given an i.i.d. sample X 1 , . . . , X n ∈ R d and a kernel K such that K(x, y) = k(x -y), the kernel density estimator has the form

In the previous section we have considered the case of a Gaussian kernel with bandwidth 1/ √ α. We now observe that, up to the scaling parameter, the kernel density estimator is the empirical version of

In this section we provide a more robust estimator of the quantity in equation (A.2). Let X be a subset of some separable Hilbert space endowed with the unknown probability distribution P ∈ M 1 + (X ). We assume that the kernel K is such that

where K is a symmetric positive semi-definite kernel on X . Our goal is to estimate

from an i.i.d. sample X 1 , . . . , X n ∈ X drawn according to P. By the Moore-Aronszajn theorem ( Proposition C.14) K defines a reproducing kernel Hilbert space H and a feature map φ :

Thus, the quantity in equation (A.3) rewrites as

which we denote by N (φ(x)) in analogy to the notation used in chapter 1. We can see that N is the quadratic form associated with the Gram operator G : H → H defined by

A.3. Kernel density estimation 123

Let N (θ) be the estimator of the Gram operator defined in Proposition 1.31 on page 47, with θ restricted to be in φ(X ). We can estimate

by N φ(x) with the following bound on the estimation error.

Proposition A.3. With the above notation, let X ∈ X be a random vector of law P. We define, for any

where K and c are defined in Proposition A.2 and

With probability at least 1 -2 , for any x ∈ X , max{E K(x, X) , σ}

Appendix B Orthogonal Projectors

In this appendix we introduce some results on orthogonal projectors.

Let P, Q : R d → R d be two orthogonal projectors. We denote by S d the unit sphere of R d . By definition, P -Q ∞ = sup

where, without loss of generality, we can take the supremum over the normalized eigenvectors of P -Q.

A good way to describe the geometry of P -Q is to consider the eigenvectors of P + Q.

Lemma B.1. Let x ∈ S d be an eigenvector of P + Q with eigenvalue λ.

1. In the case when λ = 0, then P x = Qx = 0, so that x ∈ ker(P ) ∩ ker(Q);

2. in the case when λ = 1, then P Qx = QP x = 0, so that

x ∈ ker(P ) ∩ Im(Q) ⊕ Im(P ) ∩ ker(Q);

3. in the case when λ = 2, then x = P x = Qx, so that x ∈ Im(P ) ∩ Im(Q);

x -P x = 0, and P x, x -P x is an orthogonal basis of span x, (P -Q)x .

Proof. The operator P + Q is symmetric, positive semi-definite, and P + Q ≤ 2, so that there is a basis of eigenvectors and all eigenvalues are in the intervall [0, 2].

In case 1, 0 = P x + Qx, x = P x 2 + Qx 2 , so that P x = Qx = 0.

In case 2, P Qx = P (x -P x) = 0 and similarly QP x = Q(x -Qx) = 0.

In case 3,

In case 4, remark that P Qx = P (λx -P x) = (λ -1)P x and similarly QP x = Q(λx -Qx) = (λ -1)Qx. Consequently

Therefore (P -Q)x is an eigenvector of P + Q with eigenvalue 2 -λ = λ, so that x, (P -

Therefore x -P x 2 = x 2 -P x 2 > 0, proving that x -P x = 0. Similarly, since P and Q play symmetric roles, Qx < x and x -Qx = 0.

As P is an orthogonal projector, (P x, x -P x) is an orthogonal pair of non-zero vectors. Moreover x = x -P x + P x ∈ span{P x, x -P x}

therefore, (P x, x -P x) is an orthogonal basis of span{x, (P -Q)x}.

Lemma B.2.

There is an orthonormal basis (x i ) d i=1 of eigenvectors of P + Q with corresponding eigenvalues {λ i , i = 1, . . . d} and indices 2m ≤ p ≤ q ≤ s, such that 1.

3. x 2m+1 , . . . , x p ∈ Im(P ) ∩ ker(Q) , and

4.

x p+1 , . . . , x q ∈ Im(Q) ∩ ker(P ) , and Proof. There exists a basis of eigenvectors of P + Q (as already explained at the beginning of proof of Lemma B.1). From the previous lemma, we learn that all eigenvectors in the kernel of P + Q are in ker(P ) ∩ ker(Q), as on the other hand obviously ker(P ) ∩ ker(Q) ⊂ ker(P + Q) we get that ker(P + Q) = ker(P ) ∩ ker(Q).

In the same way the previous lemma proves that the eigenspace corresponding to the eigenvalue 2 is equal to Im(P ) ∩ Im(Q). It also proves that the eigenspace correponding to the eigenvalue 1 is included in and consequently is equal to Im(P )∩ker(Q) ⊕ ker(P )∩ Im(Q) , so that we can form an orthonormal basis of this eigenspace by taking the union of an orthonormal basis of Im(P ) ∩ ker(Q) and an orthonormal basis of ker(P ) ∩ Im(Q).

Consider now an eigenspace corresponding to an eigenvalue λ ∈]0, 1[∪]1, 2[ and let x, y be two orthonormal eigenvectors in this eigenspace. Remark that (still from the previous lemma) (P -Q)x, (P -Q)y = (P -Q) 2 x, y = (2 -λ)λ x, y = 0.

Therefore, if x 1 , . . . , x k is an orthonormal basis of the eigenspace V λ corresponding to the eigenvalue λ, then (P -Q)x 1 , . . . , (P -Q)x k is an orthogonal system in V 2-λ . If this system was not spanning V 2-λ , we could add to it an orthogonal unit vector y k+1 ∈ V 2-λ so that x 1 , . . . , x k , (P -Q)y k+1 would be an orthogonal set of non-zero vectors in V λ , which would contradict the fact that x 1 , . . . , x k was supposed to be an orthonormal basis of V λ . Therefore,

is an orthonormal basis of V 2-λ . Doing this construction for all the eigenspaces V λ such that λ ∈]0, 1[ achieves the construction of the orthonormal basis described in the lemma. This proves the claim of the lemma concerning P . Since P and Q play symmetric roles, this proves also the claim concerning Q, mutatis mutandis.

Lemma

Lemma B.4. The projectors P and Q have the same rank if and only if

p -2m = q -p.

Lemma B.5. Assume that rank(P ) = rank(Q).

Then

Proof. As P -Q is a symmetric operator, we have sup

Remark that the basis described in Lemma B.2 is also a basis of eigenvectors of (P -Q) 2 . More precisely, according to Lemma B.1

If q -2m > 0, then P -Q ∞ = 1, and q -p > 0, according to Lemma B.4, so that (P -Q)x p+1 = 1, where x p+1 ∈ Im(Q). If q = 2m and m > 0, there is i ∈ {1, . . . , m}

Since x i and x m+i are two eigenvectors of (P -Q) 2 corresponding to this eigenvalue, all the non-zero vectors in span{x i , x m+i } (including Qx i ) are also eigenvectors of the same eigenspace.

and therefore that sup θ∈S d (P -Q)θ is reached on Im(Q). Finally, if q = 0, then P -Q is the null operator, so that sup θ∈S d (P -Q)θ is reached everywhere, including on Im(Q) ∩ S d .

Appendix C

Reproducing Kernel Hilbert Spaces C.1 Operators on Hilbert spaces

For more details on the results presented in this section we refer to [START_REF] Reed | Methods of modern mathematical physics[END_REF].

Let H be a Hilbert space. We denote by L(H) the space of bounded linear operators from H to H. We are now going to introduce some criterium for determining when a given operator is compact.

Let H be a separable Hilbert space and let {e n } n be an orthonormal basis of H. For any positive operator T ∈ L(H), we define the trace of T as

Combining Proposition C.4 and Proposition C.9 we obtain that the function K can be written as

where ψ n are the eigenfunctions of T and the sum converges in L 2 (X × X, µ ⊗ µ).

C.2 Definitions and main properties

Definition C.10. Let X be a set. A (real) reproducing kernel Hilbert space on X is a Hilbert space K such that 1. the element of K are functions defined on X;

2. (reproducing property) for all x ∈ X there exists

The reproducing property is equivalent to requiring that, for any x ∈ X, there exists a positive constant C x such that

The reproducing kernel corresponding to K is the function K : X × X → R defined by

The kernel K is of positive type, that is, it is symmetric and for any x 1 , . . . , x n ∈ X and c

We now present some example of reproducing kernels:

1) Linear kernels:

Next result provides a characterization of reproducing kernel Hilbert spaces.

Lemma C.11. Let K be a reproducing kernel Hilbert space. We have

Another way to characterize reproducing kernel Hilbert spaces is via the so called feature map.

Definition C.12. Let X be a set and H be a (real) Hilbert space, a feature map is a map φ

We now show that there is a one to one correspondence between reproducing kernel Hilbert spaces and feature maps. The first result shows that a feature map φ : X → H univocally determines a reproducing kernel Hilbert space. We define the set

The proposition below shows that K endowed with the product

is a reproducing kernel Hilbert space.

Proposition C.13. Given φ : X → H a feature map, the set K defined in equation (C.3) with the inner product (C.4) is a reproducing kernel Hilbert space with reproducing kernel K(x, x ) = φ(x ), φ(x) H .

Moreover, the map W : H → K defined by

is a unitary operator and in particular K x = f φ(x) .

The space K is univocally determined by the kernel K, in the sense that, given φ : X → H another feature map such that K(x, x ) = φ (x ), φ (x) H , then {f v | v ∈ H } and K coincides as Hilbert spaces.

Vice versa, any reproducing kernel Hilbert space K defines a feature map φ : X → K putting φ(x) = K x , x ∈ X, whose feature operator is the identity, by the reproducing property.

We conclude the section stating the Moore-Aronszajn theorem.

Proposition C.14. (Moore-Aronszajn theorem) Let K : X × X → R be a kernel of positive type. There exists a unique reproducing kernel Hilbert space with K as reproducing kernel.

We then conclude that a kernel K : X × X → R is of positive type if and only if there exist a Hilbert space H and a function φ : X → H such that, for any x, x ∈ X, K(x, x ) = φ(x ), φ(x) H .

C.3 The Mercer theorem

In this section we introduce the Mercer theorem, which characterizes continuous kernels on compact domains.

Let φ : X → H be a feature map and K the corresponding reproducing kernel Hilbert space with reproducing kernel K. We are going to make the following assumptions:

1. X is a compact space which satisfies the second axiom of countability, i.e. its topology has a countable basis;

2. the reproducing kernel K is continuous;

3. µ ∈ M 1 + (X) is a probability measure such that supp(µ) = X.

We consider the continuous operator S : K → L 2 (X, µ) defined by

and we denote by S * : L 2 (X, µ) → K its adjoint. In order to state some properties of these operators, we introduce the operator K x ⊗ K x : K → K defined as