Mohamed El

Amine Seddik

Mohamed Tamaazousti
email: mohamed.tamaazousti@cea.frromain.couillet@gipsa-lab.grenoble-inp.fr

Romain Couillet

A KERNEL RANDOM MATRIX-BASED APPROACH FOR SPARSE PCA

In this paper, we present a random matrix approach to recover sparse principal components from n p-dimensional vectors. Specifically, considering the large dimensional setting where n, p → ∞ with p/n → c ∈ (0, ∞) and under Gaussian vector observations, we study kernel random matrices of the type f ( Ĉ), where f is a three-times continuously differentiable function applied entry-wise to the sample covariance matrix Ĉ of the data. Then, assuming that the principal components are sparse, we show that taking f in such a way that f (0) = f (0) = 0 allows for powerful recovery of the principal components, thereby generalizing previous ideas involving more specific f functions such as the soft-thresholding function.

INTRODUCTION

Principal component analysis (PCA) is extensively used in data analysis and machine learning applications. It is a dimension reduction technique that aims to project a given dataset onto principal subspaces spanned by the leading eigenvectors of the sample covariance matrix [START_REF] Wold | Principal component analysis[END_REF], which represent the principal modes of variance. Basically, the statistical interpretation of PCA lies in the fact that most of the variance in the data is captured by these modes. Consequently, PCA reduces the dimension of the feature space while keeping most of the information in the data. It is well-known [START_REF] Wilbur | Asymptotic theory for principal component analysis[END_REF] that PCA performs efficiently in the traditional data setting where the number of features is small and the number of samples is large.

Consider a data matrix Y ∈ R p×n consisting of n centered samples, each sample having p features. The standard PCA method requires the computation of the sample covariance matrix Ĉ = Y Y /n and estimates the first principal components u 1 , u 2 , . . . (i.e., the successive dominant eigenvectors of C = E[Y Y /n]) by the ordered eigenvectors û1 , û2 , . . . of Ĉ. [START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF] demonstrated that, in the high dimensional regime where n, p → ∞ with p/n → c > 0, the principal component û1 estimated by standard PCA is inconsistent. Essentially, if p/n 0 then û1u 1 2 0 in the high-dimensional asymptotic regime. This phenomenon is well investigated within the field of random matrix theory for covariance models of the form

Ĉ = 1 n Σ 1/2 p XX Σ 1/2
p , where Σ p is a positive semi-definite matrix and X is a p × n matrix with random i.i.d. entries. One of the main results from random matrix theory concerns the so-called spiked models, where Σ p is a low-rank perturbation of the identity matrix, namely Σ p = I p + k i=1 ω i u i u i with k fixed with respect to p, n. [START_REF] Baik | Phase transition of the largest eigenvalue for nonnull complex sample covariance matrices[END_REF] and [START_REF] Paul | Asymptotics of sample eigenstructure for a large dimensional spiked covariance model[END_REF] notably exhibited a phase transition phenomenon: as p/n → c, if ω i < √ c the estimated principal component ûi using standard PCA is (almost surely) asymptotically orthogonal to the true principal component u i (i.e., û i u i → 0); on the other hand, if ω i > √ c, lim inf n |û i u i | > 0. This phase transition phenomenon has attracted recently much attention within the random matrix community [START_REF] Benaych | The eigenvalues and eigenvectors of finite, low rank perturbations of large random matrices[END_REF][START_REF] Capitaine | The largest eigenvalues of finite rank deformation of large wigner matrices: convergence and nonuniversality of the fluctuations[END_REF][START_REF] Féral | The largest eigenvalue of rank one deformation of large wigner matrices[END_REF][START_REF] Knowles | The isotropic semicircle law and deformation of wigner matrices[END_REF].

The inconsistency of standard PCA in high dimensions motivated the idea to look for more structural information on the principal components. In particular, considering that the principal components are sparse in an appropriate basis (e.g., in the wavelet domain), a large body of works have emerged and proposed improved PCA approaches that account for sparsity. One of the most consistent sparse PCA methods in the literature is the covariance thresholding (CT) algorithm [START_REF] Krauthgamer | Do semidefinite relaxations solve sparse pca up to the information limit?[END_REF]. Based on the intuition that the small entries of the empirical covariance matrix Ĉ induce noise in its principal components, this method consists in applying the popular soft-thresholding function (with threshold τ > 0); soft(• ; τ ) : t → sign(t) • (|t|τ ) + , entry-wise to the empirical covariance matrix Ĉ and performing PCA on the resulting matrix. [START_REF] Deshpande | Sparse pca via covariance thresholding[END_REF][START_REF] Deshpande | Sparse pca via covariance thresholding[END_REF] have theoretically demonstrated that the covariance thresholding algorithm recovers the sought-for principal components with high probability under controlled growth rates between p, n and the sparsity level. In this paper, we show that the soft-thresholding method in fact falls within a broader class of kernel-based1 PCA algorithms that are particularly suited to sparse PCA recovery. This method consists in considering the matrix f ( Ĉ) instead of Ĉ where f is a function applied entry-wise. By imposing some constraints on f , most importantly that f (0) = f (0) = 0, sparse PCA can be performed with provably high accuracy for sufficiently large n, p.

The rest of the paper is organized as follows. In Section 2, we present the related work on sparse PCA. We recall some necessary concentration of measure tools and notions about sparse matrices in Section 3. Our main theoretical results are then provided in Section 4. Section 5 discusses the practical aspects and provides experimental results. Section 6 concludes the article.

Notation: In following, the notation [n] denotes the set {1, . . . , n}, a denotes the integer part of a. Given a vector x ∈ R n , the 2 -norm of x is denoted as

x 2 = n i=1 x 2 i . Given an n × n matrix M , M ij or [M ] ij
denote the entry of the matrix M at line i and column j. [M ] •,j denotes the j-th column vector of M and [M ] i,• its i-th line vector. The Frobenius norm M F of the matrix M is defined as

M 2 F = n i,j=1 M 2 ij
, and the operator norm M op of M is defined as

M op = max x =1 M x . Finally, denotes the Hadamard product, with [M N ] ij = M ij N ij .

RELATED WORK

The problem of sparse PCA has been tackled with a large range of techniques. Mainly, three classes of approaches emerge in the literature. Most popular techniques are optimization-based algorithms [START_REF] Alexandre D'aspremont | A direct formulation for sparse pca using semidefinite programming[END_REF][START_REF] Moghaddam | Spectral bounds for sparse pca: Exact and greedy algorithms[END_REF][START_REF] Zass | Nonnegative sparse pca[END_REF][START_REF] Zou | Sparse principal component analysis[END_REF][START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF], where the idea is to see the problem of sparse PCA through an optimization perspective, and to propose methods to solve the latter by either considering a different formulation e.g. semi-definite programming (SDP) or convex relaxations -or adding penalties to the original optimization problem such as a LASSO regularization. The second class of approaches covers matrix decomposition-based techniques [START_REF] Asteris | Nonnegative sparse pca with provable guarantees[END_REF][START_REF] Papailiopoulos | Sparse pca through lowrank approximations[END_REF][START_REF] Shen | Sparse principal component analysis via regularized low rank matrix approximation[END_REF], where sparse principal components are extracted through solving a low rank matrix approximation problem based on Singular Value Decomposition. Finally, most consistent sparse PCA methods adopt thresholding-based approaches: initial heuristics used factor rotation techniques and thresholding of eigenvectors to obtain sparsity [START_REF] Cadima | Loading and correlations in the interpretation of principle compenents[END_REF]. Based on the well-known power method, [START_REF] Yuan | Truncated power method for sparse eigenvalue problems[END_REF] introduced an efficient sparse PCA approximation to obtain the exact level of required sparsity, by truncating to zero the principal components iteratively except for their largest entries. A step further, under a spiked covariance model (see Section 4.1), [START_REF] Ma | Sparse principal component analysis and iterative thresholding[END_REF] proposed a very efficient iterative thresholding approach for estimating principal subspaces in the sparse setting. Similarly, assuming a single-spike model, [START_REF] Krauthgamer | Do semidefinite relaxations solve sparse pca up to the information limit?[END_REF] proved that, when the sparsity level s ≥ Ω( √ n), a standard SDP approach cannot recover consistently the sparse spike; in particular, the authors presented empirical results suggesting that for s

= O( √ n),
recovery is possible by a simple covariance thresholding algorithm. More recently, [START_REF] Deshpande | Sparse pca via covariance thresholding[END_REF] analyzed and theoretically proved, under a spiked model, that indeed the covariance thresholding algorithm [START_REF] Krauthgamer | Do semidefinite relaxations solve sparse pca up to the information limit?[END_REF] succeeds with high probability under controlled growth rates between p, n and s.

In this work, while restricting ourselves to a setting where p and n grow at a controlled joint rate, we provide an elementary argument, based on a matrix-wise Taylor expansion controlled through a concentration of measure approach, that generalizes the CT method to a large family of kernel-based methods, by means of a kernel random matrix approach [START_REF] El | The spectrum of kernel random matrices[END_REF]a). Concretely, we study kernel random matrices of the form f (Y Y /n) where Y = Σ 1/2 p X and X is a random matrix with N (0, 1) i.i.d. entries. [START_REF] El | The spectrum of kernel random matrices[END_REF] studied kernel matrices of the form f (Y Y /n) (i.e., the socalled inner-product kernel matrices), which is equivalent to the case Σ p = I p when considering the form f (Y Y /n). In particular, we elaborate from El Karoui's study by Taylor expanding f (Y Y /n) in the vicinity of Σ p entry-wise and controlling the resulting matrices via concentration arguments.

PRELIMINARIES

Before introducing our model setting we recall some definitions and notions of the concentration of measure theory [START_REF] Ledoux | The concentration of measure phenomenon[END_REF]) that are at the heart of our main results. Furthermore, we recall a definition, introduced by [START_REF] El | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF], of sparse matrices in the large-dimensional context that will also be exploited in this paper.

CONCENTRATION OF MEASURE RESULTS

We start by a definition of the notion of concentration for a real random variable.

Definition 1 (Concentration of a Random Variable). Given a function δ : R + → R + , a random variable Z is said to be δ-concentrated (around its mean) and we write Z ∈ δ, if for all t > 0, P {|Z -EZ| ≥ t} ≤ δ(t). In particular, Z is said to be normally (resp., exponentially) concentrated when δ

(t) = Ce -c t 2 (resp., δ(t) = Ce -c t ) and we write Z ∈ CN (c •) (resp., Z ∈ CE(c •))
, where C, c > 0 are some absolute constants.

In particular, δ-concentration remains stable by application of Lipschitz functions:

Proposition 1 (Concentration of Lipschitz Functions). Given a λ-Lipschitz function f : R → R and a concentrated random variable Z ∈ δ, we have f (Z) ∈ δ (•/λ).
As a consequence, linear combinations of δ-concentrated random variables remain concentrated. However, products of δ-concentrated random variables are more technical to handle, but we still have the following proposition in the case of normally concentrated random variables and which will be essential in this article.

Proposition 2 (Square of Normally Concentrated Random Variables). Given Z ∈ CN (c •), the random variable Z 2 is exp-normally concentrated, precisely

