GENERATIVE PRINCIPAL COMPONENT ANALYSIS
Year: 2022

Abstract:
In this paper, we study the problem of principal component analysis with generative modeling assumptions, adopting a general model for the observed matrix that encompasses notable special cases, including spiked matrix recovery and phase retrieval. The key assumption is that the first principal eigenvector lies near the range of anL-Lipschitz continuous generative model with bounded k-dimensional inputs. We propose a quadratic estimator, and show that it enjoys a statistical rate of order √ k logL m , where m is the number of samples. Moreover, we provide a variant of the classic power method, which projects the calculated data onto the range of the generative model during each iteration. We show that under suitable conditions, this method converges exponentially fast to a point achieving the above-mentioned statistical rate. This rate is conjectured in (Aubin et al., 2019; Cocola et al., 2020) to be the best possible even when we only restrict to the special case of spiked matrix models. We perform experiments on various image datasets for spiked matrix and phase retrieval models, and illustrate performance gains of our method to the classic power method and the truncated power method devised for sparse principal component analysis.

text: of order √
k logL m , where m is the number of samples. Moreover, we provide a
variant of the classic power method, which projects the calculated data onto the range of the generative model during each iteration. We show that under suitable conditions, this method converges exponentially fast to a point achieving the above-mentioned statistical rate. This rate is conjectured in (Aubin et al., 2019; Cocola et al., 2020) to be the best possible even when we only restrict to the special case of spiked matrix models. We perform experiments on various image datasets for spiked matrix and phase retrieval models, and illustrate performance gains of our method to the classic power method and the truncated power method devised for sparse principal component analysis.

heading: 1 INTRODUCTION
text: Principal component analysis (PCA) is one of the most popular techniques for data processing and dimensionality reduction (Jolliffe, 1986), with an abundance of applications such as image recognition (Hancock et al., 1996), gene expression data analysis (Alter et al., 2000), and clustering (Ding & He, 2004; Liu & Tan, 2019). PCA seeks to find the directions that capture maximal variances in vector-valued data. In more detail, letting x1,x2, . . . ,xm be m realizations of a random vector x ∈ Rn with a population covariance matrix Σ̄ ∈ Rn×n, PCA aims to reconstruct the top principal eigenvectors of Σ̄. The first principal eigenvector can be computed as follows:
u1 = arg max w∈Rn
wTΣw s.t. ‖w‖2 = 1, (1)
where the empirical covariance matrix is defined as Σ := 1m ∑m i=1(xi − c)(xi − c)T , with c := 1 m ∑m i=1 xi. In addition, subsequent principal eigenvectors can be estimated by similar optimization problems subject to being orthogonal to the previous vectors.
PCA is consistent in the conventional setting where the dimension of the data n is relatively small compared to the sample size m (Anderson, 1962), but leads to rather poor estimates in the highdimensional setting wherem n. In particular, it has been shown in various papers that the empirical principal eigenvectors are no longer consistent estimates of their population counterparts (Nadler, ∗Corresponding authors.
2008; Johnstone & Lu, 2009; Jung & Marron, 2009; Birnbaum et al., 2013). In order to tackle the curse of dimensionality, a natural approach is to impose certain structural constraints on the principal eigenvectors. A common assumption is that the principal eigenvectors are sparse, and this gives rise to the problem of sparse principal component analysis (SPCA) (Zou et al., 2006). In particular, for recovering the top principal eigenvector, the optimization problem of SPCA is given by
u1 = arg max w∈Rn
wTΣw s.t. ‖w‖2 = 1, ‖w‖0 ≤ s, (2)
where ‖w‖0 = |{i : wi 6= 0}| denotes the number of non-zero entries of w, and s ∈ N represents the sparsity level. In addition to reducing the effective number of parameters, the sparsity assumption also enhances the interpretability (Zou et al., 2006).
Departing momentarily from the PCA problem, recent years have seen tremendous advances in deep generative models in a wide variety of real-world applications (Foster, 2019). This has motivated a new perspective of the related problem of compressed sensing (CS), in which the standard sparsity assumption is replaced by a generative modeling assumption. That is, the underlying signal is assumed to lie near the range of a (deep) generative model (Bora et al., 2017). The authors of (Bora et al., 2017) characterized the number of samples required to attain an accurate reconstruction, and also presented numerical results on image datasets showing that compared to sparsity-based methods, generative priors can lead to large reductions (e.g., a factor of 5 to 10) in the number of measurements needed to recover the signal up to a given accuracy. Additional numerical and theoretical results concerning inverse problems using generative models have been provided in (Van Veen et al., 2018; Dhar et al., 2018; Heckel & Hand, 2019; Jalal et al., 2020; Liu & Scarlett, 2020a; Ongie et al., 2020; Whang et al., 2020; Jalal et al., 2021; Nguyen et al., 2021), among others.
In this paper, following the developments in both PCA/SPCA and inverse problems with generative priors, we study the use of generative priors in principal component analysis (GPCA), which gives a generative counterpart of SPCA in (2), formulated as follows:
u1 = arg max w∈Rn
wTΣw s.t. w ∈ Range(G), (3)
where G is a (pre-trained) generative model, which we assume has a range contained in the unit sphere of Rn.1 Similarly to SPCA, the motivation for this problem is to incorporate prior knowledge on the vector being recovered (or alternatively, a prior preference), and to permit meaningful recovery and theoretical bounds even in the high-dimensional regime m n.

heading: 1.1 RELATED WORK
text: In this subsection, we summarize some relevant works, which can roughly be divided into (i) the SPCA problem, and (ii) signal recovery with generative models.
SPCA: It has been proved that the solution of the SPCA problem in (2) attains the optimal statistical rate √ s log n/m (Vu & Lei, 2012), where m is the number of samples, n is the ambient dimension, and s is the sparsity level of the first principal eigenvector. However, due to the combinatorial constraint, the computation of (2) is intractable. To address this computational issue, in recent years, an extensive body of practical approaches for estimating sparse principal eigenvectors have been proposed in the literature, including (d’Aspremont et al., 2007; Vu et al., 2013; Chang et al., 2016; Moghaddam et al., 2006; d’Aspremont et al., 2008; Jolliffe et al., 2003; Zou et al., 2006; Shen & Huang, 2008; Journée et al., 2010; Hein & Bühler, 2010; Kuleshov, 2013; Yuan & Zhang, 2013; Asteris et al., 2011; Papailiopoulos et al., 2013), just to name a few.
Notably, statistical guarantees for several approaches have been provided. The authors of (Yuan & Zhang, 2013) propose the truncated power method (TPower), which adds a truncation operation to the power method to ensure the desired level of sparsity. It is shown that this approach attains the optimal statistical rate under appropriate initialization. Most approaches for SPCA only focus on estimating the first principal eigenvector, with a certain deflation method (Mackey, 2008) being leveraged to reconstruct the rest. However, there are some exceptions; for instance, an iterative thresholding approach is proposed in (Ma, 2013), and is shown to attain a near-optimal statistical rate
1Similarly to (Liu et al., 2020; 2021a), we assume that the range of G is contained in the unit sphere for convenience. Our results readily transfer to general (unnormalized) generative models by considering its normalized version. See Remark 1 for a detailed discussion.
when estimating multiple individual principal eigenvectors. In addition, the authors of (Cai et al., 2013) propose a regression-type method that gives an optimal principal subspace estimator. Both works (Ma, 2013; Cai et al., 2013) rely on the assumption of a spiked covariance model to ensure a good initial vector. To avoid the spiked covariance model assumption, the work (Wang et al., 2014) proposes a two-stage procedure that attains the optimal subspace estimator in polynomial time.
Signal recovery with generative models: Since the seminal work (Bora et al., 2017), there has been a substantial volume of papers studying various inverse problems with generative priors. One of the problems more closely related to PCA is spectral initialization in phase retrieval, which amounts to solving an eigenvalue problem. Phase retrieval with generative priors has been studied in (Hand et al., 2018; Hyder et al., 2019; Jagatap & Hegde, 2019; Wei et al., 2019; Shamshad & Ahmed, 2020; Aubin et al., 2020; Liu et al., 2021b). In particular, the work (Hand et al., 2018) models the underlying signal as being in the range of a fully-connected ReLU neural network with no offsets, and all the weight matrices of the ReLU neural network are assumed to have i.i.d. zero-mean Gaussian entries. In addition, the neural network needs to be sufficiently expansive in the sense that ni ≥ Ω(ni−1 log ni−1), where ni is the width of the i-th layer. Under these assumptions, the authors establish favorable global optimization landscapes for the corresponding objective, and derive a near-optimal sample complexity upper bound. They minimize the objective function directly over the latent space in Rk using gradient descent, which may suffer from local minima in general optimization landscapes (Hyder et al., 2019; Shah & Hegde, 2018).
In (Aubin et al., 2020), the assumptions on the neural network are similar to those in (Hand et al., 2018), relaxing to general activation functions (beyond ReLU) and ni ≥ Ω(ni−1). The authors focus on the high dimensional regime where n,m, k →∞with the ratiom/n being fixed, and assume that the input vector in Rk is drawn from a separable distribution. They derive sharp asymptotics for the information-theoretically optimal performance and for the associated approximate message passing (AMP) algorithm. Both works (Hand et al., 2018; Aubin et al., 2020) focus on noiseless phase retrieval. When only making the much milder assumption that the generative model is Lipschitz continuous, with no assumption on expansiveness, Gaussianity, and offsets, a spectral initialization step (similar to that of sparse phase retrieval) is typically required in order to accurately reconstruct the signal (Netrapalli et al., 2015; Candès et al., 2015). The authors of (Liu et al., 2021b) propose an optimization problem similar to (3) for the spectral initialization for phase retrieval with generative models. It was left open in (Liu et al., 2021b) how to solve (or approximate sufficiently accurately) the optimization problem in practice.
Understanding the eigenvalues of spiked random matrix models has been a central problem of random matrix theory, and spiked matrices have been widely used in the statistical analysis of SPCA. Recently, theoretical guarantees concerning spiked matrix models with generative priors have been provided in (Aubin et al., 2019; Cocola et al., 2020). In particular, in (Aubin et al., 2019), the assumptions are similar to those in (Aubin et al., 2020), except that the neural network is assumed to have exactly one hidden layer. The Bayes-optimal performance is analyzed, and it is shown that the AMP algorithm can attain this optimal performance. In addition, the authors of (Aubin et al., 2019) propose the linearized approximate message passing (LAMP) algorithm, which is a spectral algorithm specifically designed for single-layer feedforward neural networks with no bias terms. The authors show its superiority to classical PCA via numerical results on the Fashion-MNIST dataset. In (Cocola et al., 2020), the same assumptions are made as those in (Hand et al., 2018) on the neural network, and the authors demonstrate the benign global geometry for a nonlinear least squares objective. Similarly to (Hand et al., 2018), the objective is minimized over Rk using a gradient descent algorithm, which can get stuck in local minima for general global geometries.

