Year: 2021

Abstract:
We present a novel view on principal component analysis (PCA) as a competitive game in which each approximate eigenvector is controlled by a player whose goal is to maximize their own utility function. We analyze the properties of this PCA game and the behavior of its gradient based updates. The resulting algorithm—which combines elements from Oja’s rule with a generalized GramSchmidt orthogonalization—is naturally decentralized and hence parallelizable through message passing. We demonstrate the scalability of the algorithm with experiments on large image datasets and neural network activations. We discuss how this new view of PCA as a differentiable game can lead to further algorithmic developments and insights.

heading: 1 INTRODUCTION
text: The principal components of data are the vectors that align with the directions of maximum variance. These have two main purposes: a) as interpretable features and b) for data compression. Recent methods for principal component analysis (PCA) focus on the latter, explicitly stating objectives to find the k-dimensional subspace that captures maximum variance (e.g., (Tang, 2019)), and leaving the problem of rotating within this subspace to, for example, a more efficient downstream singular value (SVD) decomposition step1. This point is subtle, yet critical. For example, any pair of twodimensional, orthogonal vectors spans all of R2 and, therefore, captures maximum variance of any two-dimensional dataset. However, for these vectors to be principal components, they must, in addition, align with the directions of maximum variance which depends on the covariance of the data. By learning the optimal subspace, rather than the principal components themselves, objectives focused on subspace error ignore the first purpose of PCA. In contrast, modern nonlinear representation learning techniques focus on learning features that are both disentangled (uncorrelated) and low dimensional (Chen et al., 2016; Mathieu et al., 2018; Locatello et al., 2019; Sarhan et al., 2019).
It is well known that the PCA solution of the d-dimensional dataset X ∈ Rn×d is given by the eigenvectors of X>X or equivalently, the right singular vectors of X . Impractically, the cost of computing the full SVD scales with O(min{nd2, n2d})-time and O(nd)-space (Shamir, 2015; Tang, 2019). For moderately sized data, randomized methods can be used (Halko et al., 2011). Beyond this, stochastic—or online—methods based on Oja’s rule (Oja, 1982) or power iterations (Rutishauser, 1971) are common. Another option is to use streaming k-PCA algorithms such as Frequent Directions (FD) (Ghashami et al., 2016) or Oja’s algorithm2 (Allen-Zhu and Li, 2017) with storage complexity O(kd). Sampling or sketching methods also scale well, but again, focus on the top-k subspace (Sarlos, 2006; Cohen et al., 2017; Feldman et al., 2020).
In contrast to these approaches, we view each principal component (equivalently eigenvector) as a player in a game whose objective is to maximize their own local utility function in controlled competition with other vectors. The proposed utility gradients are interpretable as a combination of Oja’s rule and a generalized Gram-Schmidt process. We make the following contributions:
• A novel formulation of PCA as finding the Nash equilibrium of a suitable game, • A sequential, globally convergent algorithm for approximating the Nash on full-batch data,
1After learning the top-k subspace V ∈ Rd×k, the rotation can be recovered via an SVD of XV . 2FD approximates the top-k subspace; Oja’s algorithm approximates the top-k eigenvectors.
• A decentralized algorithm with experiments demonstrating the approach as competitive with modern streaming k-PCA algorithms on synthetic and real data, • In demonstration of the scaling of the approach, we compute the top-32 principal components
of the matrix of RESNET-200 activations on the IMAGENET dataset (n ≈ 106, d ≈ 20 ·106).
Each of these contributions is important. Novel formulations often lead to deeper understanding of problems, thereby, opening doors to improved techniques. In particular, k-player games are in general complex and hard to analyze. In contrast, PCA has been well-studied. By combining the two fields we hope to develop useful analytical tools. Our specific formulation is important because it obviates the need for any centralized orthonormalization step and lends itself naturally to decentralization. And lastly, theory and experiments support the viability of this approach for continued research.