Z 2 ∈ K C E c 2 • + K C N c 16 E[Z] 2 • , (1)
where

K C > 0 is a constant depending only on C.
The extension of the notion of concentration to random vectors Z ∈ R p demands that R p → R Lipschitz functions are concentrated random variables.

Definition 2 (Concentration of a Random Vector). Given a function δ : R + → R + and a normal space (E, . ), a random vector

Z ∈ E is said to be δ-concentrated if for any 1-Lipschitz function f : E → R, the random variable f (Z) is δ-concentrated. We note again Z ∈ δ.
In particular, we have the concentration of Gaussian random vectors in the sense of Definition 2 in the following proposition (Tao, 2012, Theorem 2.1.12).

Proposition 3 (Normal Concentration of Gaussian Random Vectors). A Gaussian vector Z ∈ R p , with independent and identically distributed N (0, 1) entries, is normally concentrated independently on the dimension p. Furthermore, Z ∈ 2N (•/2). Remark 1. According to Definition 2, given a Lipschitz application F : R p → R q for q ∈ N * , Proposition 3 provides the normal concentration of all the random vectors F (Z). In particular, note that our results are extensible to this family of vectors and random vectors with independent entries.

ε-SPARSE MATRICES

When considering a large-dimensional random matrix setting, the notion of sparsity for such matrices is particularly attached to the choice of the matrix norm.2 [START_REF] El | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF] introduced a definition (ε-sparsity) for sparsity of matrices that is compatible with spectral analysis, and specifically adapted to the operator norm. The ε-sparsity definition requires some notions from graph theory that we present in the following: to each p × p symmetric matrix M , we define its corresponding adjacency matrix as

A(M ) = {1 Mij =0 } p i,j=1
, which corresponds to a graph G p with p vertices. A walk is said to be closed on this graph if it starts and finishes at the same vertex and the number of edges traversed by a walk defines the length of this walk. Denote C p (k) the set of closed walks of length k on G p .

Definition 3 (ε-sparse matrices [START_REF] El | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF], Definition 1)). A sequence of covariance matrices {Σ p } ∞ p=1 is said to be ε-sparse if the sequence of their associated graphs

{G p } ∞ p=1 satisfies, for all k ∈ 2N, |C p (k)| ≤ C k p ε(k-1)+1 where ε ∈ [0, 1], C k > 0 independent of p and |S| denotes the cardinality of the set S.
The ε-sparsity is both useful and convenient to out study for the following reasons: 1) it is adapted to the analysis of the operator norm of large sparse matrices (as we give concentration results on the operator norm); 2) it is also more general than other sparsity notions such as in [START_REF] Peter | Covariance regularization by thresholding[END_REF]. In the latter, the authors developed a natural permutation-invariant notion of sparsity which is more specific than Definition 3 as pointed out in the introduction of their article. Furthermore, note that both sparsity notions (Definition 3 and the one in [START_REF] Peter | Covariance regularization by thresholding[END_REF]) provide equivalent bounds for ε < 1 2 and when considering the large dimensional p ∼ n setting (see subsection 2.4 in [START_REF] Peter | Covariance regularization by thresholding[END_REF]); this is precisely the setting considered in Corollary 2 introduced subsequently (cf. µ > 0).

Remark 2. As Definition 3 is based on a graph defined by its corresponding adjacency matrix, we have the following property: given an ε-sparse matrix M and a function f such that f (0) = 0 and f (x) = x =0 0, the matrix f (M ), resulting from the application of f entry-wise to M , remains ε-sparse; this is simply a consequence of A(M ) = A(f (M )).

MAIN RESULTS

In this section, we first present the setting of the article. Then, we provide an asymptotic equivalent to the matrix f ( Ĉ). Finally, we treat as a special case the application of our result to the context of sparse PCA.

GENERAL SETTING AND MAIN RESULTS

Consider a data matrix Y ∈ R p×n defined as

Y ≡ Σ 1/2 p X = (I p + P ) 1/2 X, (2)
where X ∈ R p×n is a random matrix with i.i.d. N (0, 1) entries,

P = k i=1 ω i u i u i and U = [u 1 , . . . , u k ] ∈ R p×k is isometric.
Here, k refers to the number of principal components (or eigenvectors) u 1 , . . . , u k ∈ R p to be evaluated, with ω 1 > . . . > ω k > 0 the corresponding eigenvalues respectively. We define the quantity3

β p ≡ max i [Σ 1/2 p ] •,i . Assumptions: There exists B > 0 independent of p, n such that max ij |[Σ p ] ij | < B.
Besides, there exists > 0 such that β p ≤ B n 1 4 -for all p, n and for some absolute constant B > 0. Under these assumptions, our main technical result is as follows:

Theorem 1 (Asymptotic Equivalent). For f a three-times continuously differentiable function, define the matrices F and F respectively by4

F ≡ f 1 n Y Y ij p i,j=1 , F ≡ f (Σp) + 2 k=1 f (k) (Σp) k! Σ 1/2 p 1 n XX -Ip Σ 1/2 p k .
Then for η > 0 and for an absolute constant C > 0, we have with probability at least 1η

||F -F || op ≤ C β 6 p p n 3/2 √ η . ( 3
)
For a general smooth function f , the kernel random matrix f ( Ĉ) is particularly difficult to analyze through the usual tools of random matrix theory, such as the moment or Stieltjes transform-based methods [START_REF] Tao | Topics in random matrix theory[END_REF]. Rather than directly analyzing such a kernel random matrix, Theorem 1 gives an asymptotic equivalent to it, in operator norm, that has mainly two properties. First, the approximation matrix F contains "simple" objects that have already been analyzed in random matrix theory -in particular, the term (XX /n -I p ) in the expression of F . Second, the approximation in operator norm implies (by Weyl's inequality (Eisenstat & Ipsen, 1998, Theorem 4.1)) that, when F -F op → 0, F and F have the same eigenvalues and same "isolated" eigenvectors asymptotically (see Corollary 2 subsequently).

Sketch of Proof of Theorem 1. The main idea of the proof relies on the following intuition: for large n, the entries of XX /n -I p and its successive Hadamard products tend to zero at controllable rate. The concentration of measure framework then allows for the control of non-linear functions of the entries of XX /n -I p . Of utmost importance to this end is the following lemma.

Lemma 1 (A Concentration Result). For all i, j ∈ [p], the bilinear form g ij (X) ≡ [Σ 1/2 p ] i,• 1 n XX -I p [Σ 1/2 p ] •,j satisfies g ij (X) ∈ KE c 1 n β 2 p • + KN c 2 n β 4 p • , (4)
for some absolute constants c 1 , c 2 , K > 0.

Proof. Denoting by v i the i-th column vector of the matrix Σ 1/2 p , we have by the polarization identity, for all M Hermitian,

v i M v j = 1 4 [(v i + v j ) M (v i + v j ) -(v i -v j ) M (v i -v j )].
It thus suffices to prove the result for the quadratic form g

(X) = v 1 n XX -I p v where v ∈ R p . Noticing that v XX v = v X 2 and E 1 n v XX v = v v,
we need to prove the concentration of the random variable v X 2 . In fact, since v X is a Gaussian vector, by Proposition 3, v X ∈ 2N A Taylor expansion of F around f (Σ p ) then leads to controlling the operator norm of

f (3) (ξ n ) [Σ 1/2 p (XX /n -I p )Σ 1/2 p ] 3 for ξ n a matrix with entries in the set [[Y Y /n] ij , [Σ p ] ij ] (or [[Σ p ] ij , [Y Y /n] ij ]
). This follows precisely from exploiting Lemma 1 twice, to control the fluctuations of the entries of both ξ n (by the conditions on

max ij |[Σ p ] ij | and β p ) and [Σ 1/2 p (XX /n - I p )Σ 1/2
p ] 3 , with the bound provided in the theorem statement, thereby completing the proof.

A detailed proof of Theorem 1 is provided in Section A.2 of the Appendix. From now on, to simplify our arguments, we make the following assumptions:

Assumptions: As n → ∞, A1 p/n → c ∈ (0, ∞), A2 lim sup n max i ω i < ∞; specifically lim sup n β p < ∞.
Under this setting, we have the following important corollary to Theorem 1.

Corollary 1. Define the matrices F and F as in Theorem 1 and Assumptions A1 and A2 hold. Then, for η > 0

F = F + O η (n -1 2 ), (5)
where the notation X = O m η (n -α ) stands for the fact that P X op ≥ C n -α η -1 2m ≤ η for some absolute constant C > 0 and non-negative integer m.

As a consequence of Corollary 1, we have, by the sin(Θ) theorem of [START_REF] Davis | The rotation of eigenvectors by a perturbation. iii[END_REF], the corollary below concerning the eigenvectors of the matrices F and F .

Corollary 2. Let v 1 , . . . , v k and ṽ1 , . . . , ṽk denote respectively the k principal eigenvectors of F and F . Denote by

∆ i = ω i -ω i+1 for i ∈ [k -1].
Then for η > 0, we have

max i∈[k] min s∈{+1,-1} ∆ 2 i v i -sṽ i 2 = O η (n -1 ). (6)

SPECIAL CASE: SPARSE PCA

To get an insight on our coming results, consider the scenario where U contains finitely many non-zero entries. In this case, the perturbation matrix P in Eq. equation 11 contains finitely many non-zero entries (say s) on each line and a simple enumeration shows that |C p (k)| ≤ p s k-1 , thus P is 0-sparse in the sense of Definition 3. Similarly, I p is 0-sparse and by the additive stability5 of the ε-sparsity notion, Σ p remains 0-sparse. More generally, if we assume that it exists ε ∈ [0, 1 2 ) such that the population covariance matrix Σ p is ε-sparse, we have the following set of consequences. By Corollary 1, choosing f in such a way that f (0) = f (0) = 0 ensures that the terms f (Σ p ) . . . and f (Σ p ) . . . vanish in the expression of F . Indeed, for k ∈ {1, 2} (i) Only finitely entries of f (k) (Σ p ) do not vanish, precisely by Remark 2, since A(f (k) (Σ p )) = A(Σ p ),6 the matrix f (k) (Σ p ) is also (almost) ε-sparse.

(ii) The matrix

F (k) = Σ 1/2 p 1 n XX -I p Σ 1/2 p k has entries of order O(n -k/2
). As a result,7 we have for η > 0 and for all m > 0, max i,j |F

(k) ij | = O m η n -k 2 + 1 m .
Since in addition the operator norm of Σ

1/2 p (XX /n -I p )Σ 1/2 p
is typically of order O(1) (see e.g., [START_REF] Bai | No eigenvalues outside the support of the limiting spectral distribution of large-dimensional sample covariance matrices[END_REF])), it is then easily seen that, for each k ≥ 1, the operator norm of the Hadamard product f (k) (Σ p ) F (k) vanishes (see Lemma 5 in Appendix A). In particular, note that the non-zero entries of Σ p are controlled through the maximum entry of F (k) which is vanishing asymptotically, as mentioned in item (ii) above. On the opposite f (Σ p ) does not vanish since it has entries bounded away from zero (as long of course as f = 0). We precisely have the following result.

