A KERNEL RANDOM MATRIX-BASED APPROACH FOR SPARSE PCA
Year: 2018

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.

heading: 1 INTRODUCTION
text: 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 ∈ Rp×n consisting of n centered samples, each sample having p features. The standard PCA method requires the computation of the sample covariance matrix Ĉ = Y Y ᵀ/n and estimates the first principal components u1, u2, . . . (i.e., the successive dominant eigenvectors of C = E[Y Y ᵀ/n]) by the ordered eigenvectors û1, û2, . . . of Ĉ. (Johnstone & Lu, 2009) demonstrated that, in the high dimensional regime where n, p→∞ with p/n→ c > 0, the principal component û1 estimated by standard PCA is inconsistent. Essentially, if p/n 9 0 then ‖û1 − u1‖2 9 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 Ĉ = 1nΣ 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 = Ip + ∑k i=1 ωiuiu ᵀ 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 ûi using standard PCA is (almost surely) asymptotically orthogonal to the true principal component ui (i.e., û ᵀ i ui → 0); on the other hand, if ωi > √ c, lim infn |ûᵀi ui| > 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
∗Couillet’s work is supported by the GSTATS UGA IDEX Datascience chair and the ANR RMT4GRAPH (ANR-14-CE28-0006).
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 Ĉ induce noise in its principal components, this method consists in applying the popular soft-thresholding function (with threshold τ > 0); soft(· ; τ) : t 7→ sign(t) · (|t| − τ)+, entry-wise to the empirical covariance matrix Ĉ 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(Ĉ) instead of Ĉ 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}, bac denotes the integer part of a. Given a vector x ∈ Rn, the `2-norm of x is denoted as ‖x‖2 = ∑n i=1 x 2 i . Given an n × n matrix M , Mij 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‖2F = ∑n i,j=1M 2 ij , and the operator norm ‖M‖op of M is defined as ‖M‖op = max‖x‖=1 ‖Mx‖. Finally, denotes the Hadamard product, with [M N ]ij = MijNij .

heading: 2 RELATED WORK
text: 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
1We 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.
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/2p 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 = Ip 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.

heading: 3 PRELIMINARIES
text: 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.

heading: 3.1 CONCENTRATION OF MEASURE RESULTS
text: 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 Z2 is exp-normally concentrated, precisely
Z2 ∈ KCE ( c 2 · ) +KCN (
c 16E[Z]2 · ) , (1)
where KC > 0 is a constant depending only on C.
The extension of the notion of concentration to random vectors Z ∈ Rp demands that Rp → 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 ∈ Rp, 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 : Rp → Rq 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) = {1Mij 6=0}pi,j=1, which corresponds to a graph Gp 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 Cp(k) the set of closed walks of length k on Gp. 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 {Gp}∞p=1 satisfies, for all k ∈ 2N, |Cp(k)| ≤ Ck pε(k−1)+1 where ε ∈ [0, 1], Ck > 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 ε < 12 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) 6=x 6=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)).

heading: 4 MAIN RESULTS
text: In this section, we first present the setting of the article. Then, we provide an asymptotic equivalent to the matrix f(Ĉ). Finally, we treat as a special case the application of our result to the context of sparse PCA.

heading: 4.1 GENERAL SETTING AND MAIN RESULTS
text: Consider a data matrix Y ∈ Rp×n defined as
Y ≡ Σ1/2p X = (Ip + P )1/2X, (2)
where X ∈ Rp×n is a random matrix with i.i.d. N (0, 1) entries, P = ∑ki=1 ωiuiu ᵀ i and U = [u1, . . . , uk] ∈ Rp×k is isometric. Here, k refers to the number of principal components (or eigenvectors) u1, . . . , uk ∈ Rp to be evaluated, with ω1 > . . . > ωk > 0 the corresponding eigenvalues respectively. We define the quantity3 βp ≡ maxi ‖[Σ1/2p ]·,i‖. Assumptions: There exists B > 0 independent of p, n such that maxij |[Σ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:
2Considering the identity matrix (which is a sparse matrix), ‖Ip‖op = 1 while ‖Ip‖F = √ p→∞.
3The role of βp is to ensure the concentration of the quadratic form in equation 4 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.
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/2p ( 1 n XXᵀ − Ip ) Σ1/2p ] k .
Then for η > 0 and for an absolute constant C > 0, we have with probability at least 1− η
||F − F̃ ||op ≤ C β6p p
n3/2 √ η . (3)
For a general smooth function f , the kernel random matrix f(Ĉ) 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− Ip) 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− Ip 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− Ip. Of utmost importance to this end is the following lemma. Lemma 1 (A Concentration Result). For all i, j ∈ [p], the bilinear form gij(X) ≡ [Σ 1/2 p ]i,· ( 1 nXX ᵀ − Ip ) [Σ 1/2 p ]·,j satisfies
gij(X) ∈ KE ( c1 n β2p · ) +KN ( c2 n β4p · ) , (4)
for some absolute constants c1, c2,K > 0.
Proof. Denoting by vi the i-th column vector of the matrix Σ 1/2 p , we have by the polarization identity, for all M Hermitian, vᵀiMvj = 1 4 [(vi + vj)
ᵀM(vi + vj)− (vi − vj)ᵀM(vi − vj)]. It thus suffices to prove the result for the quadratic form g(X) = vᵀ ( 1 nXX ᵀ − Ip ) v where v ∈ Rp. Noticing that
vᵀXXᵀv = ‖vᵀX‖2 and E [ 1 nv ᵀ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 (
· 2‖v‖2
)
by Remark 1 and by Definition 2 since M 7→ vᵀM and u 7→ ‖u‖ are respectively ‖v‖-Lipschitz and 1-Lipschitz functions. We get the final result by Proposition 2.
A Taylor expansion of F around f(Σp) then leads to controlling the operator norm of f (3)(ξn) [Σ 1/2 p (XXᵀ/n − Ip)Σ1/2p ] 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 maxij |[Σp]ij | and βp) and [Σ1/2p (XXᵀ/n− Ip)Σ 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 supn maxi ωi <∞; specifically lim supn βp <∞. Under this setting, we have the following important corollary to Theorem 1.
4f and f (k) are applied entry-wise and k stands for the element-wise k-th power.
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 = Omη (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 v1, . . . , vk and ṽ1, . . . , ṽ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}
∆2i ‖vi − sṽi‖2 = Oη(n−1). (6)

heading: 4.2 SPECIAL CASE: SPARSE PCA
text: 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 |Cp(k)| ≤ p sk−1, thus P is 0-sparse in the sense of Definition 3. Similarly, Ip 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, 12 ) 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, sinceA(f (k)(Σp)) = A(Σp),6 the matrix f (k)(Σp) is also (almost) ε-sparse.
(ii) The matrix F (k) = [ Σ 1/2 p ( 1 nXX ᵀ − Ip ) Σ 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, maxi,j |F (k)ij | = Omη ( n− k 2+ 1 m ) .
Since in addition the operator norm of Σ1/2p (XXᵀ/n − Ip)Σ1/2p 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 6= 0). We precisely have the following result. Theorem 2. Let µ > 0 and suppose Σp is a 12+µ -sparse matrix. For f a three-times continuously differentiable function and for η > 0, we have for all ∈ (0, µ2(3+2µ) )
F = f(Σp) +Ob1/ cη ( n −µ 2(2+µ) + (2− 12+µ ) ) s.t. f ′(0) = f ′′(0) = 0. (7)
Proof. See Section A.3 in Appendix A. (See Corollary 1 for the notation Omη (.)).
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), F = f(Σp) +Ob1/ cη (n− 1 2+2 ).
5See Fact .1 in (El Karoui, 2008). 6Given M ∈ Rp×p, its corresponding adjacency matrix is defined as A(M) = {1Mij 6=0} p i,j=1. 7See proof of Lemma 5 in Appendix A for a proof of this result.
One may then perform a PCA on F for some function f with f ′(0) = f ′′(0) = 0 (we denote by f1,2(0) = 0 these two conditions in the following). But, while Σp = Ip+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 Ip + P is such that f1,2(0) = 0 and f(t) = t for all |t| > τ . Such a convenient choice is
f(t) = t(1− e−at2), (8) for some a > 0. This function notably satisfies
f ′(t) = 1 + e−at 2 (2at2 − 1)⇒ f ′(0) = 0, f ′′(t) = −2ate−at2 ( 2at2 − 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.

heading: 5 EXPERIMENTS
text: 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 = limn p/n.
However, thanks to sparsity, the spectrum of F = f(Ĉ) 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
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 Ĉ; namely beyond the phase transition √ p/n ≥ ω1. In fact, the standard PCA result is too noisy
compared with the one when considering f(Ĉ), 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 ûi) 2
and E = ‖UUᵀ − Û Ûᵀ‖F , where U = [u1, . . . , uk] are the ground truth principal components and Û = [û1, . . . , û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.

heading: 6 CONCLUSION
text: In 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 Ĉ and a smooth function f , we gave in this paper sufficient conditions on f to consistently estimate the principal components through the matrix f(Ĉ). 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).

heading: A PROOFS AND FURTHER EXPERIMENTS
text: 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 Johnstone et al. (Johnstone & Lu, 2009) (Section A.4).
For convenience, we make the present appendix self-contained by recalling the preliminaries and the results presented in the main paper.

heading: A.1 PRELIMINARIES
text: Proposition 2 (Square of Normally Concentrated Random Variables). Given Z ∈ CN (c ·), the random variable Z2 is exp-normally concentrated, precisely
Z2 ∈ KCE ( c 2 · ) +KCN (
c 16E[Z]2 · ) , (9)
where KC > 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 ∈ Rp, 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 : Rp → Rq for q ∈ N∗, Theorem 3 provides the normal concentration of all the random vectors F (Z). 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 graphs8 {Gp}∞p=1 satisfies, for all k ∈ 2N |Cp(k)| ≤ Ck pε(k−1)+1, (10) where ε ∈ [0, 1], Ck > 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) 6=x 6=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
8Defined through the corresponding adjacency matrix to Σp; given an p× p real symmetric matrix M , its corresponding adjacency matrix is defined as A(M) = {1Mij 6=0} p i,j=1.

heading: A.2 PROOF OF THEOREM 1
text: Setting: Consider a data matrix Y ∈ Rp×n defined as
Y ≡ Σ1/2p X = (Ip + P )1/2X, (11)
where X ∈ Rp×n with i.i.d. N (0, 1) entries, P = ∑ki=1 ωiuiu ᵀ i , U = [u1, . . . , uk] ∈ Rp×k is isometric, k refers to the number of principal components and ω1 > . . . > ωk > 0. Let βp ≡ maxi ‖[Σ1/2p ]·,i‖. Assumptions: There exists B > 0 independent of p, n such that maxij |[Σ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, we have
Theorem 1. For f a three-times continuously differentiable function, define the matrices F and F̃ respectively by
F ≡ f ( 1
n Y Y ᵀ
) = { f ([ 1
n Y Y ᵀ ] ij )}p i,j=1
F̃ ≡ f(Σp) + 2∑
k=1
f (k)(Σp) k! [ Σ1/2p ( 1 n XXᵀ − Ip ) Σ1/2p ] k .
Then for η > 0 and for an absolute constant C > 0, we have with probability at least 1− η
||F − F̃ ||op ≤ C β6p p
n3/2 √ η . (12)
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 gij(X) ≡ [Σ 1/2 p ]i,· ( 1 nXX ᵀ − Ip ) [Σ 1/2 p ]·,j satisfies
gij(X) ∈ KE ( c1 n β2p · ) +KN ( c2 n β4p · ) , (13)
for some absolute constants c1, c2,K > 0.
Proof. Denoting by vi the i-th column vector of the matrix Σ 1/2 p , we have by the polarization identity, for all M Hermitian, vᵀiMvj = 1 4 [(vi + vj)
ᵀM(vi + vj)− (vi − vj)ᵀM(vi − vj)]. It thus suffices to prove the result for the quadratic form g(X) = vᵀ ( 1 nXX ᵀ − Ip ) v where v ∈ Rp. Noticing that
vᵀXXᵀv = ‖vᵀX‖2 and E [ 1 nv ᵀ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 (
· 2‖v‖2
)
by Remark 1 and by Definition 2 since u 7→ ‖u‖ and M 7→ vᵀM are respectively 1-Lipschitz and ‖v‖-Lipschitz functions. We get the final result by Proposition 2.
Lemma 2 (A Moment Result). For gij(X) ≡ [Σ1/2p ]i,· ( 1 nXX ᵀ − Ip ) [Σ 1/2 p ]·,j , we have, for all k ∈ N and for some absolute constant Ck > 0, E|gij(X)|2k ≤ Ck
β4kp nk . (14)
Proof. Given a random variable Z, we have
∀m > 0, E|Z|m = ∫ ∞
0
mtm−1P{|Z| ≥ t}dt,
whenever the right hand side is finite. Applying this identity to the random variable gij(X) with m = 2k and exploiting the concentration property in Lemma 1 yields the result.
The proof starts by a Taylor expansion of Fij in the vicinity of [Σp]ij , i.e.,
Fij =
2∑
k=0
f (k)(σij)
k! F
(k) ij +
f (3)(ξnij)
6 F
(3) ij
where σij = [Σp]ij , ξnij ∈ b[Y Y ᵀ/n]ij , σijc,9 and F (k) is the matrix with entries
F (k) ij ≡ [Σ1/2p (n−1XXᵀ − Ip)Σ1/2p ]kij = gij(X)k.
We have by Lemma 1 that [Y Y ᵀ/n]ij concentrates around σij , so that ξnij is bounded by σij + ε, for all ε > 0, with high probability10 (note that the condition maxij |σ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 { |ξnij | ≥ σij + ε } ≤ P {|gij(X)| ≥ ε} ≤ Ke − n β2p min(c1 ε, c2 ε 2 β2p )
≤ Ke−K′ n 1 2 +2 min(c1ε,K ′c2ε 2n− 1 2 +2 ) ≡ pn → 0,
where K ′ > 0. And since f (3) is continuous, we deduce that f (3)(ξnij) is in particular bounded by
A ≡ max x∈[σij−ε,σij+ε] |f (3)(x)|,
with probability 1− pn. 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)‖2F ≤ A2‖F (3)‖2F . (15) By Lemma 2, for all k ∈ N
E‖F (k)‖2F = p∑
i,j=1
E [ |gij(X)|2k ] ≤ Ck p2β4kp nk ,
for some absolute constant Ck > 0. Thus, by Markov’s inequality, we have for all η > 0
P { ‖F (k)‖F ≥ p β2kp
n k 2 √ Ck η } ≤ η.
Recalling Eq. equation 15, we have with probability at least 1− η
‖f (3)(ξn) F (3)‖F ≤ C pβ6p
n 3 2 √ η .

heading: A.3 PROOF OF THEOREM 2

heading: Assumptions: As n→∞,
text: A1 p/n→ c ∈ (0,∞). A2 lim supn maxi ωi <∞; specifically lim supn β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 = Omη (n−α) stands for the fact that P { ‖X‖op ≥ C n−α η− 1 2m } ≤ η for
some absolute constant C > 0 and non-negative integer m. 9The notation ba, bc stands for the interval [a, b] if a < b or [b, a] otherwise.
10For 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).
Theorem 2. Let µ > 0 and suppose Σp is a 12+µ -sparse matrix. For f a three-times continuously differentiable function and for η > 0, we have for all ∈ (0, µ2(3+2µ) )
F = f(Σp) +Ob1/ cη ( n −µ 2(2+µ) + (2− 12+µ ) ) s.t. f ′(0) = f ′′(0) = 0. (17)
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 notion11
Lemma 3. Given an ε-sparse p×p real symmetric matrix M and calling m = maxij |Mij |, we have, for all k ∈ 2N ‖M‖op ≤ trace(Mk)1/k = O(mpε(1−1/k)+1/k). (18) Lemma 4. Given two real symmetric matrices M and N with |Mij | ≤ Nij . Then, we have ‖M‖op ≤ ‖N‖op. 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 12+µ -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+µ)−22(3+2µ) )
‖f (k)(Σp) F (k)‖op = Ob1/ cη ( n 2−k(2+µ) 2(2+µ) + (2− 12+µ ) ) .
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 (k)ij |2m = E|gij(X)|2km = O(n−km),
thus applying Markov’s inequality to the random variable |F (k)ij |2m yields to the following tail control.
P{|F (k)ij | ≥ t} ≤ E|F (k)ij |2m
t2m ≤ C n−km t−2m,
where C is an absolute constant. Recalling Assumption A1 and by the union bound, we have
P{max ij |F (k)ij | ≥ t} ≤
p∑
i,j=1
P{|F (k)ij | ≥ t} ≤ p2 P{|F (k) ij | ≥ t} ≤ C n2−km t−2m,
which implies for η > 0 and for all m > 0
max ij |F (k)ij | = Omη
( n− k 2+ 1 m ) (19)
Besides, let M be the matrix defined as M ≡ maxij |F (k)ij | · f (k)(Σp), we have
|[f (k)(Σp) F (k)]ij | ≤Mij , thus, one has by Lemma 4
‖f (k)(Σp) F (k)‖op ≤ ‖M‖op = max ij |F (k)ij | · ‖f (k)(Σp)‖op.
In particular, since f (k)(Σp) is 12+µ -sparse (by Remark 2), we have by Lemma 3 and by equation 19, for some η > 0
‖f (k)(Σp) F (k)‖op = O2mη ( n 1 2+µ (1− 12m )+ 12m− k2+ 12m ) ,
choosing = 12m < 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.
11Through the identity trace(Mk) ≤ maxij |Mij |k · |Cp(k)|.

heading: A.4 FURTHER EXPERIMENTS

heading: A.4.1 HIGHER RANK CASE
text: 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.
A.4.2 OTHER CHOICES OF THE KERNEL FUNCTION f
In this section, we consider functions of the form f(t) = αt3 +βt2 + γ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 α 6= 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) = αt3 with α 6= 0 is not optimal as the obtained signal is deformed (due to the unverified linearity condition), compared to the GT one.
α = 0, β = 0, γ 6= 0
α = 0, β 6= 0, γ = 0
α = 0, β 6= 0, γ 6= 0

Review summaries:
This paper discusses the sparse PCA problem from the random matrix perspective. First, it establishes a theorem that connect a f(sample covariance) to f(actual covariance) in Theorem 1, and shows that when f if three-times continuous differentiable,  f(sample covariance) can be written in terms of actual covariance and $1/n XX^T-I$. Then, based on this theorem, it shows that if f'(0)=f''(0)=0, then f(sample covariance) is well approximated by f(actual covariance).