heading: 2 PCA AS AN EIGEN-GAME
text: We adhere to the following notation. Vectors and matrices meant to approximate principal components (equivalently eigenvectors) are designated with hats, v̂ and V̂ respectively, whereas true principal components are v and V . Subscripts indicate which eigenvalue a vector is associated with. For example, vi is the ith largest eigenvector. In this work, we will assume each eigenvalue is distinct. By an abuse of notation, vj<i refers to the set of vectors {vj |j ∈ {1, . . . , i− 1}} and are also referred to as the parents of vi (vi is their child). Sums over indices should be clear from context, e.g.,∑ j<i = ∑i−1 j=1. The Euclidean inner product is written 〈u, v〉 = u>v. We denote the unit sphere by Sd−1 and simplex by ∆d−1 in d-dimensional ambient space.
Outline of derivation As argued in the introduction, the PCA problem is often mis-interpreted as learning a projection of the data into a subspace that captures maximum variance (equiv. maximizing the trace of a suitable matrix R introduced below). This is in contrast to the original goal of learning the principal components. We first develop the intuition for deriving our utility functions by (i) showing that only maximizing the trace of R is not sufficient for recovering all principal components (equiv. eigenvectors), and (ii) showing that minimizing off-diagonal terms in R is a complementary objective to maximizing the trace and can recover all components. We then consider learning only the top-k and construct utilities that are consistent with findings in (i) and (ii), equal the true eigenvalues at the Nash of the game we construct, and result in a game that is amenable to analysis.
Derivation of player utilities. The eigenvalue problem for a symmetric matrix X>X = M ∈ Rd×d is to find a matrix of d orthonormal column vectors V (implies V is full-rank) such that MV = V Λ with Λ diagonal. Given a solution to this problem, the columns of V are known as eigenvectors and corresponding entries in Λ are eigenvalues. By left-multiplying by V > and recalling V >V = V V > = I by orthonormality (i.e., V is unitary), we can rewrite the equality as
V >MV = V >V Λ unitary = Λ. (1)
Let V̂ denote a guess or estimate of the true eigenvectors V and define R(V̂ ) def= V̂ >MV̂ . The PCA problem is often posed as maximizing the trace of R (equiv. minimizing reconstruction error):
max V̂ >V̂=I {∑ i Rii = Tr(R) = Tr(V̂ >MV̂ ) = Tr(V̂ V̂ >M) = Tr(M) } . (2)
Surprisingly, the objective in (2) is independent of V̂ , so it cannot be used to recover all (i.e., k = d) the eigenvectors of M—(i). Alternatively, Equation (1) implies the eigenvalue problem can be phrased as ensuring all off-diagonal terms of R are zero, thereby ensuring R is diagonal—(ii):
min V̂ >V̂=I ∑ i6=j R2ij . (3)
It is worth further examining the entries of R in detail. Diagonal entries Rii = 〈v̂i,Mv̂i〉 are recognized as Rayleigh quotients because ||v̂i|| = 1 by the constraints. Off-diagonal entries Rij = 〈v̂i,Mv̂j〉 measure alignment between v̂i and v̂j under a generalized inner product 〈·, ·〉M .
So far, we have considered learning all the eigenvectors. If we repeat the logic for the top-k eigenvectors with k < d, then by Equation (1), R must still be diagonal. V is not square, so V V > 6= I , but assuming V is orthonormal as before, we have V V > = P is a projection matrix. Left-multiplying Equation (1) by V now reads (PM)V = V Λ so we are solving an eigenvalue problem for a subspace of M .
If we only desire the top-k eigenvectors, maximizing the trace encourages learning a subspace spanned by the top-k eigenvectors, but does not recover the eigenvectors themselves. On the other hand, Equation (3) places no preference on recovering large over small eigenvectors, but does enforce the columns of V̂ to actually be eigenvectors. The preceding exercise is intended to introduce minimizing the off-diagonal terms of R as a possible complementary objective for solving top-k PCA. Next, we will use these two objectives to construct utility functions for each eigenvector v̂i.
We want to combine the objectives to take advantage of both their strengths. A valid proposal is
max V̂ >V̂=I ∑ i Rii − ∑ i6=j R2ij . (4)
However, this objective ignores the natural hierarchy of the top-k eigenvectors. For example, v̂1 is penalized for aligning with v̂k and vice versa, but v̂1, being the estimate of the largest eigenvector, should be free to search for the direction that captures the most variance independent of the locations of the other vectors. Instead, first consider solving for the top-1 eigenvector, v1, in which case R = [〈v̂1,Mv̂1〉] is a 1× 1 matrix. In this setting, Equation (3) is not applicable because there are no off-diagonal elements, so maxv̂>1 v̂1=1〈v̂1,Mv̂1〉 is a sensible utility function for v̂1.
If considering the top-2 eigenvectors, v̂1’s utility remains as before, and we introduce a new utility for v̂2. Equation (3) is now applicable, so v̂2’s utility is
max v̂>2 v̂2=1,v̂ > 1 v̂2=0
〈v̂2,Mv̂2〉 − 〈v̂2,Mv̂1〉2
〈v̂1,Mv̂1〉 (5)
where we have divided the off-diagonal penalty by 〈v1,Mv1〉 so a) the two terms in Equation (5) are on a similar scale and b) for reasons that ease analysis. Additionally note that the constraint v̂>1 v̂2 = 0 may be redundant at the optimum (v̂∗1 = v1, v̂ ∗ 2 = v2) because the second term, 〈v̂∗2 ,Mv̂∗1〉2 = 〈v2,Mv1〉2 = Λ211〈v2, v1〉2, already penalizes such deviations (Λii is the ith largest eigenvector). These reasons motivate the following set of objectives (utilities), one for each vector i ∈ {1, . . . , k}:
max v̂>i v̂i=1
{ ui(v̂i|v̂j<i) = v̂>i Mv̂i − ∑ j<i (v̂>i Mv̂j) 2 v̂>j Mv̂j = ||Xv̂i||2 − ∑ j<i 〈Xv̂i, Xv̂j〉2 〈Xv̂j , Xv̂j〉 } (6)
where the notation ui(ai|b) emphasizes that player i adjusts ai to maximize a utility conditioned on b. It is interesting to note that by incorporating knowledge of the natural hierarchy (see Figure 1), we are immediately led to constructing asymmetric utilities, and thereby, inspired to formulate the PCA problem as a game, rather than a direct optimization problem as in Equation (4).
A key concept in games is a Nash equilibrium. A Nash equilibrium specifies a variable for each player from which no player can unilaterally deviate and improve their outcome. In this case, V̂ is a (strict-)Nash equilibrium if and only if for all i, ui(v̂i|v̂j<i) > ui(zi|v̂j<i) for all zi ∈ Sd−1. Theorem 2.1 (PCA Solution is the Unique strict-Nash Equilibrium). Assume that the top-k eigenvalues of X>X are positive and distinct. Then the top-k eigenvectors form the unique strictNash equilibrium of the proposed game in Equation (6).3 The proof is deferred to Appendix L.
Solving for the Nash of a game is difficult in general. Specifically, it belongs to the class of PPADcomplete problems (Gilboa and Zemel, 1989; Daskalakis et al., 2009). However, because the game
is hierarchical and each player’s utility only depends on its parents, it is possible to construct a sequential algorithm that is convergent by solving each player’s optimization problem in sequence.