Theorem 2. Let µ > 0 and suppose Σ p is a 1 2+µ -sparse matrix. For f a three-times continuously differentiable function and for η > 0, we have for all ∈ (0,

µ 2(3+2µ) ) F = f (Σ p ) + O 1/ η n -µ 2(2+µ) + (2-1 2+µ ) s.t. f (0) = f (0) = 0. ( 7
)
Proof. See Section A.3 in Appendix A. (See Corollary 1 for the notation O m η (.)).

Remark 3. Theorem 2 gives a general result concerning the estimation of ε-sparse covariance matrices (more precisely, element-wise functionals of sparse covariance matrices). In particular, the spiked model in Eq. equation 11 with U sparse corresponds to the particular case when µ → ∞; in this case, for η > 0 and for all ∈ (0, 1/4), We considered ω 1 = 5 and thus the phase transition for standard PCA occurs at p/n = 5, thereby suggesting another phase transition for f -PCA. Curves obtained from 500 realizations of X.

F = f (Σ p ) + O 1/ η (n -1 2 +2 ).
One may then perform a PCA on F for some function f with f (0) = f (0) = 0 (we denote by f 1,2 (0) = 0 these two conditions in the following). But, while Σ p = I p +P is a low rank perturbation of the identity (therefore having only k eigenvalues strictly greater than 1), f (Σ p ) is likely more complex and not a mere low rank deformation of the identity. Now, if Σ p has all its non-zero entries greater than a certain threshold τ , an appropriate choice for f that avoids the deformation of I p + P is such that f 1,2 (0) = 0 and f (t) = t for all |t| > τ .

Such a convenient choice is

f (t) = t(1 -e -at 2 ), (8)
for some a > 0. This function notably satisfies

f (t) = 1 + e -at 2 (2at 2 -1) ⇒ f (0) = 0, f (t) = -2ate -at 2 2at 2 -3 ⇒ f (0) = 0. -2 0 2 -2 0 2 t f (t) f (t) f (t)
The figure above depicts the function f along with its derivatives for a = 1. Note that a compromise in the choice of a must be made that both maintains a close approximation of the identity by f on a large range and rather small values of f in the vicinity of zero. Interestingly, it can be verified that the extrema of f are independent of a but are found at ± 3 2a and thus smaller values of a create sharper f in the vicinity of zero. Similarly, the extrema of f are found at

± 3± √ 6 2a ∝ 1/ √
a, and precisely given by the four values 2 3a(3

± √ 6)e -1 2 (3± √ 6) ∝ √ a.
Thus smaller a induce larger maxima for f but no sharper slope.

EXPERIMENTS

In this section, we provide some experiments in the context of sparse PCA, where we consider the spiked model presented in Section 4.1. Precise setting given in caption of Figure 1. The spectrum of the sample covariance matrix (in gray) is quite different from that of Σ p . One instead observes a "bulk" of eigenvalues spread in the vicinity of 1. Furthermore, one observes a gap between the true spike and the estimated spike (in red) through the sample covariance matrix. This phenomenon is well-understood in random matrix theory. In particular, the extreme eigenvalue in our setting converges almost surely to the quantity (1 + ω 1 ) 1 + c ω 1 , where we recall that c = lim n p/n. However, thanks to sparsity, the spectrum of F = f ( Ĉ) closely matches that of Σ p , as suggested by Theorem 2. In particular, the extreme eigenvalue, which corresponds to the principal component, is consistently estimated. Figure 1 (right) depicts the alignment between the estimated principal 0 1,000 2,000

PCA n = 128 0 1,000 2,000 n = 256 0 1,000 2,000 n = 1024 0 1,000 2,000 f -PCA (ours) 0 1,000 2,000 0 1,000 2,000
Figure 2: Principal component recovery (in orange) by standard PCA (up) and our method (down) for the "Three Peak" example of [START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF]. The signal is sparse in the "Symmlet 8" wavelet basis. We use p = 2048, ω 1 = 5 for the strength of the spike and different values of n.

component and ground truth, by standard PCA (in blue) and our method (in black), in terms of p/n.

Our method retrieves the principal component even when the spike is not visible in the spectrum of Ĉ; namely beyond the phase transition p/n ≥ ω1. In fact, the standard PCA result is too noisy compared with the one when considering f ( Ĉ), as depicted in Figure 2. Further detailed examples are provided in Section A.4 of Appendix A, that confirm the consistency of the proposed method.

In terms of complexity, as our method consists in computing the sparse eigenvectors of a p × p matrix which can be done by power method, the complexity of estimating the principal component is about O(ps) where s is the sparsity level. And regarding the performance w.r.t. state-of-theart methods, Figure 3 depicts the performances of standard PCA, different state-of-the-art sparse PCA methods and our method, in terms of total projections score (left) and total projections error (right), for different values of the amplitudes ω i 's. We refer, in this figure, to standard PCA as PCA, TpowPCA for the method in [START_REF] Yuan | Truncated power method for sparse eigenvalue problems[END_REF], ITSPCA for the method in [START_REF] Ma | Sparse principal component analysis and iterative thresholding[END_REF], CT refers to the method in [START_REF] Deshpande | Sparse pca via covariance thresholding[END_REF]) and finally we refer to our method as f -PCA. The projections S and error E are given respectively by S = 1 k k i=1 (u i ûi)

2 and E = U U -Û Û F , where U = [u 1 , . . . , u k ] are the ground truth principal components and Û = [û 1 , . . . , ûk ] are the estimated ones.

As suggested theoretically and verified experimentally, our proposed method strongly attenuates the "noise component" of the sample covariance matrix and thus consistently estimates the principal components. In particular, in term of total projections score, PCA is the most inconsistent. In general, ITSPCA, CT and our method give equivalent results. The same holds when considering the total projections error as a metric, except that TpowPCA performs inconsistently, compared to PCA, for small values of amplitudes due to the initialization step from the PCA eigenvectors. The mostly used concurrent methods to PCA in a sparse context are iterative truncated power methods (such as the TPower [START_REF] Yuan | Truncated power method for sparse eigenvalue problems[END_REF] algorithm or the ITSPCA [START_REF] Ma | Sparse principal component analysis and iterative thresholding[END_REF] approach). These algorithms, despite great observed performances, as compared to standard PCA, suffer from two limitations. First, they are usually initialized from the PCA eigenvectors themselves and may not converge to good estimates. For weak signals, PCA is so impacted by noise that the mentioned initialization limitation may lead to non convergent or dramatically erroneous outcomes of the method. The proposed approach deals precisely with this limitation by strongly attenuating the "noise component" of the sample covariance matrix. In particular, our approach gives equivalent results to the CT method while generalizing it to the class of smooth functions f such that f (0) = f (0) = 0, in the considered regime. The second limitation concerns the choice of the hyper-parameters; in fact, TPower and ITSPCA need to set up an arbitrary deterministic threshold value that maintains at each iteration step only most powerful components. The proposed method as well as CT need also to set up a "soft" parameter (a and τ respectively). But, on the basis of [START_REF] Cheng | The spectrum of random inner-product kernel matrices[END_REF][START_REF] Abla | Subspace kernel clustering of large dimensional data[END_REF], we believe that our present investigation can be extended to the asymptotically non-trivial setting where ω i = O(1/ √ p) (in which case the dominant eigenmodes scale at a similar rate with residual noise); this setting may likely allow to exhibit and estimate optimal hyper-parameter choices. Notably, this setting has already been used in [START_REF] Hafis | Random matrix asymptotics of inner product spectral clustering[END_REF] in a different context, for hyper-parameters estimation.

The proof starts by a Taylor expansion of F ij in the vicinity of [Σ p ] ij , i.e., andF (k) is the matrix with entries

F ij = 2 k=0 f (k) (σ ij ) k! F (k) ij + f (3) (ξ n ij ) 6 F (3) ij where σ ij = [Σ p ] ij , ξ n ij ∈ [Y Y /n] ij , σ ij , 9
F (k) ij ≡ [Σ 1/2 p (n -1 XX -I p )Σ 1/2 p ] k ij = g ij (X) k . We have by Lemma 1 that [Y Y /n] ij concentrates around σ ij , so that ξ n
ij is bounded by σ ij + ε, for all ε > 0, with high probability10 (note that the condition max ij |σ ij | < B ensures that σ ij is bounded and the condition on β p ensures the quasi-exponential concentration of [Y Y /n] ij around σ ij ; see considered Assumptions above), formally

P |ξ n ij | ≥ σ ij + ε ≤ P {|g ij (X)| ≥ ε} ≤ Ke -n β 2 p min(c1 ε, c 2 ε 2 β 2 p ) ≤ Ke -K n 1 2 +2 min(c1ε,K c2ε 2 n -1 2 +2 ) ≡ p n → 0,
where K > 0. And since f (3) is continuous, we deduce that f (3) (ξ n ij ) is in particular bounded by

A ≡ max x∈[σij -ε,σij +ε] |f (3) (x)|,
with probability 1p n . Knowing that the operator norm is bounded by the Frobenius norm, we look for a control of the Frobenius norm of the tailing term. We have

f (3) (ξ n ) F (3) 2 F ≤ A 2 F (3) 2 F . ( 15
) By Lemma 2, for all k ∈ N E F (k) 2 F = p i,j=1 E |g ij (X)| 2k ≤ C k p 2 β 4k p n k ,
for some absolute constant C k > 0. Thus, by Markov's inequality, we have for all η > 0

P F (k) F ≥ p β 2k p n k 2 C k η ≤ η.
Recalling Eq. equation 15, we have with probability at least 1η

f (3) (ξ n ) F (3) F ≤ C p β 6 p n 3 2 √ η . A.3 PROOF OF THEOREM 2 Assumptions: As n → ∞, A1 p/n → c ∈ (0, ∞). A2 lim sup n max i ω i < ∞; specifically lim sup n β p < ∞.
With these assumptions, we have the following corollary to Theorem 1.

Corollary 1. Define the matrices F and F as in Theorem 1 and let Assumptions A1 and A2 hold. Then, for η > 0

F = F + O η (n -1 2 ), (16)
where the notation X = O m η (n -α ) stands for the fact that P X op ≥ C n -α η -1 2m ≤ η for some absolute constant C > 0 and non-negative integer m.

A.4 FURTHER EXPERIMENTS

A.4.1 HIGHER RANK CASE

In this section, we provide further experiments by considering a rank three case and by using the "Three Peak", "Piece Poly" and "Step New" signals of [START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF], in the "Symmlet 8" wavelet basis, as principal components. We compare the estimated PCs by our method with the kernel function in equation 8 to the estimated ones through standard PCA and the CT method [START_REF] Deshpande | Sparse pca via covariance thresholding[END_REF]. As shown in Figure 4, the proposed method retrieves consistently the principal components compared to a standard PCA. In particular, we obtain results that are similar to the ones obtained by the CT method while generalizing it to the class of smooth functions with f (0) = f (0) = 0. ) recovery (in orange) with standard PCA (left), the CT method (middle) and our method (right) where the PCs are considered to be the "Three Peak", Poly" and "Step New" signals of [START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF], in the "Symmlet 8" wavelet basis. We use p = 2048, n = 1024 and the spikes strengths are set respectively as 1 = 100, ω 2 = 75 and ω 3 = 50. In particular, we note the similarity between the results obtained by our method and CT.