heading: 1.2 CONTRIBUTIONS
text: The main contributions of this paper are as follows:
• We study eigenvalue problems with generative priors (including GPCA), and characterize the statistical rate of a quadratic estimator similar to (3) under suitable assumptions.
• We propose a variant of the classic power method, which uses an additional projection operation to ensure that the output of each iteration lies in the range of a generative model. We refer to our method as projected power method (PPower). We further show that under appropriate conditions (most notably, assuming exact projections are possible), PPower
obtains a solution achieving a statistical rate that is conjectured to be optimal in (Aubin et al., 2019; Cocola et al., 2020) for spiked matrix models.
• For the spiked matrix and phase retrieval models, we perform numerical experiments on image datasets, and demonstrate that when the number of samples is relatively small compared to the ambient dimension, PPower leads to significantly better performance compared to the classic power method and TPower.
Compared to the above-mentioned works that use generative models, we make no assumption on expansiveness, Gaussianity, and offsets for the generative model, and we consider a data model that simultaneously encompasses both spiked matrix and phase retrieval models, among others.

heading: 1.3 NOTATION
text: We use upper and lower case boldface letters to denote matrices and vectors respectively. We write [N ] = {1, 2, · · · , N} for a positive integerN , and we use IN to denote the identity matrix in RN×N . A generative model is a function G : D → Rn, with latent dimension k, ambient dimension n, and input domain D ⊆ Rk. We focus on the setting where k n. For a set S ⊆ Rk and a generative model G : Rk → Rn, we write G(S) = {G(z) : z ∈ S}. We use ‖X‖2→2 to denote the spectral norm of a matrix X. We define the `q-ball Bkq (r) := {z ∈ Rk : ‖z‖q ≤ r} for q ∈ [0,+∞]. Sn−1 := {x ∈ Rn : ‖x‖2 = 1} represents the unit sphere in Rn. The symbols C,C ′, C ′′ are absolute constants whose values may differ from line to line.

heading: 2 PROBLEM SETUP
text: In this section, we formally introduce the problem, and overview some important assumptions that we adopt. Except where stated otherwise, we will focus on the following setting:
• We have a matrix V ∈ Rn×n satisfying V = V̄ + E, (4)
where E is a perturbation matrix, and V̄ is assumed to be positive semidefinite (PSD). For PCA and its constrained variants, V and V̄ can be thought of as the empirical and population covariance matrices, respectively.
• We have an L-Lipschitz continuous generative model G : Bk2 (r) → Rn. For convenience, similarly to that in (Liu et al., 2020), we assume that Range(G) ⊆ Sn−1. Remark 1. For a general (unnormalized) L-Lipschitz continuous generative modelG, we can instead consider a corresponding normalized generative model G̃ : D → Sn−1 as in (Liu et al., 2021b), where D := {z ∈ Bk2 (r) : ‖G(z)‖2 > Rmin} for some Rmin > 0, and G̃(z) = G(z)‖G(z)‖2 . Then, the Lipschitz constant of G̃ becomes L/Rmin. For a d-layer neural
network, we typically have L = nΘ(d) (Bora et al., 2017). Thus, we can set Rmin to be as small as 1/nΘ(d) without changing the scaling laws, which makes the dependence on Rmin very mild.
• We aim to solve the following eigenvalue problem with a generative prior:2
v̂ := max w∈Rn
wTVw s.t. w ∈ Range(G). (5)
Note that since Range(G) ⊆ Sn−1, we do not need to impose the constraint ‖w‖2 = 1. Since V is not restricted to being an empirical covariance matrix, (5) is more general than GPCA in (3). However, we slightly abuse terminology and also refer to (5) as GPCA.
• To approximately solve (5), we use a projected power method (PPower), which is described by the following iterative procedure:3
w(t+1) = PG ( Vw(t) ) , (6)
2To find the top r rather than top one principal eigenvectors that are in the range of a generative model, we may follow the common approach to use the iterative deflation method for PCA/SPCA: Subsequent principal eigenvectors are derived by recursively removing the contribution of the principal eigenvectors that are calculated already under the generative model constraint. See for example (Mackey, 2008).
3In similar iterative procedures, some works have proposed to replace V by V + ρIn for some ρ ∈ R to improve convergence, e.g., see Deshpande et al. (2014).
Algorithm 1 A projected power method for GPCA (PPower) Input: V, number of iterations T , pre-trained generative model G, initial vector w(0) Procedure: Iterate w(t+1) = PG ( Vw(t) ) for t ∈ {0, 1, . . . , T − 1}, and return w(T )
where PG(·) is the projection function onto G(Bk2 (r)),4 and the initialization vector w(0) may be chosen either manually or randomly, e.g., uniform over Sn−1. Often the initialization vector w(0) plays an important role and we may need a careful design for it. For example, for phase retrieval with generative models, as mentioned in (Liu et al., 2021b, Section V), we may choose the column corresponding to the largest diagonal entry of V as the starting point. See also (Yuan & Zhang, 2013, Section 3) for a discussion on the initialization strategy for TPower devised for SPCA. We present the algorithm corresponding to (6) in Algorithm 1. Remark 2. To tackle generalized eigenvalue problems encountered in some specific applications, there are variants of the projected power method, which combine a certain power iteration with additional operations to ensure sparsity or enforce other constraints. These applications include but not limited to sparse PCA (Journée et al., 2010; Yuan & Zhang, 2013), phase synchronization (Boumal, 2016; Liu et al., 2017), the hidden clique problem (Deshpande & Montanari, 2015), the joint alignment problem (Chen & Candès, 2018), and coneconstrained PCA (Deshpande et al., 2014; Yi & Neykov, 2020). For example, under the simple spiked Wigner model (Perry et al., 2018) for the observed data matrix V with the underlying signal being assumed to lie in a convex cone, the authors of (Deshpande et al., 2014) show that cone-constrained PCA can be computed efficiently via a generalized projected power method. In general, the range of a Lipschitz-continuous generative model is non-convex and not a cone. In addition, we consider a matrix model that is more general than the spiked Wigner model.
Although it is not needed for our main results, we first state a lemma (proved in Appendix A) that establishes a monotonicity property with minimal assumptions, only requiring that V is PSD; see also Proposition 3 of Yuan & Zhang (2013) for an analog in sparse PCA. By comparison, our main results in Section 4 will make more assumptions, but will also provide stronger guarantees. Note that the PSD assumption holds, for example, when E = 0, or when V is a sample covariance matrix. Lemma 1. For any x ∈ Rn, let Q(x) = xTVx. Then, if V is PSD, the sequence {Q(w(t))}t>0 for w(t) in (6) is monotonically non-decreasing.

heading: 3 SPECIALIZED DATA MODELS AND EXAMPLES
text: In this section, we make more specific assumptions on V = V̄ + E, starting with the following. Assumption 1 (Assumption on V̄). Assume that V̄ is PSD with eigenvalues λ̄1 > λ̄2 ≥ . . . ≥ λ̄n ≥ 0. We use x̄ (a unit vector) to represent the eigenvector of V̄ that corresponds to λ̄1.
In the following, it is useful to think of x̄ is being close to the range of the generative modelG. In the special case of (3), lettingm be the number of samples, it is natural to derive that the upper bound of ‖E‖2→2 grows linearly in (n/m)b for some positive constant b such as 12 or 1 (with high probability; see, e.g., (Vershynin, 2010, Corollary 5.35)). In the following, we consider general scenarios with V depending on m samples (see below for specific examples). Similarly to (Yuan & Zhang, 2013), we may consider a restricted version of ‖E‖2→2, leading to the following. Assumption 2 (Assumption on E). Let S1, S2 be two (arbitrary) finite sets in Rn satisfying m = Ω(log(|S1| · |S2|)). Then, we have for all s1 ∈ S1 and s2 ∈ S2 that∣∣sT1 Es2∣∣ ≤ C √ log(|S1| · |S2|)
m · ‖s1‖2 · ‖s2‖2, (7)
where C is an absolute constant. In addition, we have ‖E‖2→2 = O(n/m).5
4That is, for any x ∈ Rn, PG(x) := arg minw∈Range(G) ‖w − x‖2. We will implicitly assume that the projection step can be performed accurately, e.g., (Deshpande et al., 2014; Shah & Hegde, 2018; Peng et al., 2020), though in practice approximate methods might be needed, e.g., via gradient descent (Shah & Hegde, 2018) or GAN-based projection methods (Raj et al., 2019).
5For the spectral norm of E, one often expects an even tighter bound O( √ n/m), but we use O(n/m) to simplify the analysis of our examples. Moreover, at least under the typical scaling where L is polynomial
The following examples show that when the number of measurements is sufficiently large, the data matrices corresponding to certain spiked matrix and phase retrieval models satisfy the above assumptions with high probability. Short proofs are given in Appendix B for completeness. Example 1 (Spiked covariance model). In the spiked covariance model (Johnstone & Lu, 2009; Deshpande & Montanari, 2016), the observed vectors x1,x2, . . . ,xm ∈ Rn are of the form
xi = r∑ q=1 √ βquq,isq + zi, (8)
where s1, . . . , sr ∈ Rn are orthonormal vectors that we want to estimate, while zi ∼ N (0, In) and uq,i ∼ N (0, 1) are independent and identically distributed. In addition, β1, . . . , βr are positive constants that dictate the signal-to-noise ratio (SNR). To simplify the exposition, we focus on the rank-one case and drop the subscript q ∈ [r]. Let
V = 1
m m∑ i=1 (xix T i − In), (9)
and V̄ = E[V] = βssT .6 Then, V̄ satisfies Assumption 1 with λ̄1 = β > 0, λ̄2 = . . . = λ̄n = 0, and x̄ = s. In addition, letting E = V − V̄, the Bernstein-type inequality (Vershynin, 2010, Proposition 5.10) for the sum of sub-exponential random variables yields that for any finite sets S1, S2 ⊂ Rn, when m = Ω ( log(|S1| · |S2|)), with probability 1− e−Ω(log(|S1|·|S2|)), E satisfies (7) in Assumption 2. Moreover, standard concentration arguments give ‖E‖2→2 = O(n/m) with probability 1− e−Ω(n). Remark 3. We can also consider the simpler spiked Wigner model (Perry et al., 2018; Chung & Lee, 2019) where V = βssT + 1√
n H, with the signal s being a unit vector, β > 0 being an SNR
parameter, and H ∈ Rn×n being a symmetric matrix with entries drawn i.i.d. (up to symmetry) from N (0, 1). In this case, when m = n is sufficiently large, with high probability, V̄ := E[V] = βssT and E := V − V̄ similarly satisfy Assumptions 1 and 2 respectively. Example 2 (Phase retrieval). Let A ∈ Rm×n be a matrix having i.i.d. N (0, 1) entries, and let aTi be the i-th row of A. For some unit vector s, suppose that the observed vector is y = |As|, where the absolute value is applied element-wise.7 We construct the weighted empirical covariance matrix as follows (Zhang et al., 2017; Liu et al., 2021b):
V = 1
m m∑ i=1 ( yiaia T i 1{l<yi<u} − γIn ) , (10)
where u > l > 1 are positive constants, and for g ∼ N (0, 1), γ := E [ |g|1{l<|g|<u} ] . Let V̄ =
E[V] = βssT , where β := E [( |g|3 − |g| ) 1{l<|g|<u} ] . Then, V̄ satisfies Assumption 1 with λ̄1 = β > 0, λ̄2 = . . . = λ̄n = 0, and x̄ = s. In addition, letting E = V − V̄, we have similarly to Example 1 that E satisfies Assumption 2 with high probability.

heading: 4 MAIN RESULTS
text: The following theorem concerns globally optimal solutions of (5). The proof is given in Appendix D. Theorem 1. Let V = V̄+E with Assumptions 1 and 2 being satisfied by V̄ and E respectively, and let xG := PG(x̄) = arg minw∈Range(G) ‖w − x̄‖2. Suppose that v̂ is a globally optimal solution to (5). Then, for any δ ∈ (0, 1), we have
‖v̂v̂T − x̄x̄T ‖F = O
(√ k log Lrδ m ) λ̄1 − λ̄2 +O (√ δn/m λ̄1 − λ̄2 ) +O √ (λ̄1 + n)‖x̄− xG‖2 λ̄1 − λ̄2  , (11) where n = O (√k log Lrδ m ) .
in n (Bora et al., 2017), the upper bound for ‖E‖2→2 can be easily relaxed to O ( (n/m)b ) for any positive constant b, without affecting the scaling of our derived statistical rate. 6To avoid non-essential complications, β is typically assumed to be known (Johnstone & Lu, 2009). 7Without loss of generality, we assume that s is a unit vector. For a general signal s, we may instead focus
on estimating s̄ = s/‖s‖2, and simply use 1m ∑m i=1 yi to approximate ‖s‖2.
We have stated this result as an upper bound on ‖v̂v̂T − x̄x̄T ‖F, which intuitively measures the distance between the 1D subspaces spanned by v̂ and x̄. Note, however, that for any two unit vectors w1,w2 with wT1 w2 ≥ 0, the distances ‖w1−w2‖2 and ‖w1wT1 −w2wT2 ‖F are equivalent up to constant factors, whereas if wT1 w2 ≤ 0, then a similar statement holds for ‖w1 + w2‖2. See Appendix C for the precise statements.
In Theorem 1, the final term quantifies the effect of representation error. When there is no such error (i.e., x̄ ∈ Range(G)), under the scaling λ̄1 − λ̄2 = Θ(1), L = nΩ(1), and δ = O(1/n), Theorem 1 simplifies to ‖v̂v̂T − x̄x̄T ‖F = O( √ k logL m ). This provides a natural counterpart to the statistical
rate of order √
s logn m for SPCA mentioned in Section 1.1.
Before providing our main theorem for PPower described in (6), we present the following important lemma, whose proof is presented in Appendix E. To simplify the statement of the lemma, we fix δ to be O(1/n), though a more general form analogous to Theorem 1 is also possible. Lemma 2. Let V = V̄ + E with Assumptions 1 and 2 being satisfied by V̄ and E respectively, and further assume that x̄ ∈ Range(G). Let γ̄ = λ̄2/λ̄1 with λ̄1 = Θ(1). Then, for all s ∈ Range(G) satisfying sT x̄ > 0, we have
‖PG(Vs)− x̄‖2 ≤ 2γ̄‖s− x̄‖2
sT x̄ +
C
sT x̄
√ k log(nLr)
m , (12)
where C is an absolute constant. Remark 4. The assumption sT x̄ > 0 will be particularly satisfied when the range of G only contains nonnegative vectors. As mentioned in various works studying nonnegative SPCA (Zass & Shashua, 2007; Sigg & Buhmann, 2008; Asteris et al., 2014), for several practical fields such as economics, bioinformatics, and computer vision, it is natural to assume that the underlying signal has no negative entries. More generally, the assumption sT x̄ > 0 can be removed if we additionally have that −x̄ is also contained in the range of G. For this case, when sT x̄ < 0, we can instead derive an upper bound for ‖ŝ + x̄‖2.
Based on Lemma 2, we have the following theorem, whose proof is given in Appendix F. Theorem 2. Suppose that the assumptions on the data model V = V̄ + E are the same as those in Lemma 2, and assume that there exists t0 ∈ N such that x̄Tw(t0) = 2γ̄+ν with 2γ̄+ν ≤ 1−τ , where γ̄ = λ̄2/λ̄1 ∈ [0, 1), and ν, τ are both positive and scale as Θ(1). Let µ0 = 2γ̄x̄Tw(t0) = 2γ̄ 2γ̄+ν < 1, and in addition, suppose that m ≥ Cν,τ · k log(nLr) with Cν,τ > 0 being large enough. Then, we have after ∆0 = O ( log ( m k log(nLr) )) iterations of PPower (beyond t0) that
‖w(t) − x̄‖2 ≤ C
(1− µ0)ν
√ k log(nLr)
m , (13)
i.e., this equation holds for all t ≥ T0 := t0 + ∆0. Moreover, if γ̄ = 0 then ∆0 ≤ 1, whereas if γ̄ = Θ(1), we have exponentially fast convergence via the following contraction property: There exists a constant ξ ∈ (0, 1) such that for t ∈ [t0, T0), it holds that
‖w(t+1) − x̄‖2 ≤ (1− ξ)‖w(t) − x̄‖2. (14)
Regarding the assumption x̄Tw(t0) ≥ 2γ̄ + ν, we note that when t0 = 0, this condition can be viewed as having a good initialization. For both Examples 1 and 2, we have γ̄ = 0. Thus, for the spiked covariance and phase retrieval models corresponding to these examples, the assumption x̄Tw(t0) ≥ 2γ̄+ν reduces to x̄Tw(t0) ≥ ν for a sufficiently small positive constant ν, which results in a mild assumption. Such an assumption is also required for the projected power method devised for cone-constrained PCA under the simple spiked Wigner model, with the underlying signal being assumed to lie in a convex cone; see (Deshpande et al., 2014, Theorem 3). Despite using a similar assumption on the initialization, our proof techniques are significantly different from Deshpande et al. (2014); see Appendix G for discussion.
When L is polynomial in n, Theorem 2 reveals that we have established conditions under which PPower in (6) converges exponentially fast to a point achieving the statistical rate of order √
k logL m .
Based on the minimax rates for SPCA (Vu & Lei, 2012; Birnbaum et al., 2013) and the informationtheoretic lower bounds for CS with generative models (Liu & Scarlett, 2020b; Kamath et al., 2020), the optimal rate for GPCA is naturally conjectured to be of the same order √
k logL m . We highlight
that Theorem 2 partially addresses the computational-to-statistical gap (e.g., see (Wang et al., 2016; Hand et al., 2018; Aubin et al., 2019; Cocola et al., 2020)) for spiked matrix recovery and phase retrieval under a generative prior, though closing it completely would require efficiently finding a good initialization and addressing the assumption of exact projections.
Perhaps the main caveat to Theorem 2 is that it assumes the projection step can be performed exactly. However, this is a standard assumption in analyses of projected gradient methods, e.g., see (Shah & Hegde, 2018), and both gradient-based projection and GAN-based projection have been shown to be highly effective in practice (Shah & Hegde, 2018; Raj et al., 2019).

heading: 5 EXPERIMENTS
text: In this section, we experimentally study the performance of Algorithm 1 (PPower). We note that these experiments are intended as a simple proof of concept rather than seeking to be comprehensive, as our contributions are primarily theoretical. We compare with the truncated power method (TPower) devised for SPCA proposed in (Yuan & Zhang, 2013, Algorithm 1) and the vanilla power method (Power) that performs the iterative procedure w(t+1) = (Vw(t))/‖Vw(t)‖2. For a fair comparison, for PPower, TPower, and Power, we use the same initial vector. Specifically, as mentioned in (Liu et al., 2021b, Section V), we choose the initialization vector w(0) as the column of V that corresponds to its largest diagonal entry. For all three algorithms, the total number of iterations T is set to be 30. To compare the performance across algorithms, we use the scale-invariant
Cosine Similarity metric defined as Cossim ( x̄,w(T ) ) :=
〈x̄,w(T )〉 ‖x̄‖2‖w(T )‖2 , where x̄ is the ground-truth
signal to estimate, and w(T ) denotes the output vector of the algorithm.
The experiments are performed on the MNIST (LeCun et al., 1998), Fashion-MNIST (Xiao et al., 2017) and CelebA (Liu et al., 2015) datasets, with the numerical results for the Fashion-MNIST and CelebA datasets being presented in Appendix H and I. The MNIST dataset consists of 60, 000 images of handwritten digits. The size of each image is 28 × 28, and thus n = 784. To reduce the impact of local minima, we perform 10 random restarts, and choose the best among these. The cosine similarity is averaged over the test images, and also over these 10 random restarts. The generative model G is set to be a pre-trained variational autoencoder (VAE) model with latent dimension k = 20. We use the VAE model trained by the authors of (Bora et al., 2017) directly, for which the encoder and decoder are both fully connected neural networks with two hidden layers, with the architecture being 20− 500− 500− 784. The VAE is trained by the Adam optimizer with a minibatch size of 100 and a learning rate of 0.001. The projection step PG(·) is solved by the Adam optimizer with a learning rate of 0.03 and 200 steps. In each iteration of TPower, the calculated entries are truncated to zero except for the largest q entries, where q ∈ N is a tuning parameter. Since for TPower, q is usually selected as an integer larger than the true sparsity level, and since it is unlikely that the image of the MNIST dataset can be well approximated by a k-sparse vector with k = 20, we choose a relatively large q, namely q = 150. Similarly to (Bora et al., 2017) and other related works, we only report the results on a test set that is unseen by the pre-trained VAE model, i.e., the training of G and the PPower computations do not use common data.8
1. Spiked covariance model (Example 1): The numerical results are shown in Figures 1 and 2. We observe from Figure 1 that Power and TPower attain poor reconstructions, and the generative prior based method PPower attains significantly better reconstructions. To illustrate the effect of the sample size m, we fix the SNR parameter β = 1 and vary m in {100, 200, 300, 400, 500}. In addition, to illustrate the effect of the SNR parameter β, we fix m = 300, and vary β in {0.6, 0.7, 0.8, 0.9, 1, 2, 3, 4}. From Figure 2, we observe that for these settings of m and β, PPower always leads to a much higher cosine similarity compared to Power and TPower, which is natural given the more precise modeling assumptions used.
8All experiments are run using Python 3.6 and Tensorflow 1.5.0, with a NVIDIA GeForce GTX 1080 Ti 11GB GPU. The corresponding code is available at https://github.com/liuzq09/ GenerativePCA.
2. Phase retrieval (Example 2): The results are shown in Figure 2 (Right) and Figure 3. Again, we can observe that PPower significantly outperforms Power and TPower. In particular, for sparse phase retrieval, when performing experiments on image datasets, even for the noiseless setting, solving an eigenvalue problem similar to (5) can typically only serve as a spectral initialization step, with a subsequent iterative algorithm being required to refine the initial guess. In view of this, it is notable that for phase retrieval with generative priors, PPower can return meaningful reconstructed images for m = 200, which is small compared to n = 784.

heading: 6 CONCLUSION
text: We have proposed a quadratic estimator for eigenvalue problems with generative models, and we showed that this estimator attains a statistical rate of order √
k logL m . We provided a projected power
method to efficiently solve (modulo the complexity of the projection step) the corresponding optimization problem, and showed that our method converges exponentially fast to a point achieving a statistical rate of order √
k logL m under suitable conditions.

heading: APPENDIX (GENERATIVE PRINCIPAL COMPONENT ANALYSIS,
text: LIU/LIU/GHOSH/HAN/SCARLETT, ICLR 2022)

heading: A PROOF OF LEMMA 1 (NON-DECREASING PROPERTY OF Q)
text: Since w(t+1) = PG ( Vw(t) ) and w(t) ∈ Range(G), we have
‖Vw(t) −w(t+1)‖2 ≤ ‖Vw(t) −w(t)‖2, (15)
and since ‖w(t+1)‖2 = ‖w(t)‖2 = 1, expanding the square gives〈 Vw(t),w(t+1) 〉 ≥ Q(w(t)). (16)
Then, we obtain Q(w(t+1)) = 〈 Vw(t+1),w(t+1) 〉 (17)
= Q(w(t+1) −w(t)) + 2 〈 V(w(t+1) −w(t)),w(t) 〉 +Q(w(t)) (18)
≥ 2 〈 V(w(t+1) −w(t)),w(t) 〉 +Q(w(t)) (19)
≥ Q(w(t)), (20)
where (18) follows by writing w(t+1) = w(t) + (w(t+1) −w(t)) and expanding, (19) follows from the assumption that V is PSD, and (20) follows from (16).

heading: B PROOFS FOR SPIKED MATRIX AND PHASE RETRIEVAL EXAMPLES
text: Before proceeding, we present the following standard definitions. Definition 1. A random variable X is said to be sub-Gaussian if there exists a positive constant C such that (E [|X|p])1/p ≤ C√p for all p ≥ 1. The sub-Gaussian norm of a sub-Gaussian random variable X is defined as ‖X‖ψ2 := supp≥1 p−1/2 (E [|X|p])
1/p. Definition 2. A random variableX is said to be sub-exponential if there exists a positive constantC such that (E [|X|p]) 1 p ≤ Cp for all p ≥ 1. The sub-exponential norm of X is defined as ‖X‖ψ1 := supp≥1 p −1 (E [|X|p]) 1 p .
The following lemma states that the product of two sub-Gaussian random variables is subexponential, regardless of the dependence between them. Lemma 3. (Vershynin, 2018, Lemma 2.7.7) Let X and Y be sub-Gaussian random variables (not necessarily independent). Then XY is sub-exponential, and satisfies
‖XY ‖ψ1 ≤ ‖X‖ψ2‖Y ‖ψ2 . (21)
The following lemma provides a useful concentration inequality for the sum of independent subexponential random variables. Lemma 4. (Vershynin, 2010, Proposition 5.16) Let X1, . . . , XN be independent zero-mean subexponential random variables, and K = maxi ‖Xi‖ψ1 . Then for every α = [α1, . . . , αN ]T ∈ RN and ≥ 0, it holds that
P (∣∣∣ N∑
i=1
αiXi ∣∣∣ ≥ ) ≤ 2 exp(−c ·min( 2 K2‖α‖22 , K‖α‖∞ )) , (22)
where c > 0 is an absolute constant. In particular, with α = [
1 N , . . . , 1 N
]T , we have
P (∣∣∣ 1 N N∑ i=1 Xi ∣∣∣ ≥ ) ≤ 2 exp(−c ·min(N 2 K2 , N K )) . (23)
The widely-used notion of an -net is introduced as follows. Definition 3. Let (X , d) be a metric space, and fix > 0. A subset S ⊆ X is said be an -net of X if, for all x ∈ X , there exists some s ∈ S such that d(s, x) ≤ . The minimal cardinality of an -net of X , if finite, is denoted C(X , ) and is called the covering number of X (at scale ).
The following lemma provides a useful upper bound for the covering number of the unit sphere. Lemma 5. (Vershynin, 2010, Lemma 5.2) The unit Euclidean sphere SN−1 equipped with the Euclidean metric satisfies for every > 0 that
C(SN−1, ) ≤ ( 1 + 2 )N . (24)
The following lemma provides an upper bound for the spectral norm of a symmetric matrix. Lemma 6. (Vershynin, 2010, Lemma 5.4) Let X be a symmetric N × N matrix, and let C be an -net of SN−1 for some ∈ [0, 1/2). Then,
‖X‖2→2 = sup r∈SN−1 |〈Xr, r〉| ≤ (1− 2 )−1 sup r∈C |〈Xr, r〉|. (25)
With the above auxiliary results in place, we provide the proofs of Assumption 2 holding for the two examples described in Section 3.
B.1 SPIKED COVARIANCE MODEL (EXAMPLE 1)
As per Assumption 2, fix two finite signal sets S1 and S2. For r = 1, we have xi = √ βuis + zi and a direct calculation gives E[V] = V̄ = βssT . Recall also that ‖s‖2 = 1, ui ∼ N (0, 1), and zi ∼ N (0, In). It follows that for any s1 ∈ S1, we have that sT1 xi = √ βuis T s1 + z T i s1 is sub-Gaussian,
with the sub-Gaussian norm being upper bounded by C( √ β + 1)‖s1‖2. Similarly, we have for any
s2 ∈ S2 that ‖sT2 xi‖ψ2 ≤ C( √ β + 1)‖s2‖2. Applying Lemma 3, we deduce that (sT1 xi)(sT2 xi) is
sub-exponential, with the sub-exponential norm being upper bounded by C2( √ β + 1)2‖s1‖2‖s2‖2. In addition, from (9) and V̄ = βssT , we have
sT1 Es2 = s T 1 (V − V̄)s2 (26)
= 1
m m∑ i=1 ( (xTi s1)(x T i s2)− ( (sT1 s2) + β(s T s1)(s T s2) )) , (27)
and we observe that E[(xTi s1)(xTi s2)] = (sT1 s2)+β(sT s1)(sT s2). Then, from Lemma 4, we obtain that for any t > 0 satisfying m = Ω(t), the following holds with probability 1− e−Ω(t) (recall that C may vary from line to line):∣∣∣∣∣ 1m m∑ i=1 ( (xTi s1)(x T i s2)− ( (sT1 s2) + β(s T s1)(s T s2) ))∣∣∣∣∣ ≤ C(√β + 1)2‖s1‖2‖s2‖2 · √ t√ m , (28)
where we note that the assumption m = Ω(t) ensures that the first term is dominant in the minimum in (23). Taking a union bound over all s1 ∈ S1 and s2 ∈ S2, and setting t = log(|S1| · |S2|), we obtain with probability 1− e−Ω(log(|S1|·|S2|)) that (7) holds (with β being a fixed positive constant). Next, we bound |rTEr| for fixed r ∈ Sn−1, but this time consider t > 0 (different from the above t) satisfying t = Ω(m). In this case, we can follow the above analysis (with s1 and s2 both replaced by r), but the assumption t = Ω(m) means that when applying Lemma 4, the second term in the minimum in (23) is now the dominant one. As a result, for any t > 0 satisfying t = Ω(m), and any r ∈ Sn−1, we have with probability 1− e−Ω(t) that
|rTEr| = ∣∣∣∣∣ 1m m∑ i=1 ( (xTi r) 2 − ( 1 + β(sT r)2 ))∣∣∣∣∣ (29) ≤ C( √ β + 1) · t
m . (30)
From Lemma 5, there exists an (1/4)-net C 1 4
of Sn−1 satisfying log ∣∣C 1
4 ∣∣ ≤ n log 9. Taking a union bound over all r ∈ C 1
4 , and setting t = Cn, we obtain with probability 1− e−Ω(n) that
sup r∈C 1
4
∣∣rTEr∣∣ = O ( n m ) . (31)
Then, from Lemma 6, we have
‖E‖2→2 ≤ 2 sup r∈C 1
4
∣∣rTEr∣∣ = O ( n m ) . (32)
B.2 PHASE RETRIEVAL (EXAMPLE 2) Let W = 1m ∑m i=1 yiaia T i 1{l<yi<u} and W̄ = βss
T + γIn. It is shown in (Liu et al., 2021b, Lemma 8) that E [W] = W̄, (33) which implies E[V] = βssT = V̄. (34) Then, for any s1 ∈ S1 and s2 ∈ S2, we have
sT1 Es2 = s T 1 (V − V̄)s2 = sT1 (W − W̄)s2 (35)
= 1
m m∑ i=1 ( yi(a T i s1)(a T i s2)1{l<yi<u} − ( β(sT s1)(s T s2) + γ(s T 1 s2) )) . (36)
Since each ai has i.i.d. N (0, 1) entries, we observe that yi(aTi s1)(aTi s2)1{l<yi<u} is subexponential with the sub-exponential norm being upper bounded by Cu‖s1‖2‖s2‖2. In addition, from (33), we have E[yi(aTi s1)(aTi s2)1{l<yi<u}] = β(sT s1)(sT s2) + γ(sT1 s2). Then, from Lemma 4, we obtain that for any t > 0 satisfying m = Ω(t), with probability 1− e−Ω(t),∣∣∣∣∣ 1m m∑ i=1 ( yi(a T i s1)(a T i s2)1{l<yi<u} − (β(s T s1)(s T s2) + γ(s T 1 s2)) )∣∣∣∣∣ ≤ Cu‖s1‖2‖s2‖2 · √ t√ m . (37) Taking a union bound over all s1 ∈ S1 and s2 ∈ S2, and setting t = log(|S1| · |S2|), we obtain that with probability 1− e−Ω(log(|S1|·|S2|)), (7) holds as desired (with u being a fixed positive constant). In addition, similarly to (32), we have with probability 1− e−Ω(n) that ‖E‖2→2 = O ( n m ) .

heading: C EQUIVALENCE OF DISTANCES
text: The following lemma gives a useful equivalence between two distances.
Lemma 7. For any pair of unit vectors w1,w2 with wT1 w2 ≥ 0, we have
‖w1 −w2‖22 ≤ ‖w1wT1 −w2wT2 ‖2F ≤ 2‖w1 −w2‖22. (38)
Moreover, if wT1 w2 < 0, then the same holds with ‖w1 −w2‖2 replaced by ‖w1 + w2‖2.
Proof. When wT1 w2 ≥ 0, we have
‖w1wT1 −w2wT2 ‖2F = tr((w1wT1 −w2wT2 )T (w1wT1 −w2wT2 )) (39) = 2 ( 1− (wT1 w2)2 ) (40)
≥ 2 ( 1−wT1 w2 ) (41)
= ‖w1 −w2‖22, (42)
where (40) follows by expanding the product and writing tr(w1wT1 w1w T 1 ) = tr(w T 1 w1w T 1 w1) = (wT1 w1) 2 = 1 and handling the other terms similarly, and (41) follows since wT1 w2 ∈ (0, 1). In
addition, we have
‖w1wT1 −w2wT2 ‖2F = 2 ( 1− (wT1 w2)2 ) (43)
= 2 ( 1−wT1 w2 ) ( 1 + wT1 w2 ) (44)
≤ 4 ( 1−wT1 w2 ) (45)
= 2‖w1 −w2‖22, (46)
which gives the desired inequality. The case wT1 w2 < 0 is handled similarly

heading: D PROOF OF THEOREM 1 (GUARANTEE ON THE GLOBAL OPTIMUM)
text: Let the singular value decomposition (SVD) of V̄ be
V̄ = ŪD̄ŪT , (47)
where D̄ = Diag([λ̄1, . . . , λ̄n]), and Ū ∈ Rn is an orthonormal matrix with the first column being x̄. For i > 1, let the i-th column of Ū be ūi. Then, we have
v̂T V̄v̂ = λ̄1 ( x̄T v̂ )2 + ∑ i>1 λ̄i ( ūTi v̂ )2 (48)
≤ λ̄1 ( x̄T v̂ )2 + λ̄2 ∑ i>1 ( ūTi v̂ )2 (49)
= λ̄1 ( x̄T v̂ )2 + λ̄2 ( 1− ( x̄T v̂ )2) , (50)
where we use ( x̄T v̂ )2 + ∑ i>1 ( ūTi v̂ )2 = 1 in (50).
In addition, for any A ∈ Rn×n and any s1, s2 ∈ Rn, we have
sT1 As1 − sT2 As2 = (
s1 + s2 2 + s1 − s2 2
)T A ( s1 + s2
2 + s1 − s2 2 ) − (
s1 + s2 2 − s1 − s2 2
)T A ( s1 + s2
2 − s1 − s2 2
) (51)
= 2
( s1 + s2
2
)T A ( s1 − s2
2
) + 2 ( s1 − s2
2
)T A ( s1 + s2
2
) . (52)
In particular, when A is symmetric, we obtain
sT1 As1 − sT2 As2 = (s1 + s2)TA(s1 − s2). (53)
Let M be a (δ/L)-net of Bk2 (r); from (Vershynin, 2010, Lemma 5.2), we know that there exists such a net with
log |M | ≤ k log 4Lr δ . (54)
Since G is L-Lipschitz continuous, we have that G(M) is a δ-net of Range(G) = G(Bk2 (r)). We write v̂ = (v̂ − x̃) + x̃, (55) where x̃ ∈ G(M) satisfies ‖v̂− x̃‖2 ≤ δ. Suppose that x̄T v̂ ≥ 0; if this is not the case, we can use analogous steps to obtain an upper bound for ‖x̄ + v̂‖2 instead of ‖x̄− v̂‖2. We have
λ̄1 − λ̄2 2 · ‖x̄− v̂‖22 (56)
= (λ̄1 − λ̄2) ( 1− x̄T v̂ ) (57)
≤ (λ̄1 − λ̄2) ( 1− ( x̄T v̂ )2) (58)
= λ̄1 − ( λ̄1 ( x̄T v̂ )2 + λ̄2 ( 1− ( x̄T v̂ )2)) (59)
= x̄T V̄x̄− ( λ̄1 ( x̄T v̂ )2 + λ̄2 ( 1− ( x̄T v̂ )2)) (60)
≤ x̄T V̄x̄− v̂T V̄v̂ (61) = xTGV̄xG + (x̄ + xG)
T V̄(x̄− xG)− v̂T V̄v̂ (62) ≤ xTGV̄xG + 2λ̄1‖x̄− xG‖2 − v̂T V̄v̂ (63) = xTG(V −E)xG + 2λ̄1‖x̄− xG‖2 − v̂T (V −E)v̂ (64) ≤ v̂TEv̂ − xTGExG + 2λ̄1‖x̄− xG‖2 (65)
= x̃TEx̃ + 2
( v̂ − x̃
2
)T E ( v̂ + x̃
2
) + 2 ( v̂ + x̃
2
)T E ( v̂ − x̃
2
) − xTGExG + 2λ̄1‖x̄− xG‖2
(66)
≤ x̃TEx̃ + 2δ‖E‖2→2 − xTGExG + 2λ̄1‖x̄− xG‖2 (67)
= 2
( x̃ + xG
2
)T E ( x̃− xG
2
) + 2 ( x̃− xG
2
)T E ( x̃ + xG
2
) + 2δ‖E‖2→2 + 2λ̄1‖x̄− xG‖2
(68)
≤ 2C
√ k log 4Lrδ
m · ‖x̃− xG‖2 + 2δ‖E‖2→2 + 2λ̄1‖x̄− xG‖2 (69)
≤ 2C
√ k log 4Lrδ
m · (‖x̃− v̂‖2 + ‖v̂ − x̄‖2 + ‖x̄− xG‖2) + 2δ‖E‖2→2 + 2λ̄1‖x̄− xG‖2 (70)
≤ 2C
√ k log 4Lrδ
m · ‖v̂ − x̄‖2 +O
( δn
m
) +O ( (λ̄1 + n)‖x̄− xG‖2 ) , (71)
where:
• (57)–(58) follow from ‖x̄‖2 = ‖v̂‖2 = 1 and hence |x̄T v̂| ≤ 1;
• (60) follows since (λ̄1, x̄) are an eigenvalue-eigenvector pair for V̄ with ‖x̄‖2 = 1;
• (61) follows from (50);
• (62) follows from (53) with V̄ being symmetric and setting s1 = x̄, s2 = xG;
• (63) follows from ‖x̄ + xG‖2 ≤ 2 and ‖V̄‖2→2 = λ̄1;
• (64) follows since V̄ = V −E;
• (65) follows since v̂ is a globally optimal solution to (5) and xG ∈ Range(G);
• (66) follows from (52) with s1 = v̂ and s2 = x̃;
• (67) follows from (55) along with ‖v̂ − x̄‖2 ≤ δ and ‖v̂ + x̃‖2 ≤ 2;
• (68) follows from (52);
• (69) follows from Assumption 2 (with S1 = S2 being G(M) shifted by xG) and (54);
• (70) follows from the triangle inequality;
• (71) follows by substituting ‖v̂ − x̃‖2 ≤ δ, along with the assumptions ‖E‖2→2 = O(n/m), m = Ω ( k log Lrδ ) , and n = O (√k log Lrδ m ) .
From (71), we have the following when v̂T x̄ ≥ 0:
‖v̂ − x̄‖2 = O
(√ k log Lrδ m ) λ̄1 − λ̄2 +O (√ δn/m λ̄1 − λ̄2 ) +O √ (λ̄1 + n)‖x̄− xG‖2 λ̄1 − λ̄2  . (72)
As mentioned earlier, if v̂T x̄ < 0, we have the same upper bound as in (72) for ‖v̂+x̄‖2. Therefore, we obtain
‖v̂v̂T − x̄x̄T ‖F = √ 2 ( 1− (x̄T v̂)2 )
(73) = √ 2 (1− x̄T v̂) (1 + x̄T v̂) (74) ≤ √ 2 min{‖v̂ − x̄‖2, ‖v̂ + x̄‖2} (75)
=
O (√ k log Lrδ m ) λ̄1 − λ̄2 +O (√ δn/m λ̄1 − λ̄2 ) +O √ (λ̄1 + n)‖x̄− xG‖2 λ̄1 − λ̄2  , (76) where (73) follows from (40), (75) follows from ‖v̂ ± x̄‖22 = 2(1 ± x̄T v̂), and (76) follows from (72).

heading: E PROOF OF LEMMA 2 (AUXILIARY RESULT FOR PPOWER ANALYSIS)
text: By the assumption Range(G) ⊆ Sn−1, for any x ∈ Rn and a > 0, we have
PG(x) = PG(ax), (77)
which is seen by noting that when comparing ‖x−a‖2 with ‖x−b‖2 (in accordance with projection mapping to the closest point), as long as ‖a‖2 = ‖b‖2, the comparison reduces to comparing 〈x,a〉 with 〈x,b〉, so is invariant to positive scaling of x. Let η̄ = 1/λ̄1 > 0 and ŝ = PG(Vs). Then, we have ŝ = PG(Vs) = PG(η̄Vs). Since x̄ ∈ Range(G), we have
‖η̄Vs− ŝ‖2 ≤ ‖η̄Vs− x̄‖2. (78)
This is equivalent to ‖(η̄Vs− x̄) + (x̄− ŝ)‖22 ≤ ‖η̄Vs− x̄‖22, (79)
and expanding the square gives
‖x̄− ŝ‖22 ≤ 2〈η̄Vs− x̄, ŝ− x̄〉. (80)
Note also that from V̄x̄ = λ̄1x̄, we obtain x̄ = η̄V̄x̄, which we will use throughout the proof.
For δ > 0, let M be a (δ/L)-net of Bk2 (r); from Lemma 5, there exists such a net with
log |M | ≤ k log 4Lr δ . (81)
By the L-Lipschitz continuity of G, we have that G(M) is a δ-net of Range(G) = G(Bk2 (r)). We write
s = (s− s0) + s0, ŝ = (ŝ− s̃) + s̃, (82)
where s̃ ∈ G(M) satisfies ‖ŝ− s̃‖2 ≤ δ, and s0 ∈ G(M) satisfies ‖s− s0‖2 ≤ δ. Then, we have
〈η̄Vs− x̄, ŝ− x̄〉 = 〈η̄V̄(s− x̄), ŝ− x̄〉+ 〈η̄Es, ŝ− x̄〉, (83)
which follows from V = V̄ + E and x̄ = η̄V̄x̄. In the following, we control the two terms in (83) separately.
1. The term 〈η̄V̄(s−x̄), ŝ−x̄〉: We decompose s = αx̄+βt and ŝ = α̂x̄+ β̂t̂, where ‖t‖2 = ‖t̂‖2 = 1 and tT x̄ = t̂T x̄ = 0. Since ‖s‖2 = ‖ŝ‖2 = 1, we have α2 + β2 = α̂2 + β̂2 = 1. In addition, we have α = sT x̄ and α̂ = ŝT x̄. Recall that in (47), we write the SVD of V̄ as V̄ = ŪD̄ŪT . Since tT x̄ = 0, we can write t as t = ∑ i>1 hiūi. In addition, since
‖t‖2 = 1, we have ∑ i>1 h 2 i = 1. Hence, by the Cauchy-Schwarz inequality, we have
|〈V̄t, t̂〉| ≤ ‖V̄t‖2 (84)
= ∥∥∥∥∥∑ i>1 λ̄ihiūi ∥∥∥∥∥ 2
(85)
= √∑ i>1 λ̄2ih 2 i (86)
≤ √ λ̄22 ∑ i>1 h2i (87)
= λ̄2. (88)
Therefore, we obtain
|〈η̄V̄(s− x̄), ŝ− x̄〉| = |〈(α− 1)x̄ + η̄βV̄t, (α̂− 1)x̄ + β̂t̂〉| (89) = |(α− 1)(α̂− 1) + η̄ββ̂〈V̄t, t̂〉| (90) ≤ (1− α)(1− α̂) + η̄|ββ̂|λ̄2 (91)
= (1− α)(1− α̂) + γ̄ √ 1− α2 √ 1− α̂2, (92)
where (89) uses ηV̄x̄ = x̄, and (90) uses ‖x̄‖2 = 1 and 〈x̄, t〉 = 0. 2. The term 〈η̄Es, ŝ− x̄〉: We have
|〈η̄Es, ŝ− x̄〉| = 〈η̄E((s− s0) + s0), ŝ− x̄〉 (93) = 〈η̄E(s− s0), ŝ− x̄〉+ 〈η̄Es0, (ŝ− s̃) + (s̃− x̄)〉 (94)
≤ η̄‖E‖2→2δ‖ŝ− x̄‖2 + η̄‖E‖2→2δ +O √k log Lrδ m  · ‖s̃− x̄‖2 (95)
≤ O (δ‖E‖2→2) +O √k log Lrδ m  · ‖ŝ− x̄‖2, (96) where (95) follows from Assumption 2 and (81), and (96) follows from η̄ = 1/λ̄1, along with the fact that we assumed λ̄1 = Θ(1).
Note that ‖x̄− ŝ‖22 = 2(1− ŝT x̄) = 2(1− α̂). Hence, and using (80), (83), (92), and (96), we obtain 2(1− α̂) ≤ 2 ( (1− α)(1− α̂) + γ̄ √ 1− α2 √ 1− α̂2 )
+O (δ‖E‖2→2) +O √k log Lrδ m  ·√2(1− α̂). (97) Using 2(1− α̂)− 2(1− α)(1− α̂) = 2α(1− α̂), √ 1− α2 = √ 1− α √ 1 + α ≤ √ 2(1− α), and similarly √ 1− α̂2 ≤ √ 2(1− α̂), we obtain from (97) that
2α(1− α̂) ≤ 2γ̄ √ 2(1− α) √ 2(1− α̂) +O √k log Lrδ m  ·√2(1− α̂) +O (δ‖E‖2→2) . (98) Since ‖ŝ− x̄‖22 = 2(1− α̂) and ‖s− x̄‖22 = 2(1− α), this is equivalent to
α‖ŝ− x̄‖22 ≤ 2γ̄‖s− x̄‖2 +O √k log Lrδ
m
 ‖ŝ− x̄‖2 +O (δ‖E‖2→2) . (99)
This equation is of the form az2 ≤ bz + c (where z = ‖ŝ− x̄‖2), and using a simple application of the quadratic formula, we obtain
‖ŝ− x̄‖2 ≤ 2γ̄‖s− x̄‖2
sT x̄ +O  1 sT x̄ √k log Lrδ m + √ sT x̄ · δ‖E‖2→2  (100) ≤ 2γ̄‖s− x̄‖2
sT x̄ +O
( 1
sT x̄
√ k log(nLr)
m
) , (101)
where we use α = sT x̄ in (100), and we use the assumption ‖E‖2→2 = O(n/m) and set δ = 1/n in (101).

heading: F PROOF OF THEOREM 2 (MAIN THEOREM FOR PPOWER)
text: Suppose for the time being that (13) holds for at least one index t ≥ t0 (we will later verify that this is the case), and let T0 ≥ t0 be the smallest such index. Thus, we have
‖w(T0) − x̄‖2 ≤ C
(1− µ0)ν
√ k log(nLr)
m . (102)
Note that according to the theorem statement, 1−µ0 is bounded away from zero. Using ‖w(T0)‖2 = ‖x̄‖2 = 1 and the assumption that m ≥ Cν,τ · k log(nLr) with Cν,τ > 0 being large enough, we deduce from (102) that ‖w(T0) − x̄‖2 is sufficiently small such that
x̄Tw(T0) ≥ 1− τ. (103)
Next, using the assumption 2γ̄ + ν ≤ 1− τ , we write
2γ̄
(1− µ0)ν(1− τ) +
1
1− τ = 2γ̄ + (1− µ0)ν (1− µ0)ν(1− τ) ≤ 2γ̄ + ν (1− µ0)ν(1− τ) ≤ 1 (1− µ0)ν . (104)
Then, from Lemma 2, we obtain
‖w(T0+1) − x̄‖2 ≤ 2γ̄
x̄Tw(T0) · ‖w(T0) − x̄‖2 +
C
x̄Tw(T0)
√ k log(nLr)
m (105)
≤ 2γ̄ 1− τ · C (1− µ0)ν
√ k log(nLr)
m +
C
1− τ
√ k log(nLr)
m (106)
≤ C (1− µ0)ν
√ k log(nLr)
m , (107)
where (106) follows from (102)–(103), and (107) follows from (104). Thus, we have transferred (102) from T0 to T0 + 1, and proceeding by induction, we obtain
‖w(t) − x̄‖2 ≤ C
(1− µ0)ν
√ k log(nLr)
m (108)
for al t ≥ T0. Next, we consider t ∈ [t0, T0). Again using Lemma 2 (with ŝ = w(t0+1) = PG(Vw(t0))), we have
‖w(t0+1) − x̄‖2 ≤ µ0‖w(t0) − x̄‖2 + C 2γ̄ + ν · √ k log(nLr) m , (109)
where we recall that µ0 = 2γ̄x̄Tw(t0) = 2γ̄ 2γ̄+ν < 1, and note that the denominator in the second term of (109) follows since x̄Tw(t0) = 2γ̄ + ν. Supposing that t0 < T0 (otherwise, the above analysis for t ≥ T0 alone is sufficient), we have that (13) is reversed at t = t0:
‖w(t0) − x̄‖2 > C
(1− µ0)ν
√ k log(nLr)
m . (110)
This means that we can upper bound the second term in (109) by Cν · √ k log(nLr) m < (1−µ0)‖w
(t0)− x̄‖2, which gives ‖w(t0+1) − x̄‖2 < ‖w(t0) − x̄‖2. (111) Squaring both sides, expanding, and canceling the norms (which all equal one), we obtain
x̄Tw(t0+1) > x̄Tw(t0), (112)
and by induction, we obtain that {x̄Tw(t)}t∈[t0,T0) is monotonically increasing.
Recall that we assume that λ̄1 = Θ(1), and that T0 = t0 + ∆0 is the smallest integer such that (13) holds. To verify that T0 is finite and upper bound ∆0, we consider the following three cases:
• µ0 = 0 (or equivalently, γ̄ = λ̄2 = 0): In this case, (109) gives T0 = t0 +1 (or T0 = t0, which we already addressed above). Thus, we have ∆0 ≤ 1, as stated in the theorem. • µ0 = o(1) (or equivalently, γ̄ = o(1) and λ̄2 = o(1)): Since {x̄Tw(t)}t∈[t0,T0) is monotonically increasing, for any positive integer ∆ with t0 +∆ ≤ T0, by applying Lemma 2 (or (109)) multiple times, we obtain9
‖w(t0+∆) − x̄‖2 ≤ µ∆0 ‖w(t0) − x̄‖2 + 1− µ∆0 1− µ0 · C 2γ̄ + ν
· √ k log(nLr)
m (113)
≤ µ∆0 ‖w(t0) − x̄‖2 + 1 1− µ0 · C 2γ̄ + ν · √ k log(nLr) m . (114)
Then, if we choose ∆0 ∈ N such that
µ∆0−10 ≤ C 2ν · √ k log(nLr) m , (115)
we obtain from (114) that
‖w(t0+∆0) − x̄‖2 ≤ µ∆00 ‖w(t0) − x̄‖2 + 1 1− µ0 · C 2γ̄ + ν · √ k log(nLr) m (116)
≤ 2µ∆00 + 1 1− µ0 · C 2γ̄ + ν · √ k log(nLr) m (117)
< Cµ0 ν(1− µ0) · √ k log(nLr) m + 1 1− µ0 · C 2γ̄ + ν · √ k log(nLr) m (118)
= C (1− µ0)ν · √ k log(nLr) m , (119)
where (117) follows from ‖w(t0) − x̄‖2 ≤ 2, (118) follows from (115) and 1 < 11−µ0 , and (119) follows from µ0 = 2γ̄2γ̄+ν , which implies µ0 ν + 1 2γ̄+ν = µ0(2γ̄+ν)+ν ν(2γ̄+ν) = 2γ̄+ν ν(2γ̄+ν) = 1 ν . Observe that (119) coincides with (13), and since µ0 = o(1), we obtain from (115) that ∆0 = O ( log ( m k log(nLr) )) as desired.
• µ0 = Θ(1) (or equivalently, γ̄ = Θ(1) and λ̄2 = Θ(1)): Recall that we only need to focus on the case T0 > t0. This means that (110) holds, implying that we can upper bound the second term in (109) by (1−µ0)ν2γ̄+ν · ‖w (t0) − x̄‖2, yielding
‖w(t0+1) − x̄‖2 < µ0‖w(t0) − x̄‖2 + (1− µ0)ν
2γ̄ + ν · ‖w(t0) − x̄‖2 (120)
= 2γ̄ + (1− µ0)ν
2γ̄ + ν · ‖w(t0) − x̄‖2 (121)
= (1− ξ)‖w(t0) − x̄‖2, (122) 9In simpler notation, if zt+1 ≤ azt+b, then we get zt+2 ≤ a2zt+(1+a)b, then zt+3 ≤ a3zt+(1+a+a2)b,
and so on, and then we can apply 1 + a+ . . .+ ai−1 = 1−a i
1−a for a 6= 1.
where ξ = µ0ν2γ̄+ν = µ0(1− µ0) = Θ(1). With the distance to x̄ shrinking by a constant factor in each iteration according to (122), and the initial distance ‖w(t0) − x̄‖2 being at most 2 due to the vectors having unit norm, we deduce that ∆0 = O ( log ( m k log(nLr) )) iterations suffice to ensure that (13) holds for t = t0 + ∆0.

heading: G COMPARISON OF ANALYSIS TO DESHPANDE ET AL. (2014)
text: As mentioned in Section 4, our analysis is significantly different from that of Deshpande et al. (2014) despite using a similar assumption on the initialization. We highlight the differences as follows:
1. Perhaps the most significant difference is that the proof of (Deshpande et al., 2014, Theorem 3) is highly dependent on the Moreau decomposition, which is only valid for a closed convex cone (see (Deshpande et al., 2014, Definition 1.2)). In particular, the Moreau decomposition needs to be used at the beginning of the proof of (Deshpande et al., 2014, Theorem 3), such as Eqs. (18) and (19) in the supplementary material therein. We do not see a way for the proof to proceed without the Moreau decomposition, and our Range(G) may be very different from a convex cone.
2. We highlight that one key observation in our proof of Lemma 2 (and thus Theorem 2) is that for a generative model G with Range(G) ⊆ Sn−1, and any x ∈ Rn and a > 0, we have PG(ax) = PG(x) (Eq. (77)). This enables us to derive the important equation ŝ = PG(Vs) = PG(η̄Vs). We are not aware of a similar idea being used in the proof of (Deshpande et al., 2014, Theorem 3).
3. In the PPower method in (Deshpande et al., 2014), the authors need to add ρIn with ρ > 0 to the observed data matrix V to improve the convergence. In particular, they mention in the paragraph before the statement of Theorem 3 that “the memory term ρvt is necessary for our proof technique to go through”. In contrast, our proof of Theorem 2 does not require adding such terms, even when our data model is restricted to the spiked Wigner model considered in (Deshpande et al., 2014).
4. We consider a matrix model that is significantly more general than the spiked Wigner model studied in (Deshpande et al., 2014).

heading: H NUMERICAL RESULTS FOR THE FASHION-MNIST DATASET
text: The Fashion-MNIST dataset consists of Zalando’s article images with a training set of 60, 000 examples and a test set of 10, 000 examples. The size of each image in the Fashion-MNIST dataset is also 28× 28, and thus n = 784. The generative model G is set to be a boundary-seeking generative adversarial network (BEGAN). The BEGAN architecture is summarized as follows.10 The generator has latent dimension k = 62 and four layers. The first two are fully connected layers with the architecture 62 − 1024 − 6272, and with ReLU activation functions. The output of the second layer, reshaped to 128 × 7 × 7, is forwarded to a deconvolution layer with kernel size 4 and stride 2. The third layer uses ReLU activations and has output size 64× 14× 14, where 64 is the number of channels. The fourth layer is a deconvolution layer with kernel size 4 and strides 2, and it uses ReLU activations and has output size 1× 28× 28, where the number of channels is 1. The BEGAN is trained with a mini-batch size of 256, a learning rate of 0.0002, and 100 epochs. The other parameters are the same as those for the MNIST dataset. We perform two sets of experiments, considering the spiked covariance and phase retrieval models separately. The corresponding results are reported in Figures 4, 5, 6, and 7. From these figures, again, we can observe clear superiority of PPower to Power and TPower. We note that for the Fashion-MNIST dataset, some of the images are not sparse in the natural basis, but we observe from Figures 4 and 6 that even for the sparsest images (sandals), PPower also significantly outperforms TPower.
10Further details of the architecture can be found at https://github.com/hwalsuklee/ tensorflow-generative-model-collections.

heading: I NUMERICAL RESULTS FOR THE CELEBA DATASET
text: The CelebA dataset consists of more than 200, 000 face images of celebrities, where each input image is cropped to a 64 × 64 RGB image with n = 64 × 64 × 3 = 12288. The generative model G is set to be a pre-trained Deep Convolutional Generative Adversarial Networks (DCGAN) model with latent dimension k = 100. We use the DCGAN model trained by the authors of (Bora et al., 2017) directly. We select the best estimate among 2 random restarts. The Adam optimizer with 100 steps and a learning rate of 0.1 is used for the projection operator.
For the images of the CelebA dataset, the corresponding vectors are clearly not sparse in the natural basis. To make a fairer comparison to the sparsity-based method TPower, we convert the original images to the wavelet basis, and perform TPower on these converted images. The obtained results of TPower are then converted back to the vectors in the natural basis. The corresponding method is denoted by TPowerW, with “W” referring to the conversion to images in the wavelet basis. In each iteration of TPower and TPowerW, the calculated entries are truncated to zero except for the largest q entries. For CelebA, q is set to be 2000. Other parameters are the same as those for the MNIST dataset. We also perform two sets of experiments, considering the spiked covariance and phase retrieval models separately. The corresponding results are reported in Figures 8 and 9. From these figures, we can observe the superiority of PPower to Power, TPower and TPowerW, whereas TPowerW only marginally improves over TPower.

Review summaries:
Principal component analysis (PCA) is one of the most commonly used techniques for dimension reduction in data and for explaining variability of the data in terms of its eigenvalues/vectors. Unfortunately in high-dimensions, PCA is mathematically inconsistent: this is well-known from random matrix literature. Hence the standard approach happens to be to assume that there is sparsity in the leading eigenvectors. In this paper however, motivated by Bora et. al (2017), the study is based on the conjecture/premise that sparsity assumption may be replaceable by the assumption that the leading eigenvector is close to the range of a deep generative model.   

This work studies PCA under the contraint that the principal eigenvectors lie in the range of some fixed generative model, e.g., a neural network with fixed weights.  More precisely, given an $L$-Lipschitz function $G : \mathbb R^k \mapsto R(G)$ the authors try to find a vector $w \in R(G)$ that maximizes $w^T V w$, where $V = \bar{V} + E$, is a (noisy) version of the "true" psd matrix $\bar{V}$.  The corresponding true solution $\bar{x}$ is defined to be the top eigenvector of $\bar{V}$.  They assume that the range $R(G)$ of $G$ is a subset of the unit ball of $\mathbb R^n$ and that the error matrix $E$ is not too large in the sense that for any sets $S_1, S_2$ with $m = \Omega( \log( |S_1| |S_2|))$ it holds that $|s_1^T E s_2| \leq \sqrt{ \log(|S_1| |S_2|)/ m} \|s_1\|_2 \|s_2\|_2$.  This error bound is satisfied (via standard concentration inequalities) by the spiked covariance and phase retrieval model assuming that the number of samples $m$ is large enough.  The authors show two main results.  First, they show that under the above assumptions the global maximizer given the noisy matrix $V$, i.e., the maximizer $\hat{v}$ of $w^T V w$ constrained on $R(G)$ has error roughly $\|\hat{v} - \bar{x}\|_2 \leq \sqrt{k \log L/m}$.  Next, they show that the power iteration algorithm together with a projection step onto the contraint set $R(G)$ will also eventually converge to a vector with similar error.  They also give experimental results showing that adding this projection step in the power iteration leads to improved recovery for the spiked covariance and phase retrieval problems for the MNIST and Fashion-MNIST datasets in the case where the generative model $G$ is a pre-trained variational autoencoder of small latent dimension $k$.

This paper studies the problem of PCA with a generative model setup. More precisely, they study the following setting  Given $\overline V = V + E$ where $V$ is PSD and $E$ is a perturbation matrix, find $\text{max}_w~ w^T V w$ where $w \in \text{Range}(G)$ where $G : B_2^k(r) \rightarrow \mathbb S^{n-1}$ is $L$-lipschitz.    The algorithm they analyze is a projected power method (PPower), similar to iterative hard thresholding for the case when $w$ is sparse. They demonstrate that if $w^{(i)} $ for $i \leq t$ are the $w$'s at each iteration, then for $Q(x) = x^T V x$, $Q(w^{(i)})$ is monotonically non-decreasing.   The subsequent results rely on the following two assumptions, which the authors demonstrate in the setting of the spiked wigner model, spiked covariance model and phase retrieval.    A. There is a gap between the largest and second largest eigenvalue of $\overline V$.  B. For $S_1, S_2 \subset \mathbb R^n$ satisfying $m = \Omega( \log(|S_1| \cdot |S_2|))$, for all $s_i \in S_i$, $|s_1^T E s_2| \leq C \sqrt{\frac{ \log(|S_1| \cdot |S_2|)}{m}} \cdot \|s_1\| \cdot \|s_2\|$.   The authors then provide guarantees for the optimizer in the context of these three settings, under these assumptions. More precisely, if $\overline v$ is the global optimizer, $\overline x$  is the true hidden signal, $x_G$ is the closest point in the range of the GAN to $\overline x$ and $\Delta = \overline \lambda_1 - \overline \lambda_2$ is the gap between the largest and second largest eigenvector of $\overline V$, they show,  $ \min \{\|\overline x - \overline v\|, \|\overline x + \overline v\| \} \leq O(\sqrt{k \log(Lr/\delta)/m}/\Delta) + 4(\overline \lambda_1/\Delta) \cdot \|x_G - \overline x\| + O(\delta n / m \Delta)$  Additionally, they demonstrate conditions on the initial $w_0$ which is input to PPower under which the algorithm has a linear rate of convergence, however they also need certain other assumptions which are satisfied in natural settings (such as the range of the GAN being positive).   Finally, they demonstrate experiments on MNIST for the spiked covariance model where they compare their algorithm (PPower) to the truncated power method (TPower) and the vanilla power method as a baseline. They show that the performance of PPower is much greater.