heading: 3 METHOD
text: Utility gradient. In Section 2, we mentioned that normalizing the penalty term from Equation (5) had a motivation beyond scaling. Dividing by 〈v̂j ,Mv̂j〉 results in the following gradient for player i:
∇v̂iui(v̂i|v̂j<i) = 2M [ v̂i − ∑ j<i v̂>i Mv̂j v̂>j Mv̂j
v̂j︸ ︷︷ ︸ generalized Gram-Schmidt
] = 2X> [ Xv̂i − ∑ j<i 〈Xv̂i, Xv̂j〉 〈Xv̂j , Xv̂j〉 Xv̂j ] . (7)
The resulting gradient with normalized penalty term has an intuitive meaning. It consists of a single generalized Gram-Schmidt step followed by the standard matrix product found in power iteration and Oja’s rule. Also, notice that applying the gradient as a fixed point operator in sequence (v̂i ← 1 2∇v̂iui(v̂i|v̂j<i)) on M = I recovers the standard Gram-Schmidt procedure for orthogonalization.
A sequential algorithm. Each eigenvector can be learned by maximizing its utility. The vectors are constrained to the unit sphere, a non-convex Riemannian manifold, so we use Riemmanian gradient ascent with gradients given by Equation (7). In this case, Riemannian optimization theory simply requires an intermediate step where the gradient,∇v̂i , is projected onto the tangent space of the sphere to compute the Riemannian gradient, ∇Rv̂i . A more detailed illustration can be found in Appendix J. Recall that each ui depends on v̂j<i. If any of v̂j<i are being learned concurrently, then v̂i is maximizing a non-stationary objective which makes a convergence proof difficult. Instead, for completeness, we prove convergence assuming each v̂i is
learned in sequence. Algorithm 1 learns v̂i given fixed parents v̂j<i; we present the convergence guarantee in Section 4 and details on setting ρi and α in Appendix O.
Algorithm 1 EigenGameR-Sequential Given: matrix X ∈ Rn×d, maximum error tolerance ρi, initial vector v̂0i ∈ Sd−1, learned approximate parents v̂j<i, and step size α. v̂i ← v̂0i ti = d 54 min(||∇v̂0i ui||/2, ρi)
−2e for t = 1 : ti do
rewards← Xv̂i penalties← ∑ j<i 〈Xv̂i,Xv̂j〉 〈Xv̂j ,Xv̂j〉Xv̂j
∇v̂i ← 2X> [ rewards− penalties ] ∇Rv̂i ← ∇v̂i − 〈∇v̂i , v̂i〉v̂i v̂′i ← v̂i + α∇Rv̂i v̂i ← v̂ ′ i
||v̂′i|| end for return v̂i
Algorithm 2 EigenGameR (EigenGame—update with∇v̂i instead of∇Rv̂i )
Given: stream, Xt ∈ Rm×d, total iterations T , initial vector v̂0i ∈ Sd−1, and step size α. v̂i ← v̂0i for t = 1 : T do
rewards← Xtv̂i penalties← ∑ j<i 〈Xtv̂i,Xtv̂j〉 〈Xtv̂j ,Xtv̂j〉Xtv̂j
∇v̂i ← 2X>t [ rewards− penalties ] ∇Rv̂i ← ∇v̂i − 〈∇v̂i , v̂i〉v̂i v̂′i ← v̂i + α∇Rv̂i v̂i ← v̂ ′ i ||v̂′i|| broadcast(v̂i)
end for return v̂i
A decentralized algorithm. While Algorithm 1 enjoys a convergence guarantee, learning every parent v̂j<i before learning v̂i may be unnecessarily restrictive. Intuitively, as parents approach their respective optima, they become quasi-stationary, so we do not expect maximizing utilities in parallel to be problematic in practice. To this end, we propose Algorithm 2 visualized in Figure 2.
3Unique up to a sign change; this is expected as both vi and −vi represent the same principal component.
In practice we can assign each eigenvector update to its own device (e.g. a GPU or TPU). Systems with fast interconnects may facilitate tens, hundreds or thousands of accelerators to be used. In such settings, the overhead of broadcast(v̂i) is minimal. We can also specify that the data stream is co-located with the update so v̂i updates with respect to its own Xi,t. This is a standard paradigm for e.g. data-parallel distributed neural network training. We provide further details in Section 6.
Message Passing on a DAG. Our proposed utilities enforce a strict hierarchy on the eigenvectors. This is a simplification that both eases analysis (see Appendix M) and improves convergence4, however, it is not optimal. We assume vectors are initialized randomly on the sphere and, for instance, v̂k may be initialized closer to v1 than even v̂1 and vice versa. The hierarchy shown in Figure 1 enforces a strict graph structure for broadcasting information of parents to the childrens’ utilities.
To our knowledge, our utility formulation in Equation (6) is novel. One disadvantage is that stochastic gradients of Equation (7) are biased. This is mitigated with large batch sizes (further discussion in Appendix I).