A.4.2 OTHER CHOICES OF THE KERNEL FUNCTION f

In this section, we consider functions of the form f (t) = αt 3 + βt 2 + γt where α, β, γ ∈ R are some parameters to fix in order to allow or not the conditions f (0) = f (0) = 0. In particular, we set different parameters choices for α, β and γ in order to validate these conditions. Figure 5 depicts different PC recovery using the f -PCA method with the considered class of functions. As we can observe from this figure, the "cleanest" signal recovery is obtained when α = 0, β = 0, γ = 0 (i.e., when f (0) = f (0) = 0) thereby validating our theoretical conditions on the kernel function f for a consistent sparse PCA recovery. Note that these conditions are necessary but not sufficient in the sense that f has to be linear for large values of t (In particular, this is the case for the function f given by equation 8). In fact, the outcome provided by f -PCA for f (t) = αt 3 with α = 0 is not optimal as the signal is deformed (due to the unverified linearity condition), compared to the GT one. . We consider the "Three Peak" example of (Johnstone & Lu, 2009) which is sparse in the "Symmlet 8" wavelet basis. We use p = 2048, n = 256 and ω 1 = 5. In particular, we notice that the "cleanest" signal is obtained when α = 0, β = 0, γ = 0 which validate our theoretical conditions f (0) = f (0) = 0.

• 2 v 2 by

 2 Remark 1 and by Definition 2 since M → v M and u → u are respectively v -Lipschitz and 1-Lipschitz functions. We get the final result by Proposition 2.

|Figure 1 :

 1 Figure 1: (Left) Spectrum of Ĉ (up) and f ( Ĉ) (bottom) for p = 2048 and n = 7500. Limiting Marčenko-Pastur density[START_REF] Vladimir | Distribution of eigenvalues for some sets of random matrices[END_REF] in blue versus spectrum of Σ p in black, with ω 1 = 2; estimated largest eigenvalue in red. (Right) Alignment between estimated PC and GT (the "Three Peak" example of (Johnstone & Lu, 2009) in the "Symmlet 8" wavelet basis), in terms of p/n. We considered ω 1 = 5 and thus the phase transition for standard PCA occurs at p/n = 5, thereby suggesting another phase transition for f -PCA. Curves obtained from 500 realizations of X.

Figure 3 :

 3 Figure3: Performances of standard PCA, different state-of-the-art sparse PCA methods and our method in term of total projections score (left) and total projections error (right) for different values of the amplitudes ω i . The PCs u i , for i ∈ [4] are the "Three Peak", "Piece Poly", "Step New" and "Sing" signals of[START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF]. We use p = 2048 and n = 1024. The soft-parameters a and τ (respectively for our method and CT) are selected by cross-validation using a validation set of size n. The selected parameters are a = 20 and τ = 0.1.

Figure 4 :

 4 Figure4: Multiple principal components (k = 3) recovery (in orange) with standard PCA (left), the CT method (middle) and our method (right) where the PCs are considered to be the "Three Peak", Poly" and "Step New" signals of[START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF], in the "Symmlet 8" wavelet basis. We use p = 2048, n = 1024 and the spikes strengths are set respectively as 1 = 100, ω 2 = 75 and ω 3 = 50. In particular, we note the similarity between the results obtained by our method and CT.

Figure 5 :

 5 Figure5: PC recovery (in orange) by f -PCA with the function f (t) = αt 3 + βt 2 + γt for different values of the parameters (α, β, γ) ∈ R 3 . We consider the "Three Peak" example of (Johnstone & Lu, 2009) which is sparse in the "Symmlet 8" wavelet basis. We use p = 2048, n = 256 and ω 1 = 5. In particular, we notice that the "cleanest" signal is obtained when α = 0, β = 0, γ = 0 which validate our theoretical conditions f (0) = f (0) = 0.

We use the kernel-based terminology to highlight that our work falls within the framework of kernel random matrices and should not be confused with the standard kernel PCA.

Considering the identity matrix (which is a sparse matrix), Ip op = 1 while Ip F = √ p → ∞.

The role of βp is to ensure the concentration of the quadratic form in equation

introduced subsequently. When Σp is a sparse matrix, βp plays the same role as the maximum spike strength in the bounds given in[START_REF] Deshpande | Sparse pca via covariance thresholding[END_REF] for the CT method.

f and f(k) are applied entry-wise and k stands for the element-wise k-th power.

See Fact .1 in[START_REF] El | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF].

Given M ∈ R p×p , its corresponding adjacency matrix is defined asA(M ) = {1 M ij =0 } p i,j=1 .

See proof of Lemma 5 in Appendix A for a proof of this result.

CONCLUSIONIn this paper, we tackled the problem of sparse PCA through a random matrix perspective thereby generalizing recent ideas to a broader kernel-based method. Our analysis of this problem has yielded insights into how the principal components can be consistently estimated. Namely, given a spiked covariance model Ĉ and a smooth function f , we gave in this paper sufficient conditions on f to consistently estimate the principal components through the matrix f ( Ĉ). Our methodology can be generalized to other sparse covariance matrix-based contexts, in the same vein as the works in[START_REF] Peter | Covariance regularization by thresholding[END_REF][START_REF] El | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF].

Defined through the corresponding adjacency matrix to Σp; given an p × p real symmetric matrix M , its corresponding adjacency matrix is defined asA(M ) = {1 M ij =0 } p i,j=1 .

The notation a, b stands for the interval [a, b] if a < b or [b, a] otherwise.

For a given asymptotic variable n, we say that an event En occurs with high probability when it exist a function ψ(n) quasi-exponentially decreasing in n such that P{En} ≥ 1 -ψ(n).

Through the identity trace(M k ) ≤ maxij |Mij| k • |Cp(k)|.

* Couillet's work is supported by the GSTATS UGA IDEX Datascience chair and the ANR RMT4GRAPH (ANR-14-CE28-0006).

A PROOFS AND FURTHER EXPERIMENTS

In this appendix we provide the proofs of the different results presented in the paper and some additional experiments that validate our findings. It includes the following items: (i) the proof of Theorem 1 (Section A.2); (ii) The proof of Theorem 2 concerning the analysis of the sparse case (Section A.3); and finally (iii) further experiments which confirm the consistency of our method and the necessity of the conditions f (0) = f (0) = 0 on the kernel function f , using the signals of [START_REF] Iain | On consistency and sparsity for principal components analysis in high dimensions[END_REF] (Section A.4).

For convenience, we make the present appendix self-contained by recalling the preliminaries and the results presented in the main paper.

A.1 PRELIMINARIES

Proposition 2 (Square of Normally Concentrated Random Variables). Given Z ∈ CN (c •), the random variable Z 2 is exp-normally concentrated, precisely

where K C > 0 is a constant depending only on C. Definition 2 (Concentration of a Random Vector). Given a function δ : R + → R + and a normal space (E, . ), a random vector Z ∈ E is said to be δ-concentrated if for any 1-Lipschitz function f : E → R, the random variable f (Z) is δ-concentrated. We note again Z ∈ δ. Proposition 3 (Normal Concentration of Gaussian Random Vectors (Tao, 2012, Theorem 2.1.12)).

A Gaussian vector Z ∈ R p , with independent and identically distributed N (0, 1) entries, is normally concentrated independently on the dimension p. Furthermore, Z ∈ 2N (•/2). Remark 1. According to Definition 2, given a Lipschitz application F : R p → R q for q ∈ N * , Theorem 3 provides the normal concentration of all the random vectors F (Z). Definition 3 (ε-sparse matrices [START_REF] El | Operator norm consistent estimation of large-dimensional sparse covariance matrices[END_REF], Definition 1)). A sequence of covariance matrices {Σ p } ∞ p=1 is said to be ε-sparse if the sequence of their associated graphs

where ε ∈ [0, 1], C k > 0 independent of p and |S| denotes the cardinality of the set S.

Remark 2. As Definition 3 is based on a graph defined by its corresponding adjacency matrix, we have the following property: given an ε-sparse matrix M and a function f such that f (0) = 0 and f (x) = x =0 0, the matrix f (M ), resulting from the application of f entry-wise to M , remains ε-sparse; this is simply a consequence of A(M ) = A(f (M )). 1

A.2 PROOF OF THEOREM 1

Setting: Consider a data matrix Y ∈ R p×n defined as

where X ∈ R p×n with i.i.d. N (0, 1) entries,

Besides, there exists > 0 such that β p ≤ B n 1 4 -for all p, n and for some absolute constant B > 0. Under these assumptions, we have Theorem 1. For f a three-times continuously differentiable function, define the matrices F and F respectively by

Then for η > 0 and for an absolute constant C > 0, we have with probability at least 1η

Proof. Before starting the proof, we need to introduce the following key lemmas:

Lemma 1 (A Concentration Result). For all i, j ∈ [p], the bilinear form

for some absolute constants c 1 , c 2 , K > 0.

Proof. Denoting by v i the i-th column vector of the matrix Σ 1/2 p , we have by the polarization identity, for all M Hermitian,

It thus suffices to prove the result for the quadratic form g

we need to prove the concentration of the random variable v X 2 . In fact, since v X is a Gaussian vector, by Proposition 3, v X ∈ 2N Lemma 2

p ] •,j , we have, for all k ∈ N and for some absolute constant C k > 0,

Proof. Given a random variable Z, we have

whenever the right hand side is finite. Applying this identity to the random variable g ij (X) with m = 2k and exploiting the concentration property in Lemma 1 yields the result.

