Title: A KERNEL RANDOM MATRIX-BASED APPROACH FOR SPARSE PCA

Abstract: 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.

1. 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 (Wold et al., 1987) , 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 (Anderson, 1963) 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 Ĉ. (Johnstone & Lu, 2009) 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. (Baik et al., 2005) and (Paul, 2007) 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 (Benaych-Georges & Nadakuditi, 2011; Capitaine et al., 2009; Féral & Péché, 2007; Knowles & Yin, 2013) .

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 (Krauthgamer et al., 2015) . 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. (Deshpande & Montanari, 2014; 2016) 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 .

2. 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 (d'Aspremont et al., 2005; Moghaddam et al., 2006; Zass & Shashua, 2007; Zou et al., 2006; Wright et al., 2009) , 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 (Asteris et al., 2014; Papailiopoulos et al., 2013; Shen & Huang, 2008) , 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 (Cadima & Jolliffe, 1995) . Based on the well-known power method, (Yuan & Zhang, 2013) 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), (Ma et al., 2013) proposed a very efficient iterative thresholding approach for estimating principal subspaces in the sparse setting. Similarly, assuming a single-spike model, (Krauthgamer et al., 2015) 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, (Deshpande & Montanari, 2016) analyzed and theoretically proved, under a spiked model, that indeed the covariance thresholding algorithm (Krauthgamer et al., 2015) 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 (El Karoui, 2010b; 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. (El Karoui, 2010b) 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.

3. PRELIMINARIES:

Before introducing our model setting we recall some definitions and notions of the concentration of measure theory (Ledoux, 2005 ) that are at the heart of our main results. Furthermore, we recall a definition, introduced by (El Karoui, 2008) , of sparse matrices in the large-dimensional context that will also be exploited in this paper.

3.1. 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

EQUATION

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.

3.2. ε-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 (El Karoui, 2008) 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 (El Karoui, 2008 , 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 (Bickel & Levina, 2008) . 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 (Bickel & Levina, 2008) ) provide equivalent bounds for ε < 1 2 and when considering the large dimensional p ∼ n setting (see subsection 2.4 in (Bickel & Levina, 2008) ); 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 )).

4. 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.

4.1. GENERAL SETTING AND MAIN RESULTS:

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

EQUATION

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η

EQUATION

)

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 (Tao, 2012) . 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.

EQUATION

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

EQUATION

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 (Davis & Kahan, 1970) , 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

EQUATION

4.2. 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., (Bai & Silverstein, 1998 )), 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,

EQUATION

)

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

EQUATION

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.

5. 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 (Johnstone & Lu, 2009) . 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 (Yuan & Zhang, 2013) , ITSPCA for the method in (Ma et al., 2013) , CT refers to the method in (Deshpande & Montanari, 2016 ) and finally we refer to our method as f -PCA. The total projections score 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 (Yuan & Zhang, 2013) algorithm or the ITSPCA (Ma et al., 2013) 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 (Cheng & Singer, 2013; Kammoun & Couillet, 2017) , 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 (Tiomoko Ali et al., 2018) 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.,  and F (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

EQUATION

) 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

EQUATION

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 Johnstone et al. (Johnstone & Lu, 2009) , 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 (Deshpande & Montanari, 2016) . 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", "Piece Poly" and "Step New" signals of (Johnstone & Lu, 2009) , 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 obtained 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.

:

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(Deshpande & Montanari, 2016) 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(El Karoui, 2008).

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(Bickel & Levina, 2008;El Karoui, 2008).

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)|.

Figures and Tables:

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: (Left) Spectrum of Ĉ (up) and f ( Ĉ) (bottom) for p = 2048 and n = 7500. Limiting Marčenko-Pastur density(Marčenko & Pastur, 1967) 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.

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(Johnstone & Lu, 2009). 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.

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", "Piece Poly" and "Step New" signals of(Johnstone & Lu, 2009), 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.

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.

Peer review (iclr2019):
 This paper proposes an algorithm to approximate kernel matrix based on the Taylor expansion of the element-wise functions. The authors provide a spectral norm based error bound for their method and the corresponding results for the special case, .epsilon-sparse matrix. I have some comment as follow. 1. Can you provide some comparison with Nystrom methods? It is very popular for kernel approximation and looks more efficient than the proposed algorithm. 2. The analysis relies on the Gaussian assumption on the input matrix. Can we extend it to more general case? 3. In section 5, the paper said “as our method consists in computing the sparse eigenvectors of a p .times 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”. The time complexity of the proposed algorithm is not clearly. a) Is there any bound for the sparsity level s? Why the eigenvectors of p .times p matrix is sparse? b) The convergence of power method is heavily affected by the eigen-gap of the matrix. Is there any theoretical or empirical result for the convergence behavior of power method on approximate matrix and original matrix?