heading: 4 CONVERGENCE OF EIGENGAME
text: Here, we first show that Equation (6) has a simple form such that any local maximum of ui is also a global maximum. Player i’s utility depends on its parents, so we next explain how error in the parents propagates to children through mis-specification of player i’s utility. Using the first result and accounting for this error, we are then able to give global, finite-sample convergence guarantees in the full-batch setting by leveraging recent non-convex Riemannian optimization theory.
The utility landscape and parent-to-child error propagation. Equation (6) is abstruse, but we prove that the shape of player i’s utility is simply sinusoidal in the angular deviation of v̂i from the optimum. The amplitude of the sinusoid varies with the direction of the angular deviation along the unit-sphere and is dependent on the accuracy of players j < i. In the special case where players j < i have learned the top-(i− 1) eigenvectors exactly, player i’s utility simplifies (see Lemma N.1) to
ui(v̂i, {vj<i}) = Λii − sin2(θi) ( Λii − ∑ l>i zlΛll ) (8)
where θi is the angular deviation and z ∈ ∆d−1 parameterizes the deviation direction. Note that sin2 has period π instead of 2π, which simply reflects the fact that vi and −vi are both eigenvectors.
4EigenGame sans order learns max 1 PC and sans order+normalization 5 PCs on data in Figure 3a. 5EigenGame runtimes are longer than those of EigenGameR in the synthetic experiments despite strictly requiring fewer FLOPS; apparently this is due to low-level floating point arithmetic specific to the experiments.
An error propagation analysis reveals that it is critical to learn the parents to a given degree of accuracy. The angular distance between vi and the maximizer of player i’s utility with approximate parents has tan−1 dependence (i.e., a soft step-function; see Lemma N.5 and Figure 13 in Appendix N). Theorem 4.1 (Global convergence). Algorithm 1 achieves finite sample convergence to within θtol angular error of the top-k principal components, independent of initialization. Furthermore, if each v̂i is initialized to within π4 of vi, Algorithm 1 returns the components with angular error less than θtol
in T = ⌈ O ( k [
(k−1)! θtol k∏ i=1 ( 16Λ11 gi )]2)⌉ iterations. Proofs are deferred to Appendices O.4 and O.5.
Angular error is defined as the angle between v̂i and vi: θi = sin−1( √
1− 〈vi, v̂i〉2). The first k in the formula for T appears from a naive summing of worst case bounds on the number of iterations required to learn each v̂j<k individually. The constant 16 arises from the error propagation analysis; parent vectors, v̂j<i, must be learned to under 1/16th of a canonical error threshold, gi(i−1)Λ11 , for the child v̂i where gi = Λii − Λi+1,i+1. The Riemannian optimization theory we leverage dictates that 1 ρ2 iterations are required to meet a O(ρ) error threshold. This is why the squared inverse of the error threshold appears here. Breaking down the error threshold itself, the ratio Λ11/gi says that more iterations are required to distinguish eigenvectors when the difference between them (summarized by the gap gi) is small relative to the scale of the spectrum, Λ11. The (k − 1)! term appears because learning smaller eigenvectors requires learning a much more accurate v̂1 higher up the DAG.
Lastly, the utility function for each v̂i is sinusoidal, and it is possible that we initialize v̂i with initial utility arbitrarily close to the trough (bottom) of the function where gradients are arbitrarily small. This is why the global convergence rate depends on the initialization in general. Note that Algorithm 1 effectively detects the trough by measuring the norm of the initial gradient (∇v̂0i ui) and scales the number of required iterations appropriately. A complete theorem that considers the probability of initializing v̂i within π4 of vi is in Appendix O, but this possibility shrinks to zero in high dimensions.
We would also like to highlight that these theoretical findings are strong relative to some other claims. For example, the exponential convergence guarantee for Matrix Krasulina requires the initial guess at the eigenvectors capture the top-(k − 1) subspace (Tang, 2019), unlikely when d k. A similar condition is required in (Shamir, 2016b). These guarantees are given for the mini-batch setting while ours is for the full-batch, however, we provide global convergence without restrictions on initialization.