Published as a conference paper at ICLR 2019 Theorem 2. Let µ > 0 and suppose Σ p is a 1 2+µ -sparse matrix. For f a three-times continuously differentiable function and for η > 0, we have for all ∈ (0,

Proof. The proof needs the introduction of the following two lemmas, that can be found in (El Karoui, 2008, Lemma A.1 and A.2) and which are a consequence of the ε-sparsity notion 11

Lemma 3. Given an ε-sparse p×p real symmetric matrix M and calling

Lemma 4. Given two real symmetric matrices M and N with |M ij | ≤ N ij . Then, we have

First, we show that when Σ p is ε-sparse, the Hadamard product f (k) (Σ p ) F (k) is of vanishing operator norm for k ≥ 1, precisely Lemma 5. Let µ > 0, suppose Σ p is a 1 2+µ -sparse matrix. For f a real and differentiable function, k ∈ {1, 2} such that f (k) (0) = 0 and for η > 0, we have for all ∈ (0, k(2+µ)-2 2(3+2µ) )

Proof. We start by proving that the matrix F (k) has entries of order O(n -k/2 ). In fact, we have by Lemma 2, for all m ∈ N *

E|F

where C is an absolute constant. Recalling Assumption A1 and by the union bound, we have

which implies for η > 0 and for all m > 0 max

Besides, let M be the matrix defined as

thus, one has by Lemma 4

In particular, since f (k) (Σ p ) is 1 2+µ -sparse (by Remark 2), we have by Lemma 3 and by equation 19, for some η > 0

, choosing = 1 2m < k(2+µ)-2 2(3+2µ) yields the final result.

When considering f such that f (0) = f (0) = 0, the result holds by Corollary 1 and Lemma 5. In fact, the dominant order corresponds to k = 1 in Lemma 5. Which completes the proof.


Renaud Poncelet

Anne Verroust-Blondet

Fawzi Nashashibi

Safe Geometric Speed Planning Approach for Autonomous Driving through Occluded Intersections

come

Introduction

Autonomous driving in urban areas requires behaving safely in the presence of other vehicles that can be visible or occluded by physical obstacles. If we observe human drivers, we notice that they slow down when approaching an occluded area and drive cautiously until they realize that the path is clear. In fact, they assume that at any moment a vehicle may emerge from the occluded area.

Driving in a well-defined urban road infrastructure links the vehicle trajectory strictly enough to follow the bounds which are indicated by lane markings, once the maneuver has been determined. Thus, in what follows, we assume that the vehicle's trajectory is given and we use the pathvelocity decomposition introduced by Kant and Zucker [START_REF] Kant | Toward Efficient Trajectory Planning: The Path-Velocity Decomposition[END_REF] to focus on a speed profile adaptation method that ensures the safety of the motion planning.

In what follows, we present a geometric method computing a safe speed profile for an autonomous vehicle navigating in the presence of dynamic obstacles that may be occluded. Our approach uses a visibility graph path planning method to compute a save speed profile in the Time-Path domain in which the potential risks from occlusion and limited sensing are considered. Moreover, we take into account passenger comfort by constraining the speed and the acceleration of the vehicle.

Related Research

The problem of occluded areas has been addressed in numerous recent approaches [START_REF] Hoermann | Entering crossroads with blind corners. A safe strategy for autonomous vehicles[END_REF][START_REF] Damerow | Riskbased driver assistance for approaching intersections of limited visibility[END_REF][START_REF] Bouton | Scalable decision making with sensor occlusions for autonomous driving[END_REF][START_REF] Orzechowski | Tackling Occlusions Limited Sensor Range with Set-based Safety Verification[END_REF][START_REF] Tas | Limited Visibility and Uncertainty Aware Motion Planning for Automated Driving[END_REF][START_REF] Hubmann | A POMDP Maneuver Planner For Occlusions in Urban Scenarios[END_REF][START_REF] Lin | Decision making through occluded intersections for autonomous driving[END_REF][START_REF] Narksri | Crossing Blind Intersections from a Full Stop Using Estimated Visibility of Approaching Vehicles[END_REF][START_REF] Naumann | Safe but not Overcautious Motion Planning under Occlusions and Limited Sensor Range[END_REF][START_REF] Schorner | Predictive trajectory planning in situations with hidden road users using partially observable markov decision processes[END_REF][START_REF] Yu | Occlusion-aware risk assessment for autonomous driving in urban environments[END_REF]. They can be classified in different families, according to the way they represent these regions: * 1 Inria Paris, 2 rue Simone Iff 75012 Paris FRANCE {renaud.poncelet,anne.verroust,fawzi.nashashibi}@inria.fr

Particle filters

In order to take into consideration risks from occluded areas and the limited range of sensors, some approaches use particle filters. In [START_REF] Narksri | Crossing Blind Intersections from a Full Stop Using Estimated Visibility of Approaching Vehicles[END_REF], particles are generated randomly in occluded areas with random speed. A decision to cross is made considering particle kinematics. In [START_REF] Yu | Occlusion-aware risk assessment for autonomous driving in urban environments[END_REF], the particles are also spread to represent their possible deviation from the center of the lane. Here the particles are used to compute a risk cost which combined with a speed cost allows a speed along a given path to be computed. Particle filters approaches give a probabilist representation of the environment but do not have worst case assumption.

Virtual vehicles

Other approaches model the worst case scenario by using virtual vehicles to represent all potential vehicles in occluded areas [START_REF] Damerow | Riskbased driver assistance for approaching intersections of limited visibility[END_REF][START_REF] Orzechowski | Tackling Occlusions Limited Sensor Range with Set-based Safety Verification[END_REF][START_REF] Tas | Limited Visibility and Uncertainty Aware Motion Planning for Automated Driving[END_REF][START_REF] Hubmann | A POMDP Maneuver Planner For Occlusions in Urban Scenarios[END_REF][START_REF] Schorner | Predictive trajectory planning in situations with hidden road users using partially observable markov decision processes[END_REF][START_REF] Lin | Decision making through occluded intersections for autonomous driving[END_REF]. For a risk evaluation purpose, Damerow et al. [START_REF] Damerow | Riskbased driver assistance for approaching intersections of limited visibility[END_REF] position one virtual car for each relevant lane at the occluded position closest to the upcoming intersection. Similarly, in Orzechowski et al. [START_REF] Orzechowski | Tackling Occlusions Limited Sensor Range with Set-based Safety Verification[END_REF], virtual vehicles are placed at visibility limits. A failsave maneuver is then computed, following the reachable set approach introduced in [START_REF] Althoff | Set-based prediction of traffic participants on arbitrary road networks[END_REF]. In [START_REF] Tas | Limited Visibility and Uncertainty Aware Motion Planning for Automated Driving[END_REF], virtual vehicles are used to find rules for a safe motion or a safe stop for free drive, follow drive, right-of-way maneuvers and give-way maneuvers. They consider uncertainties and the worst case scenario to solve an optimization-based planner. In [START_REF] Hubmann | A POMDP Maneuver Planner For Occlusions in Urban Scenarios[END_REF], virtual vehicles are placed at visibility limits on roads that are likely to cross the ego-vehicle's path. They are modeled with an infinite length and are considered as driving at maximum speed. In [START_REF] Schorner | Predictive trajectory planning in situations with hidden road users using partially observable markov decision processes[END_REF], probability of occurrence is computed using traffic density to position virtual vehicles in occluded areas. Virtual vehicles are supposed to arrive at the same time as the ego-vehicle at the junction as long as they keep a realistic speed.

Other models

Hoermann et al. [START_REF] Hoermann | Entering crossroads with blind corners. A safe strategy for autonomous vehicles[END_REF] introduce an unobservable regions layer in their grid-based representation of the environment to estimate whether a region is safe to drive in at maximum speed. To be safe but not to overcautious Naumann et al. [START_REF] Naumann | Safe but not Overcautious Motion Planning under Occlusions and Limited Sensor Range[END_REF] focus on conflict zones inside junctions and compute the speed by minimizing a cost. They consider the trajectories of potentially occluded objects along with the probability of their occurrence to estimate the cost caused by them.

Numerous approaches solve the problem as a Partially Observable Markov Decision Process (POMDP). In [START_REF] Bouton | Scalable decision making with sensor occlusions for autonomous driving[END_REF], possible actions are relative to acceleration: for each user, optimal belief action utility is computed, a common belief is computed for all invisible users and then the minimum belief action utility over users is computed with the sum of all actors belief action utilities. In [START_REF] Hubmann | A POMDP Maneuver Planner For Occlusions in Urban Scenarios[END_REF], an observation model is defined depending on the field of view and virtual vehicles have a specific transition model. Lin et al. [START_REF] Lin | Decision making through occluded intersections for autonomous driving[END_REF] use a belief tree search built from possible ego-vehicle actions: speed-up, keep speed or brake. Virtual vehicles are set at visibility limits of the ego-vehicle and they are considered as real vehicles. Positions, speeds, and accelerations of other vehicles are supposed to be known and all possible future paths are considered. In our work, we introduce virtual vehicles to represent occluded parts of the road. We also take into account comfort and safety [START_REF] Shalev-Shwartz | On a formal model of safe and scalable self-driving cars[END_REF] as acceleration is bounded and a safe state is planned before intersections.

Approach

In what follows we assume that a map of the environment is given with all the information on the road topology. A path s is computed off-line given the initial position and a global destination of the ego vehicle. It is assumed to be collision-free regarding static obstacles. We assume that the positions and the velocities of the other vehicles are known for a short period of time once they are visible from the ego vehicle. However, when they are close to an intersection, we assume that they can either turn right, turn left or go straight when possible. Our goal here is to adapt the speed of the ego vehicle to obtain a safe and comfortable trajectory when it follows its path in the presence of dynamic obstacles and considering occluded areas. This is done every time step in the Time-Path domain. Before presenting the whole model, we first explain our approach in the case where there are no occluded areas. Two processes are executed in parallel, with different time steps and different frame rates: 1) for each time step t i , our algorithm computes the longitudinal motion s(t) for the ego vehicle, considering the interacting vehicles, i.e. the list N of surrounding vehicles that intersect the part of the path of the ego vehicle limited to [s(t i ), s h ] during the period ∆T = [t i , t i + t h ], where t h is a time horizon and s h the visibility horizon.

2) for each time step τ j , the longitudinal motion is simulated through a car model (Section 3.3). The time interval between two consecutive time steps τ j and τ j+1 is significantly smaller than that of two consecutive time steps t i and t i+1 , as illustrated in Figure 2. In the rest of the paper, s will denote the path, v the velocity and a the longitudinal acceleration of the ego vehicle.

Computing a speed profile in the Time-Path domain

We follow [START_REF] Kant | Toward Efficient Trajectory Planning: The Path-Velocity Decomposition[END_REF] and compute a longitudinal motion s(t) for the ego vehicle during the period ∆T = [t i , t i + t h ] along the path s and work in the Time-Path domain T P = (t, s). The longitudinal motion s(t) has to be a monotonically increasing function of t with constraints on the slope and on the curvature of s(t) deduced from the constraints on the velocity v(t) and the acceleration a(t) of the ego vehicle. Moreover, in order to be a safe motion, s(t) can be deduced from the visibility graph built in T P as in [START_REF] Kant | Toward Efficient Trajectory Planning: The Path-Velocity Decomposition[END_REF], considering the set N of vehicles that may intersect s during ∆T . The computation of s(t), v(t) and a(t), t ∈ ∆T is carried out in three steps:

1. A polygonal zone, called T P obstacle, is inserted in T P for each vehicle belonging to N (Section 3.1.1)

2. A polygonal line s(t), candidate for s(t), is built in T P (this is detailed in Section 3.1.2).

3. s(t) is smoothed in order to respect the constraint on the longitudinal acceleration (described in Section 3.1.3).

Steps 2 and 3 are illustrated in Figure 3.

T P obstacles

We assume that the vehicles belonging to N have a constant velocity during the time period ∆T . Then, for each of them, a T P obstacle is built: it corresponds to the polygonal zone occupied in T P by the vehicle during the period ∆T . To ensure safety, comfort and to take into account possible changes in velocity during the time period ∆T , the occupied zone in T P by the vehicle is expanded in time of a time gap t gap and in path of a margin. In order to take into account maneuver uncertainties, we consider that the ego vehicle cannot know which maneuver another vehicle will choose when it enters a road junction before being fully engaged in one lane. Thus, for the scenario in Figure 1 in which the ego-vehicle will cross an intersection where another vehicle is coming from the right, paths a, b and c are considered and the region that may be occupied by the vehicle during ∆T is represented by the union of the T P obstacles corresponding to maneuvers a, b and c.

Visibility graph method

Once the T P obstacles have been built, a candidate path s(t) for s(t) is computed on T P . As mentioned above, it is a path composed of line segments having a positive slope bounded by the maximum speed value v M ax , joining point (0, 0) and a point belonging to the vertical line corresponding to time t i + t h or belonging to the horizontal line corresponding to s h . This path connects a subset of the vertices of the T P obstacles and it is built by a visibility graph-based method ([1]).

Achieving acceleration constraints

The candidate path s(t) for s(t) built previously respects the speed limits but needs to be smoothed as it may not satisfy the acceleration bounds. This is done by smoothing the velocity ṽ(t) with a saturation on the acceleration ã(t) of the candidate path s(t), with ṽ

(t i ) = v i : ã(t i ) = ( s(t i + dt) -s(t i ) dt -v i )/dt (1)
and, when t i < t < t i + t h : A longitudinal motion is planned every time step t i considering the period ∆T and the control of the ego vehicle is computed every τ j

ṽ(t) = s(t + dt) -s(t) dt , (2)
ã(t) = ṽ(t + dt) -ṽ(t) dt (3)
Then, the smoothed path s(t) is defined by :

s(t i ) = s 0 , v(t i ) = v i and, when t i ≤ t < t i + t h : a(t) = min(max(ã(t), a M in ), a M ax ), (4)
v(t) = min(max(ṽ(t) + a(t) × dt, 0), v M ax ), (5)
s(t + dt) = s(t) + dt × v(t) (6)
Once s(t) is computed, we have to check whether it intersects a T P obstacle, as in the example of Figure 3. In this case, a new candidate s(t) is built and new smoothed path s(t) is computed: s(t) corresponds to a polygonal line of the visibility graph that passes on the right of the intersected polygon.

Safe state before intersections

Before validating s(t) as a safe plan, we ensure that the plan does not end in an intersection to remove uncertainties about safety when crossing the intersection. We decide to let the ego vehicle drive as fast as possible to the junction while preserving a safe distance to let it brake before the intersection. In order to model this behavior in the T P domain, we reduce the visibility limit for the ego vehicle considering that the visibility horizon s h is moved just before the intersection. Moreover, to enforce the vehicle to stop when reaching the visibility limit, we add a stopgap at in the graph at :

s brake = s h -d brake (7)
where d brake is the distance required to brake when the vehicle's speed is equal to v M ax : if s(t) > s brake then s(t) = s brake as in the Fig. 3.

Car simulation

In parallel, at each time step τ j a longitudinal controller adapts the jerk of the vehicle. It uses as input ŝ(τ j ), v(τ j ), the measured values in the simulation of s(t) and v(t) at t = τ j and s(τ j ), v(τ j ), Then it computes â(τ j+1 ), v(τ j+1 ) and ŝ(τ j+1 ) for the car simulator. The initial speed considered at time t i+1 to compute the longitudinal motion for the next period in Section 3.1.3 is v i+1 = v(t i+1 ).

Tackling occlusions and limited sensor range

Let us now extend our approach to take into account the existence of occluded areas on the road. To follow what a human driver may do in this situation we introduce virtual vehicles and a virtual horizon line (Section 3.4.1) and adapt the longitudinal motion s(t) to be able to react safely without being blocked (Section 3.4.2).

Virtual vehicles and virtual horizon

To represent all the potential risks from occluded areas belonging to the lanes crossing the path of the ego vehicle, we choose to introduce virtual vehicles on the visibility borders, as in [START_REF] Damerow | Riskbased driver assistance for approaching intersections of limited visibility[END_REF][START_REF] Hoermann | Entering crossroads with blind corners. A safe strategy for autonomous vehicles[END_REF][START_REF] Orzechowski | Tackling Occlusions Limited Sensor Range with Set-based Safety Verification[END_REF][START_REF] Naumann | Safe but not Overcautious Motion Planning under Occlusions and Limited Sensor Range[END_REF][START_REF] Tas | Limited Visibility and Uncertainty Aware Motion Planning for Automated Driving[END_REF]. In order to plan for the worst case, we suppose that the virtual vehicle is of infinite length and that it drives at the maximum speed possible, except when there exists a vehicle moving in the same lane in front of the visibility boundary. In this case, we assume that the virtual vehicle has the same speed as this vehicle. Moreover, when the visibility boundary is close to a road intersection, we assume that it may follow any path, either turning left, turning right or going straight when it is possible. Thus we create the following T P obstacles: At time step t i , if a border of an occluded area belongs to a lane crossing the path of the ego vehicle during ∆T , and if this border is very close to the road junction, a T P obstacle is inserted in T P to represent the region occupied by a virtual vehicle positioned at the visibility boundary and driving at the maximum speed during ∆T or at the same speed as the previous vehicle when such a vehicle exists. As a virtual vehicle is supposed be of infinite length, the corresponding polygon is bounded by the right boundary of T S, i.e. the vertical line t = t i + ∆T (Figure 4). Areas outside the sensor range may also contain hidden obstacles. Thus, if we want to ensure a safe behavior, the ego vehicle must stop when reaching the visibility horizon value s h in its path and the horizontal line s = s h will be considered as an obstacle when computing the candidate path s(t). 3: Computation of a trajectory curve. On the extreme left, an overview of the scene. On the middle left, the candidate path s (red line) from the visibility graph method as in [START_REF] Kant | Toward Efficient Trajectory Planning: The Path-Velocity Decomposition[END_REF] and the smoothed version (black dashed line) of s intersects one T P obstacle (in red). At t = 0 the ego-vehicle is stationnary and its speed is equal to 0. On the middle right, a new candidate path s is built, which passes after the polygon in collision. On the extreme right: the visual horizon (s h ) is virtually reduced before the intersection area (in grey). A new candidate path s is built considering a safe state before the intersection thanks to the constrain s < s brake . The last smoothed version s is safe. We notice that a solution can always be found under two conditions : The initial situation have to be safe i.e the ego-vehicle can brake without entering an intersection, and obstacles can't collide the ego-vehicle when it is stopped. Figure 4: A scenario with an occluded area. Left: the ego vehicle follows a path (red dashed line) and an obstacle occludes a part of the road (the visible area corresponds to the transparent blue area) where a virtual vehicle has been created (drawn in purple). Right: the corresponding T P obstacle (note that its shape is the addition of the three possible cases seen in Figure 1), the candidate path s(t) and the safe path s(t) for the ego vehicle in T P .

Speed planning computation

Once the T P obstacles corresponding to the virtual and the visible vehicles and the visibility horizon line have been inserted in T P , we have to compute a candidate path s(t) for the ego vehicle. We have two cases to consider:

• If the limits of the occluded areas are sufficiently far away from the intersection with the path of the ego vehicle, the corresponding T P obstacles will be small. In this case the ego vehicle will have enough visibility to be able to cross the road junction safely and a polygonal line s(t) built as previously will pass through the left part of the zone corresponding to the road junction in T P (as in case b of Figure 5).

• When the visibility is not sufficient, our solution should be safe but not overcautious. The visibility horizon s h is moved to the smallest value of s of the first T P obstacle associated to a virtual vehicle to let it brake if another vehicle comes from the occluded area.

Then a corresponding smoothed version s(t) of s(t) is computed as in Section 3.1.3.

Experiments