heading: 5 RELATED WORK
text: PCA is a century-old problem and a massive literature exists (Jolliffe, 2002; Golub and Van Loan, 2012). The standard solution to this problem is to compute the SVD, possibly combined with randomized algorithms, to recover the top-k components as in (Halko et al., 2011) or with Frequent Directions (Ghashami et al., 2016) which combines sketching with SVD.
In neuroscience, Hebb’s rule (Hebb, 2005) refers to a connectionist rule that solves for the top eigenvector of a matrix M using additive updates of a vector v as v ← v + ηMv. Likewise, Oja’s rule (Oja, 1982; Shamir, 2015) refers to a similar update v ← v + η(I − vv>)Mv. In machine learning, using a normalization step of v ← v/||v|| with Hebb’s rule is somewhat confusingly referred to as Oja’s algorithm (Shamir, 2015), the reason being that the subtractive term in Oja’s rule can be viewed as a regularization term for implicitly enforcing the normalization. In the limit of infinite step size, η →∞, Oja’s algorithm effectively becomes the well known Power method. If a normalization step is added to Oja’s rule, this is referred to as Krasulina’s algorithm (Krasulina, 1969). In the language of Riemannian manifolds, v/||v|| can be recognized as a retraction and (I − vv>) as projecting the gradient Mv onto the tangent space of the sphere (Absil et al., 2009).
Many of the methods above have been generalized to the top-k components. Most generalizations involve adding an orthonormalization step after each update, typically accomplished with a QR factorization plus some minor sign accounting (e.g., see Algorithm 3 in Appendix A.1). An extension of Krasulina’s algorithm to the top-k setting, termed Matrix Krasulina (Tang, 2019), was recently proposed in the machine learning literature. This algorithm can be recognized as projecting the gradient onto the Stiefel manifold (the space of orthonormal matrices) followed by a QR step to maintain orthonormality, which is a well known retraction.
Maintaining orthonormality via QR is computationally expensive. Amid and Warmuth (2019) propose an alternative Krasulina method which does not require re-orthonormalization but instead requires inverting a k × k matrix; in a streaming setting restricted to minibatches of size 1 (Xt ∈ Rd), Sherman-Morrison (Golub and Van Loan, 2012) can be used to efficiently replace the inversion step. Raja and Bajwa (2020) develop a data-parallel distributed algorithm for the top eigenvector. Alternatively, the Jacobi eigenvalue algorithm explicitly represents the matrix of eigenvectors as a Givens rotation matrix using sin’s and cos’s and rotates M until it is diagonal (Golub and Van der Vorst, 2000).
In contrast, other methods extract the top components in sequence by solving for the ith component using an algorithm such as power iteration or Oja’s, and then enforcing orthogonality by removing the learned subspace from the matrix, a process known as deflation. Alternatively, the deflation process may be intertwined with the learning of the top components. The generalized Hebbian algorithm (Sanger, 1989) (GHA) works this way as do Lagrangian inspired formulations (Ghojogh et al., 2019) as well as our own approach. We make the connection between GHA and our algorithm concrete in Prop. K.1. Note, however, that the GHA update is not the gradient of any utility (Prop. K.2) and therefore, lacks a clear game interpretation.
Of these, Oja’s algorithm has arguably been the most extensively studied (Shamir, 2016a; Allen-Zhu and Li, 2017)6 Note that Oja’s algorithm converges to the actual principal components (Allen-Zhu and Li, 2017) and Matrix Krasulina (Tang, 2019) converges to the top-k subspace. However, neither can be obviously decentralized. GHA (Sanger, 1989) converges to the principal components asymptotically and can be decentralized (Gang et al., 2019). Each of these is applicable in the streaming k-PCA setting.