Experimental evaluations have been performed using CARLA simulator [START_REF] Dosovitskiy | CARLA: an open urban driving simulator[END_REF] on an environment consisting of several road junctions, spaced 40 m apart, with a building hiding some parts of the road. (Video example at : https://team.inria.fr/rits/membres/renaud-poncelet/) All the other vehicles have been controlled by a CARLA agent. They follow a random path guided by the lanes and they drive at maximum speed v M ax except when another vehicle is in a control area of ten meters, in which case they brake. Several scenarios have been evaluated where the ego vehicle has to cross a road junction with three different behaviors :

• Omniscient: the ego vehicle knows the position and the speed of all the vehicles present in the scene, including the occluded ones.

• Naive: the ego vehicle reacts with respect to the visible vehicles and does not take into account the occluded parts of the scene.

• SGSPA for Safe Geometric Speed Planning Approach: This is the complete approach presented in this paper.

For each behavior the motion of the ego-vehicle is simulated using a longitudinal control described in Section 3.3. The following parameters have been used during the tests:

v M ax = 8.3 m/s, a M ax = 2.5 m/s 2 , a M in = -4 m/s 2 , ∆T = 5 s, t i+1 -t i = 1/5 s, dt = 1/
20 s, t gap = 2 s, margin = 1 m. The value of s h is given by the equation:

s h = 2v M ax × t h .

Crossing an X junction

The first situation is an X junction where the road leading from the right is partially occluded by a building. The ego vehicle has to cross the intersection and go straight. Two scenarios are presented:

1. Empty occluded area. Let us first consider the motion computed by SGSPA illustrated in Figure 5:

(a) When the ego vehicle is far from the intersection (t = 0.0 s), the border of the occluded area is near the intersection and the T P obstacle associated to the virtual vehicle is very long, thus the algorithm plans to go at maximum speed and to brake when the braking distance is reached. As expected, in Figure 6, we can notice that the Omniscient and the Naive approaches both plan their speed in the same way as there is no hidden vehicle.

2. One car is on the occluded part of the road. In this scenario, the other vehicle drives at v M ax , reaches the junction approximately at the same time step as the ego-vehicle and goes straight onto the junction. The ego vehicle has to give way to the the other car.

Let us go through each approach portion by portion in Figure 7: • In portion b, SGSPA begins to brake to be able to stop before the junction while Naive goes still very fast.

• In portion c, Naive brakes hard to avoid a collision with the other vehicle, which is now visible, and Omniscient begins to speed up because its speed is low enough to let the other vehicle pass.

All passed safely the intersection, but Omniscient and SGSPA approaches computed comfortable motions whereas it was not the case for Naive.

T junction

Let us now compare the three methods on a scenario inspired by [START_REF] Tas | Limited Visibility and Uncertainty Aware Motion Planning for Automated Driving[END_REF]: the ego vehicle has to turn left at a T junction where the visibility is limited. As illustrated in Figure 8, both sides are occluded: on the left hand by a building and on the right hand by a line of stopped cars on the road. The visibility is more limited on the right hand side where the arrival of a hidden vehicle enforces the ego vehicle to give way to it. In this example, the Naive approach can't brake hard enough to avoid the collision while the other approaches reduce their speeds to let the other vehicle pass, as shown in Figure 9.

Figure 9: T junction: measured speed over curvilinear abscissa for the three methods.

X junction with ten cars

Finally, to evaluate the safety of the approach in more complex cases, we created 50 scenarios for the X junction environment of Section 4.1 with ten vehicles randomly positioned on the roads and following a random path. Ten cars may not all interact with the ego-vehicle at the same time : for example, in Fig. 10, only two cars are involved.

Note that not all the vehicles are hidden by the building. In all the scenarios, the ego vehicle had to cross the X junction but with two different behaviors :

-SGSPA, where the ego vehicle is enforced to stop before the junction if it cannot pass through it.

-GSPA for Geometric Speed Planning Approach without ensuring a safe sate before intersections. Let us compare the success rate and the time spend in the junction in both case and discuss the cases where a crash had occurred:

-For SGSPA, we obtained 100 % success and the mean time spend by the ego vehicle in the junction was 3.8 s.

-For GSPA without safe state before intersections, nine collision occurred (success rate: 82 %), among which for eight of them the ego vehicle was at rest in the middle of the junction. In fact, it is not really a passive safety [START_REF] Shalev-Shwartz | On a formal model of safe and scalable self-driving cars[END_REF], as the vehicle was inside the junction. In the last case the ego-vehicle was moving but the other vehicle in collision was behind the ego-vehicle. The mean time spent inside the junction in case of success was 4.3 s.

Conclusion

In this paper we proposed a geometric approach computing a safe and comfortable speed planning for autonomous driving in urban area with occlusions. Evaluations show that SGSPA achieves both safety and comfort in all scenarios. The approach is cautious but not over defensive, thanks to the speed profile we choose when approaching an occluded region: we decide to move quickly to the intersection in order to increase the ego vehicle visibility while keeping enough space to brake. This approach can be generalized to other dynamic obstacles, such as pedestrian, considering the possible existence and intent of hidden agents as Nager et al. [START_REF] Nager | What lies in the shadows? Safe and computation-aware motion planning for autonomous vehicles using intent-aware dynamic shadow regions[END_REF].

However, in some really limited visibility cases, we cannot ensure that the ego-vehicle will not be blocked in a "freezing state" as discussed in [START_REF] Lin | Decision making through occluded intersections for autonomous driving[END_REF]. In this case, a human driver will drive very slowly in order to increase its visibility. Work in progress extends our approach to mimic such a behavior, knowing that it is not completely satisfying regarding safety and comfort.

Figure 1 :

 1 Figure 1: The three possible maneuvers for the blue vehicle coming from the right and their corresponding T P obstacles in T P .

Figure 2 :

 2 Figure 2: The different time steps used in our approach.A longitudinal motion is planned every time step t i considering the period ∆T and the control of the ego vehicle is computed every τ j

Figure

  Figure3: Computation of a trajectory curve. On the extreme left, an overview of the scene. On the middle left, the candidate path s (red line) from the visibility graph method as in[START_REF] Kant | Toward Efficient Trajectory Planning: The Path-Velocity Decomposition[END_REF] and the smoothed version (black dashed line) of s intersects one T P obstacle (in red). At t = 0 the ego-vehicle is stationnary and its speed is equal to 0. On the middle right, a new candidate path s is built, which passes after the polygon in collision. On the extreme right: the visual horizon (s h ) is virtually reduced before the intersection area (in grey). A new candidate path s is built considering a safe state before the intersection thanks to the constrain s < s brake . The last smoothed version s is safe. We notice that a solution can always be found under two conditions : The initial situation have to be safe i.e the ego-vehicle can brake without entering an intersection, and obstacles can't collide the ego-vehicle when it is stopped.

Figure 5 :

 5 Figure 5: Empty occluded area. Different time steps of SGSPA. Virtual vehicles: possible paths and T P obstacles are drawn in pink. Ego vehicle: s(t) is drawn in red and s(t) in black.

Figure 6 :

 6 Figure 6: Empty occluded area : measured speed over curvilinear abscissa for the three approaches.

Figure 7 :

 7 Figure 7: Measured speed over curvilinear abscissa for the three methods when a vehicle is occluded

Figure 8 :

 8 Figure 8: T junction with limited visibility with no consideration of the occluded parts of the scene. Naive approach does not anticipate the arrival of the hidden vehicle on the right and the collision cannot be avoided.

Figure 10 :

 10 Figure 10: X Junction with two other cars arriving at the intersection at the same time . The ego vehicle is stationary before the intersection. The blue TP corresponds to the car arriving on the left side of the ego vehicle. The TP corresponding to the real car coming from the right is masked by the virtual TP that represents the possibility of the presence of a car hidden by it on the front track.


Mohamed El

Amine Seddik

Mohamed Tamaazousti

Romain Couillet

KERNEL RANDOM MATRICES OF LARGE CONCENTRATED DATA: THE EXAMPLE OF GAN-GENERATED IMAGES

Keywords: Kernel methods, spectral clustering, random matrix theory, concentration of measure, GANs

published or not. The documents may come

INTRODUCTION

The big data paradigm involves the ability of performing classification or regression tasks on large dimensional and numerous datasets (i.e., the so-called "large p", "large n" regime of random matrix theory). Generally, the used methods for achieving these tasks are based on non-linear approaches including neural networks [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF][START_REF] Louart | A random matrix approach to neural networks[END_REF] and algorithms that are based on kernel methods, such as kernel-based support vector machines [START_REF] Liao | Random matrices meet machine learning: A large dimensional analysis of ls-svm[END_REF], semi-supervised classification [START_REF] Mai | A random matrix analysis and improvement of semi-supervised learning for large dimensional data[END_REF], kernel-based PCA [START_REF] Seddik | A kernel random matrix-based approach for sparse PCA[END_REF] and spectral clustering [START_REF] Hafiz Tiomoko | Random matrix-improved kernels for large dimensional spectral clustering[END_REF][START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF]. Due to their nonlinear design, these methods are particularly difficult to analyze theoretically. For practical large and numerous data, the study of kernel-based methods relies on the characterization of kernel matrices K ∈ R n×n in the large dimensional regime (i.e., p/n → c 0 as n → ∞). Under asymptotically non-trivial growth rate assumptions on the data statistics (i.e., maintaining a feasible get not too easy problem), the entries K ij = f (x i x j /p) or K ij = f ( x i -x j 2 /p) of K tend to a limiting constant independently of the data classes -the between and within class vectors are "equidistant" in highdimension. This observation allows one to study K through Couillet's work is supported by the GSTATS UGA IDEX Datascience chair and the ANR RMT4GRAPH (ANR-14-CE28-0006).

a Taylor expansion, thereby giving access to the characterization of functionals of K or its (informative) eigenvalueeigenvector pairs in the large dimensional regime.

Indeed, such an analysis was initiated in [START_REF] El | The spectrum of kernel random matrices[END_REF] where it has been shown that K has a linear behavior in the large p, n asymptotics. Under a k-class Gaussian mixture model, it has been shown in [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF] that the normalized Laplacian matrix associated with K behaves asymptotically as a so-called spiked random matrix, where some of the isolated eigenvalues and eigenvectors contain clustering information. In particular, the authors in [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF] demonstrated that the obtained theoretical model agrees with empirical results using the popular MNIST dataset [START_REF] Lecun | The mnist database of handwritten digits[END_REF], thereby suggesting a sort of universality of spectral clustering regarding the underlying data distribution.

The aim of this paper is to confirm this observation by relaxing the Gaussianity assumption to a wide range of distributions. In fact, most of real world data (e.g., images or CNN representations that are commonly used in computer vision [START_REF] Tamaazousti | Learning more universal representations for transfer-learning[END_REF]) belong to complex manifolds, and therefore are unlikely close to Gaussian. However, due to recent advances in generative models since the arrival of Generative Adversarial Networks [START_REF] Goodfellow | Generative adversarial nets[END_REF], it is now possible to generate complex data structures by applying successive Lipschitz operations to Gaussian vectors. On the other hand, the concentration of measure phenomenon tells us that Lipschitz transformations of Gaussian vectors satisfy a concentration property [START_REF] Tao | Topics in random matrix theory[END_REF]Thm 2.1.12]. Precisely, defining a concentrated vector X ∈ E through the real concentration of F(X), for any Lipschitz function F : E → R, defines a larger class of distributions [START_REF] Louart | Large sample covariance matrices of concentrated random vectors[END_REF]. This suggests that making the aforementioned concentration assumption on data is a suitable model for real world data.

In this paper, we analyze the kernel matrix K under a kclass concentration mixture model [START_REF] Louart | Large sample covariance matrices of concentrated random vectors[END_REF]. Precisely, we prove that K behaves (up to centering) asymptotically as a spiked random matrix in the large p large n regime, thereby generalizing the results of [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF] to a broader class of distributions. We particularly confirm our theoretical findings by considering the input data as CNN representations of images generated by a GAN, where the latter is trained to fit the manifold distribution of the well-known CIFAR-10 dataset [START_REF] Krizhevsky | Convolutional deep belief networks on cifar-10[END_REF]. We further consider real images for comparison.

Notation: Vectors are denoted by boldface lowercase letters and matrices by boldface uppercase letters. The notation • stands for the Euclidean norm for vectors and the operator norm for matrices. The vector 1 n ∈ R n denotes the vector filled with ones. For an integer k, [k] stands for the set {1, . . . , k}. [x 1 , . . . , x n ] ∈ R p×n denotes the concatenation of the vectors x 1 , . . . , x n ∈ R p . Given a normed space (E, • E ) and a real q, an element X ∈ E is said to be q-exponentially concentrated if for any 1-Lipschitz function

F : E → R, P{|F(X) -EF(X)| ≥ t} ≤ Ce -ct q
for all t > 0, and we shall write

X ∈ O(e -• q ) in (E, • E ). Remark 1. Let X ∈ O(e -• q ) in (E, • E ) and F n : E → F u n -Lipschitz
, where u n depends on some asymptotic variable n. Then, the concentration property on X is transferred to

F n (X), precisely F n (X) ∈ O(e -(•/un) q ) in (F, • F ).

MODEL SETTING AND ASSUMPTIONS

Consider n independent random vectors x 1 , . . . , x n ∈ R p distributed in k < ∞ classes represented by k distributions µ 1 , . . . , µ k supposedly all distinct. We consider the hypothesis of q-exponential concentration, meaning that there exists q ≥ 2 such that for all m ∈ N, any ∈ [k] and any family of independent vectors y 1 , . . . , y m following the distribution µ , we have the concentration

[y 1 , . . . , y m ] ∈ O(e -• q ) in (R p×m , • F ). (1)
For ∈ [k], we denote by m the mean of the distribution µ , C denotes its covariance matrix and n stands for the number of vectors among the x i 's following µ . Let m ∈ R p and C ∈ R p×p be respectively defined as

m ≡ k =1 n n m , C ≡ k =1 n n C (2)
We further denote m ≡ mm and C ≡ C -C. We shall consider the following set of assumptions on the data statistics and the kernel function in the large dimensional regime, meaning that both p and n grow at controlled joint rate. These assumptions notably guarantee the non-triviality of spectral clustering under the considered regime.

Assumption 1 (Growth rate). As n → ∞, consider the following conditions:

i. For c 0 ≡ p n ; 0 < lim inf n c 0 ≤ lim sup n c 0 < ∞. ii. For each ∈ [k], define c ≡ n n and c ≡ {c } k =1 ; 0 < lim inf n c ≤ lim sup n c < ∞. iii. lim sup p max m < ∞, lim sup p E xi √ p < ∞. iv. lim sup p max C < ∞, lim sup p max tr C √ p < ∞.
Assumption 2 (Kernel function). Let τ ≡ 2 p tr C and let f : R + → R + be a three-times continuously differentiable function in a neighborhood of the values taken by τ and such that lim inf n f (τ ) > 0.

Without loss of generality, for each ∈ [k], we arrange the x i 's as x n1+•••+n -1 +1 , . . . , x n1+•••+n ∼ µ . and define the kernel matrix K as the translation-invariant random matrix

K ≡ f 1 p x i -x j 2 n i,j=1
.

(3)

MAIN RESULTS

Our first and fundamental result states that the between and within class vectors are "equidistant" in the high-dimensional regime. Namely, we have the following lemma under the qexponential concentration hypothesis and Assumption 1.

Lemma 1. Denote τ ≡ 2 p tr C and let Assumption 1 hold. Then for any δ > 0, we have with probability at least 1 -δ

max 1≤i =j≤n 1 p x i -x j 2 -τ = O log( n √ δ ) 1/q √ n . (4)
From this observation, all the off-diagonal entries of the kernel matrix K tend to the same quantity f (τ ) asymptotically. Therefore, K can be Taylor expanded entry-wise and we show in the following that it asymptotically has (up to centering) the same behavior as a spiked random matrix.

Before introducing this asymptotic equivalent and for subsequent use, we introduce the following quantities

M = [ m1 , . . . , mk ] ∈ R p×k , t = tr C √ p k =1 ∈ R k J = [j 1 , . . . , j k ] ∈ R n×k , T = tr Ca Cb p k a,b=1 ∈ R k×k Z = [z 1 , . . . , z n ] ∈ R p×n , P = I n - 1 n 1 n 1 n ∈ R n×n
where j ∈ R n stands for the canonical vector of the class represented by µ , defined by (j ) i = δ xi∼µ . The vectors z i are defined as

z i ≡ (x i -m )/ √ p for each ∈ [k].
Our main technical result states that there exists a matrix K such that PKP -K → 0 asymptotically, where K has a tractable behavior from a mathematical standpoint.

Theorem 1 (Asymptotic Random Matrix Equivalent). Let Assumptions 1 and 2 hold and let K be defined as

K = -2f (τ ) [PZ ZP + UAU ] + F (τ )P A =    A 11 I k -f (τ ) 2f (τ ) t I k 0 k×k 0 k×1 -f (τ ) 2f (τ ) t 0 1×k -f (τ ) 2f (τ )    A 11 = M M -Ξ - f (τ ) 2f (τ ) [tt + 2T] U = J -1 n c √ p , PΦ, Pψ , Ξ = ma 2 + mb 2 2 k a,b=1 Φ = Z M -{Z m 1 k } k =1 ψ i = z i 2 -E z i 2 = z i 2 - 1 p tr C F (τ ) = (f (0) -f (τ ) + τ f (τ )).
For δ > 0, there exists C δ > 0 such that for all γ > 0

PKP -K ≤ C δ n -1/2+γ log(n) γ with proba. 1 -δ.
Theorem 1 shows that, up to centering by P, the kernel matrix K has asymptotically the same behavior as K. In particular, the obtained approximation in operator norm implies that K and K share the same eigenvalues (by Weyl's inequality [15, Thm 4.1]) and same isolated eigenvectors asymptotically. Therefore, the asymptotic spectral properties of K (i.e., the classification performance of algorithms involving K) may be studied through its equivalent K.

Indeed, note that K is made of a sum of a random matrix PZ ZP and a maximum (k -1)-rank matrix containing linear combinations of the class-wise canonical vectors j weighted by the inner-products between class means M M and class covariance-products and traces (through t and T). The matrix K can then be identified as a so-called spiked random matrix model [START_REF] Benaych | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF]. Note however that, unlike the standard spiked random matrices, the low-rank part of K depends statistically on the noise part and the latter is a mixture between random matrices made of concentrated vectors [START_REF] Louart | Large sample covariance matrices of concentrated random vectors[END_REF].

Random matrix theory offers a wide range of tools to analyze such spiked models. Importantly, authors in [START_REF] Louart | Large sample covariance matrices of concentrated random vectors[END_REF] have characterized the spectrum of a sample covariance matrix of concentrated vectors, and more precisely the bulk of eigenvalues of the random matrix PZ ZP. The spectrum of K is then composed of a bulk along with up to k -1 isolated eigenvalues, and the associated eigenvectors are aligned with the eigenvectors in U, therefore with linear combinations of the class canonical vectors j 1 , . . . , j k . Consequently, characterizing the asymptotic performance of spectral clustering relies on the characterization of the isolated eigenvectors of K. In fact, these eigenvectors are informative if their associated eigenvalues are far away from the main eigenvalue bulk.

In the following, we provide the conditions under which the informative eigenvalues become visible in the spectrum of K. Before introducing this phase transition phenomenon, let us introduce the following lemma which characterizes the spectrum of the random matrix PZ ZP.

Lemma 2 (Deterministic equivalent). Let Assumption 1 hold. Let z ∈ C

S with S introduced subsequently and define the resolvent matrix

Q δ ≡ ( k =1 c C 1+δ (z) -zI p ) -1 where δ (z) is the unique solution of the fixed point equa- tion δ (z) = 1 n tr(C Q δ ).
Then the spectral distribution ν n = 1 n n i=1 δ λi almost surely converges to a probability measure ν defined on a compact support S and having Stieltjes 1 transform the limit of 1 p tr Q δ as p → ∞.

1 Defined as m(z) = ν(dt) t-z , for z ∈ C + .
Having Lemma 2 we can now determine the conditions under which the spikes can be visible outside the main bulk of PZ ZP, and the result concerning the isolated eigenvectors. We however need the following technical assumption on the class-wise covariances to ensure that PZ ZP does not produce non-informative isolated eigenvalues. Now we can state the theorem that ensures the presence of informative eigenvalues in the spectrum of K, and gives the characterization of the corresponding isolated eigenvectors, which results from standard random matrix techniques [START_REF] Benaych | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF].

Theorem 2 (Spikes and isolated eigenvectors). Let Assumptions 1-3 hold and z ∈ C

S. There exists a matrix Λ z ∈ C k×k and a complex-valued function α τ (z) both depending on the data statistics such that, if λ * ∈ R S with α τ (λ * ) = 0 and Λ λ * having a zero eigenvalue of multiplicity m * , then PKP produces m * spikes asymptotically close to ρ * = -2f (τ )λ * + F (τ ). Furthermore, the eigenspace projector corresponding to the (asymptotically converging to ρ * ) isolated eigenvalues of PKP has a non-vanishing projection onto span(j 1 , . . . , j k ).

Theorem 2 gives the conditions under which the spikes can be observed in the spectrum of K, and states that the corresponding eigenvectors are aligned to some extent to the class canonical vectors j 1 , . . . , j k , which is important for spectral clustering. Note that, for lack of space, the explicit formulas of Λ z and α τ (z) shall be given in an extended version of the present paper. We refer the reader to [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF] for a detailed statement of Theorem 2 in the k-class Gaussian mixture model case.

APPLICATION TO GAN-GENERATED IMAGES

The k-class concentration mixture model is motivated, among other things, by the fact that data generated by GANs belong to this category of random vectors. To highlight this aspect, we evaluate our theoretical findings by considering x 1 , . . . , x n as CNN representations of GAN-generated images and we further consider real images for comparison.

GAN architecture and training: We train a conditional GAN [START_REF] Mirza | Conditional generative adversarial nets[END_REF] on the whole CIFAR-10 train set. Precisely, the generator G takes as input a Gaussian vector of dimension 100 and a one-hot-encoder vector corresponding to a given class, and outputs an image of size 32 × 32 × 3. In particular, the considered architecture for G is a deep convolutional network [START_REF] Radford | Unsupervised representation learning with deep convolutional generative adversarial networks[END_REF] composed of four convolutional layers, each one followed by batch-normalization and ReLU activation except for the last layer where tanh is used as recommended by [START_REF] Radford | Unsupervised representation learning with deep convolutional generative adversarial networks[END_REF]. Examples of the generated images are shown in top of The x i 's are therefore set as x i = N R g • G(ω i ) for the generated images, and as x i = N R r (I i ) for the real images, where the ω i 's are random Gaussian vectors in R 100 and the I i 's denote images from the CIFAR-10 test set. Importantly, note that the mapping N R g • G is Lipschitz since it is constructed from successive convolutions, activations (ReLU and tanh) and batch-normalizations, which are all Lipschitz operations [START_REF] Gouk | Regularisation of neural networks by enforcing lipschitz continuity[END_REF]. Therefore, by Remark 1, the x i 's are random concentrated vectors by design for the generated images.

In Fig. 1 we consider the spectral clustering of n = 3 × 1000 vectors x i of size p = 1024 belonging to k = 3 classes (Airplane, Automobile and Bird). In particular, we consider the clustering of these vectors at different training phases of the networks N g and N r (0, 5 and 40 epochs), with the generated images (top) and the CIFAR-10 test set images (bottom). Means and covariances for K are computed empirically and K is obtained using f (x) = exp(-x).

It is important to note, from the different histograms, that the spectrum of K is quite close to our theoretical approxima-tion K given by Theorem 1, mainly in the cases (b) and (c) for the generated images, and even for real images in (e) and (f). Another important aspect concerns the match between the spikes and the almost perfect match between the corresponding eigenvectors, which provide clustering information as predicted by Theorem 2. These observations notably show the universality aspect of spectral clustering, thereby confirming, under the concentration assumption, the observations of [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF]. Beyond this result, we have shown through this paper, for the first time, that the processing of real world data (e.g., GAN-generated images ≈ real images) can be theoretically analyzed through random matrix theory, this being made possible thanks to the concentration of measure phenomenon and the Lipschitz character of recent generative methods and their strong performances to generate complex data structures.

CONCLUSION

In this paper, we have analyzed large kernel matrices under a k-class concentration mixture model. The presented findings notably extend the results of [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF] to a wide class of distributions of the data vectors, including Lipschitz-ally transformed Gaussian vectors. Our results notably confirm a side of universality of spectral clustering as suggested in [START_REF] Couillet | Kernel spectral clustering of large dimensional data[END_REF]. More importantly, we have demonstrated, through this paper, that real data behave similar to concentrated vectors. Besides, random matrix theory allows for the theoretical understanding of machine learning methods for concentrated vectors, thereby demonstrating its relevance for machine learning.

Assumption 3 (

 3 Spikes control). Denote λ 1 , . . . , λ p the eigenvalues of C , for each ∈ [k]. As n → ∞, with support S , and max 1≤i≤p dist(λ i , S ) → 0.

Fig. 1 .

 1

Fig. 1 :

 1 Fig. 1: Spectral clustering on CNN representations of GAN-generated images (top) and CIFAR-10 images (bottom) at different training phases of the representations. The performance clustering is notably predictable through random matrix theory.

CNN representations:

  In order to build CNN representations as commonly used in computer vision[START_REF] Tamaazousti | Learning more universal representations for transfer-learning[END_REF], we train two CNNs with a 10-class classification problem. A network N g to classify a set of 50000 generated images (by G), and a network N r to classify the CIFAR-10 train set. The two networks have the same architecture: Six convolutional layers with ReLU activation, followed by a dense layer of 1024 units ReLU activated and a 10-units classification layer. CNN representations (denoted N R i for i ∈ {g, r}) correspond to the dense layer which are vectors of dimension p = 1024.