heading: 6 EXPERIMENTS
text: We compare our approach against GHA, Matrix Krasulina, and Oja’s algorithm7. We present both EigenGame and EigenGameR which projects the gradient onto the tangent space of the sphere each step. We measure performance of methods in terms of principal component accuracy and subspace distance. We measure principal component accuracy by the number of consecutive components, or longest streak, that are estimated within an angle of π8 from ground truth. For example, if the angular errors of the v̂i’s returned by a method are, in order, [θ1, θ2, θ3, . . .] = [ π16 , π 4 , π 10 , . . .], then the method is credited with a streak of only 1 regardless of the errors θi>2. For Matrix Krasulina, we first compute the optimal matching from v̂i to ground truth before measuring angular error. We present the longest streak as opposed to “# of eigenvectors found” because, in practice, no ground truth is available and we think the user should be able to place higher confidence in the larger eigenvectors being correct. If an algorithm returns k vectors, k2 of which are accurate components but does not indicate which, this is less helpful. We measure normalized subspace distance using 1− 1k · Tr(U ∗P ) ∈ [0, 1] where U∗ = V V † and P = V̂ V̂ † similarly to Tang (2019).
Synthetic data. Experiments on synthetic data demonstrate the viability of our approach (Figure 3a). Oja’s algorithm performs best on synthetic experiments because strictly enforcing orthogonalization with an expensive QR step greatly helps when solving for all eigenvectors. EigenGame is able to effectively parallelize this over k machines and the advantage of QR diminishes in Figure 3b. The remaining algorithms perform similarly on a linearly decaying spectrum, however, EigenGame performs better on an exponentially decaying spectrum due possibly to instability of Riemannian gradients near the equilibrium (see Appendix J for further discussion). GHA and EigenGameR are equivalent under specific conditions (see Proposition K.1).
Figure 4a shows EigenGame solves for the eigenvectors up to a high degree of accuracy π32 , i.e. the convergence results in Figure 3a are not the result of using a loose tolerance of π8 . With the lower tolerance, all algorithms take slightly more iterations to learn the eigenvectors of the linear spectrum; it is difficult to see any performance change for the exponential spectrum. Although Theorem 4.1 assumes distinct eigenvalues, Figure 4b supports the claim that EigenGame does not require distinct eigenvalues for convergence. We leave proving convergence in this setting to future work.
6See Table 1 in (Allen-Zhu and Li, 2017). 7A detailed discussion of Frequent Directions (Ghashami et al., 2016) can be found in Appendix H.
MNIST handwritten digits. We compare EigenGame against GHA, Matrix Krasulina, and Oja’s algorithm on the MNIST dataset (Figure 3b). We flatten each image in the training set to obtain a 60, 000× 784 dimensional matrix. EigenGame is competitive with Oja’s in a high batch size regime (1024 samples per mini-batch). The performance gap between EigenGame and the other methods shrinks as the mini-batch size is reduced (see Appendix I), expectedly due to biased gradients.
The principal components of RESNET-200 activations on IMAGENET are edge filters. A primary goal of PCA is to obtain interpretable low-dimensional representations. To this end we present an example of using EigenGame to compute the top-32 principal components of the activations of a pretrained RESNET-200 on the IMAGENET dataset. We concatenate the flattened activations from the output of each residual block resulting in a d ≈ 20M dimensional vector representation for each of the roughly 1.2M input images. It is not possible to store the entire 195TB matrix in memory, nor incrementally compute the Gram/covariance matrix.
We implemented a data-and-model parallel version of EigenGame in JAX (Bradbury et al., 2018) where each v̂i is assigned to it’s own TPU (Jouppi et al., 2017). Each device keeps a local copy of the RESNET parameters and the IMAGENET datastream. Sampling a mini-batch (of size 128), computing the network activations and updating v̂i are all performed locally. The broadcast(v̂i) step is handled by the pmap and lax.all_gather functions. Computing the top-32 principal components takes approximately nine hours on 32 TPUv3s.
Figure 5a shows the top principal components of the activations of the trained network organized by network topology (consisting of five residual blocks). Note that EigenGame is not applied block-wise, but on all 20M dimensions. We do not assume independence between blocks and the eigenvector has unit norm across all blocks. We observe that Block 1 (closest to input) of PC 1 has very small magnitude activations relative to the other PCs. This is because PC 1 should capture the variance which discriminates most between the classes in the dataset. Since Block 1 is mainly
concerned with learning low-level image filters, it stands to reason that although these are important for good performance, they do not necessarily extract abstract representations which are useful for classification. Conversely, we see that PC 1 has larger relative activations in the later blocks.
We visualize the average principal activation in Block 18 in Figure 5b. The higher PCs learn distinct filters (Gabor filters, Laplacian-of-Gaussian filters c.f. (Bell and Sejnowski, 1997)).
7 CONCLUSION It seems easier to train a bi-directional LSTM with attention than to compute the SVD of a large matrix. –Chris Re NeurIPS 2017 Test-of-Time Award, Rahimi and Recht
(Rahimi and Recht, 2017).
In this work we motivated PCA from the perspective of a multi-player game. This inspired a decentralized algorithm which enables large-scale principal components estimation. To demonstrate this we used EigenGame to analyze a large neural network through the lens of PCA. To our knowledge this is the first academic analysis of its type and scale (for reference, (Tang, 2019) compute the top-6 PCs of the d = 2300 outputs of VGG). EigenGame also opens a variety of research directions.
Scale. In experiments, we broadcast across all edges in Figure 1 every iteration. Introducing lag or broadcasting with dropout may improve efficiency. Can we further reduce our memory footprint by storing only scalars of the losses and avoiding congestion through online bandit or reinforcement learning techniques? Our decentralized algorithm may have implications for federated and privacy preserving learning as well (Heinze et al., 2016; Heinze-Deml et al., 2018; Bonawitz et al., 2019).
Games. EigenGame has a unique Nash equilibrium due to the fixed DAG structure, but vectors are initialized randomly so v̂k may start closer to v1 than v̂1 does. Adapting the DAG could make sense, but might also introduce spurious fixed points or suboptimal Nash. Might replacing vectors with populations accelerate extraction of the top principal components?
Core ML. EigenGame could be useful as a diagnostic or for accelerating training (Desjardins et al., 2015; Krummenacher et al., 2016); similarly, spectral normalization has shown to be a valuable tool for stabilizing GAN training (Miyato et al., 2018).
Lastly, GANs (Goodfellow et al., 2014) recently reformulated learning a generative model as a two-player zero-sum game. Here, we show how another fundamental unsupervised learning task can be formulated as a k-player game. While two-player, zero-sum games are well understood, research on k-player, general-sum games lies at the forefront in machine learning. We hope that marrying a fundamental, well-understood task in PCA with the relatively less understood domain of many player games will help advance techniques on both ends.

heading: ACKNOWLEDGEMENTS
text: We are grateful to Trevor Cai for his help scaling the JAX implementation of EigenGame to handle the large IMAGENET experiment and to Daniele Calandriello for sharing his expert knowledge of related work and advice on revising parts of the manuscript.

Review summaries:
Principal component analysis (PCA) is a well-known dimensionality reduction and feature learning technique in the literature that leads to uncorrelated features. While there are a plethora of algorithms for PCA, along with accompanying analysis, a majority of these works have been developed from an optimization perspective. This paper differs from existing works in that it motivates the $k$-PCA problem, which involves learning the $k$-dominant eigen vectors of the sample covariance matrix, as a competitive game between $k$ players in which each player is supposed to estimate one of the eigen vectors and the PCA solution is the unique strict-Nash equilibrium. The main contributions of the paper in this regard are the following:

The authors present new insights on PCA analysis by reconceiving it in terms of a Nash equilibrium among different players, related to the different components. The importance of an objective function minimizing the off-diagonal elements of R is emphasized. The insights lead to parallel algorithms and are demonstrated on large scale problems, which is nice. Overall the new insights can be very valuable and also inspiring for future work and for new developments,  from a broader perspective.