A Stochastic PCA and SVD Algorithm with an Exponential Convergence Rate

Ohad Shamir

Source-linked AI summary

A Stochastic PCA and SVD Algorithm with an Exponential Convergence Rate

Ohad Shamir

arXiv:1409.2848v5Sep 9, 2014cs.LGmath.NAmath.OCstat.ML

TL;DR

Existing PCA and SVD methods trade off cheap iterations against slow convergence or expensive data-dependent computation. The paper introduces VR-PCA, which applies variance reduction to stochastic PCA updates and proves exponential convergence under suitable assumptions, with logarithmic accuracy dependence and a favorable eigengap/runtime relationship.

Problem
Existing PCA and SVD algorithms either converge slowly or use computationally intensive iterations whose runtime scales unfavorably with data size or accuracy.
Method
VR-PCA combines stochastic PCA updates with variance reduction so the stochastic variance decays as iterates approach the optimum.
Results
VR-PCA achieves exponential convergence, with runtime O(d_s n log(1/ϵ)) when λ/r ≥ Ω(1/√n), up to logarithmic factors.
Takeaways & Limitations
Under suitable assumptions, the method offers high-accuracy PCA and SVD using runtime proportional up to logarithmic factors to a single data scan.
Takeaways & Limitations
The formal analysis begins after a burn-in period reaching a constant distance from v1, and randomly initialized behavior remains unanalyzed.

Abstract

from arXiv · show

We describe and analyze a simple algorithm for principal component analysis and singular value decomposition, VR-PCA, which uses computationally cheap stochastic iterations, yet converges exponentially fast to the optimal solution. In contrast, existing algorithms suffer either from slow convergence, or computationally intensive iterations whose runtime scales with the data size. The algorithm builds on a recent variance-reduced stochastic gradient technique, which was previously analyzed for strongly convex optimization, whereas here we apply it to an inherently non-convex problem, using a very different analysis.

1 Introduction

The paper studies recovering leading singular vectors for PCA and SVD at scale, where exact decomposition is costly and existing iterative methods face runtime or convergence limitations. It proposes VR-PCA, a stochastic method intended to combine cheap iterations with favorable accuracy dependence and eigengap scaling.

Problem setting: PCA seeks a k-dimensional subspace capturing maximum data variance, equivalently the top k left singular vectors of X.The paper mainly studies k = 1, while noting that the algorithm can extend to k > 1.
Problem setting: PCA supports dimensionality reduction, data compression, and data visualization, making efficient large-scale computation practically relevant.
Computational challenge: Full SVD is prohibitive for large matrices, while power and related iterations require repeated covariance-matrix operations over the data.A covariance-matrix multiplication can require O(d_s n) time when d_s is average input sparsity.
Computational challenge: Deterministic iterative methods can require many data passes when the eigengap λ is small, whereas stochastic methods use cheaper O(d_s) iterations but generally converge slowly with iteration counts scaling linearly in 1/ϵ.
Proposed approach: VR-PCA uses stochastic iterations with runtime logarithmic in accuracy ϵ and scaling as a sum of n and an eigengap-dependent factor rather than their product.For λ ≥ Ω(1/√n), the bound is proportional up to logarithmic factors to one data scan, d_s n.
Proposed approach: The method applies variance reduction to an inherently non-convex, concave optimization problem, requiring an analysis different from prior strongly convex applications.

2 Algorithm and Analysis

VR-PCA organizes stochastic power iterations into epochs, using variance reduction to make the stochastic error decay as iterates approach the solution. Under eigengap and initialization assumptions, the analysis establishes high-probability exponential convergence with runtime logarithmic in accuracy.

Algorithm: VR-PCA runs m stochastic iterations per epoch, sampling data points uniformly and starting each epoch from the previous epoch's vector.
Algorithm: Each iteration combines a power-iteration update with a stochastic zero-mean term whose variance is controlled by the distance between the current vector and the epoch-start vector.The variance is expected to decay as the iterates approach the optimum.
Extensions: VR-PCA extends to multiple singular vectors through orthogonalization or matrix deflation, but deflation requires a positive eigengap between all top k singular values.
Analysis: The theorem assumes a positive leading eigengap, suitable step size and epoch length, and an initialization with constant alignment to the top eigenvector.The stated probability guarantee is at least 1 − 2 log(1/ϵ)δ.
Limitations: The initialization condition requires a separate burn-in analysis because a random unit vector typically has only O(1/√d) alignment with the top eigenvector.The paper leaves this burn-in analysis to future work, although experiments did not visibly show a burn-in period.

3 Experiments

Experiments compare VR-PCA with power iterations, Oja’s algorithm, and a hybrid method on synthetic, MNIST, and CCAT data. VR-PCA exhibits exponential convergence and substantially faster convergence than the tested alternatives, while multi-vector recovery is qualitatively similar to the single-vector case.

Experimental setup: The experiments use epoch length m = n and a data-dependent step size η = 1/(r̄√n), avoiding parameter tuning and requiring no knowledge of λ.This choice balances stochastic updates with computation of the epoch reference vector.
Synthetic data: The experiments use synthetic matrices with n = 200000 and d = 10000, varying the eigengap λ and initializing all algorithms from the same random vector.The synthetic spectrum is constructed through a diagonal matrix with controlled leading eigengap.
Synthetic data: VR-PCA converges much faster than Oja’s algorithm and power iterations across all tested synthetic eigengaps, despite untuned parameters.The logarithmic y-axis indicates an exponential convergence rate for VR-PCA, whereas Oja’s convergence appears sub-exponential.
Real data: The MNIST and CCAT experiments compare the same algorithms with a hybrid method that initializes VR-PCA using n iterations of Oja’s algorithm.The hybrid uses Oja’s decaying step size for initialization and can outperform either algorithm alone.
Multiple singular vectors: The multi-vector experiments recover k = 6 singular vectors for MNIST and k = 3 for CCAT, with results qualitatively similar to the k = 1 case.The competitors use analogous generalizations for multiple singular vectors.

4 Proof of Thm. 1

The proof normalizes the data and analyzes VR-PCA under a bounded-norm assumption with a positive eigengap. Its conventions fix the spectral setting and simplify constants used throughout the argument.

Normalization: The analysis rescales the data so that max_i ∥x_i∥^2 ≤ 1, replacing the original eigengap and step size by their rescaled equivalents.The final result is recovered by substituting λ/r and ηr for the normalized parameters.
Spectral assumptions: The matrix A is assumed to have eigenvalues s1 > s2 ≥ . . . ≥ sd with eigengap s1 − s2 = λ > 0.Under normalization, all eigenvalues lie in [0, 1].
Proof conventions: The proof uses c for positive numerical constants whose values may change from one occurrence to another.This convention keeps the bounds readable without tracking each constant explicitly.

Part I: Establishing a Stochastic Recurrence Relation

The first proof part derives a one-step recurrence for the squared alignment error with the top eigenvector. It controls the stochastic perturbation through the distance between the current iterate and the epoch reference vector.

Error measure: The analysis tracks how 1 − ⟨w_t, v1⟩^2 changes during one iteration and one epoch.This quantity measures the iterate’s misalignment with the leading eigenvector.
Stochastic decomposition: The stochastic update is decomposed into fixed coefficients and zero-mean random terms induced by the sampled data point.Conditioned on w_t, the coefficients are fixed while the stochastic terms vary with the sampled x.
Expectation bounds: A second-order Taylor expansion relates the expected ratio defining the alignment to its values at expected numerator and denominator.Bounds on the Hessian and fluctuation terms control the approximation error.
Variance control: The perturbation terms are bounded by quantities proportional to η^2∥w_t − ẇ∥^2, linking stochastic variance to the iterate’s distance from the epoch reference.The argument uses bounded norms, orthonormal eigenvectors, and sufficiently small η.
One-step recurrence: For sufficiently small η ≤ cλ, the resulting lower bound yields a contraction term governed by the eigengap λ.The proof obtains the recurrence by substituting successive bounds and subtracting from one.

Part II: Solving the Recurrence Relation for a Single Epoch

The second proof part solves the single-epoch recurrence by controlling the stochastic process with concentration inequalities. With suitable initialization and parameters, the iterate remains in the analyzed region and preserves nonnegative alignment with the top eigenvector.

Recurrence process: The proof sets η = αλ and studies the process b_t representing the alignment error throughout an epoch.The step size is tied directly to the eigengap through a sufficiently small constant α.
Invariant region: The recurrence analysis is conditioned on the event that b_t remains at most 3/4 throughout the epoch.This event ensures the iterate stays within the region required by the one-step bounds.
High-probability control: A maximal Hoeffding-Azuma argument bounds the process simultaneously over all iterations with probability at least 1 − β.The resulting bound controls both the error process and the final alignment sign.
Recurrence solution: The conditional expectation of b_{t+1} contracts b_t up to an additive term proportional to α^2λ^2 times the epoch-start error.The bounded increments support concentration around this recurrence.
Epoch guarantee: Combining the recurrence lemmas and Markov’s inequality gives an epoch-level guarantee with failure probability at most β + γ.The conclusion also preserves nonnegative alignment with v1 at the epoch’s end.

Part III: Analyzing the Entire Algorithm’s Run

The analysis repeatedly contracts the squared error relative to the top eigenvector across epochs, yielding exponential convergence and logarithmic dependence on the target accuracy.

Theorem guarantee: The proof begins from an iterate whose squared error is at most 1/2 and propagates the contraction while maintaining the conditions needed by the single-epoch lemma.The initial bound is established from the assumed constant alignment with v1.
Epoch-wise convergence: The squared error contracts by a factor γ per epoch, so after T epochs it is at most γ^T with probability at least 1 − T(β + γ).The analysis obtains this by repeatedly applying the single-epoch lemma and a union bound.
Epoch-wise convergence: For target accuracy ϵ, the algorithm needs a logarithmic number of epochs in 1/ϵ.This follows by choosing T so that γ^T is below ϵ.
Theorem guarantee: The resulting accuracy guarantee is ⟨w̃_T, v1⟩^2 ≥ 1 − ϵ with probability at least 1 − 2 log(1/ϵ)δ under the stated parameter conditions.The theorem substitutes the confidence parameters and the step-size relation η = αλ into the preceding bounds.

5 Discussion

The discussion highlights VR-PCA’s logarithmic accuracy dependence and additive data-size/eigengap runtime, while identifying regimes and initialization questions that remain unresolved.

Main discussion: VR-PCA’s runtime scales as n plus an eigengap factor 1/λ^2 and logarithmically in the required accuracy ϵ.The paper contrasts this with prior iterative methods whose runtime scales either multiplicatively with n and the eigengap factor or polynomially with ϵ.
Open limitations: When n is dominated by (r/λ)^2, the runtime’s 1/λ^2 dependence is inferior to the deterministic methods discussed earlier.This identifies a moderate-data-size regime where the stochastic method is not preferable by the paper’s runtime comparison.
Open limitations: The analysis does not establish whether the 1/λ^2 factor is necessary or can be improved to 1/λ.The discussion motivates the conjecture by analogy with variance-reduced methods for strongly convex optimization but leaves it unresolved.
Open limitations: A full analysis from random initialization remains open, although experiments do not show a discernible burn-in period.The theorem instead analyzes convergence after an initial iterate has reached a constant distance from the optimum.

A Implementing Epochs in O(ds(m + n)) Amortized Runtime

The implementation represents iterates through sparse updates and auxiliary state, reducing the amortized epoch cost to O(ds(m + n)) while preserving normalization.

Runtime: An epoch costs O(ds(m + n)) by combining an O(dsn) full-data computation with m iterations implemented in O(ds) amortized time each.The sparse-data implementation improves the direct O(dm + dsn) accounting to the stated bound.
Epoch structure: The algorithm initializes w0 from the previous epoch’s normalized iterate and samples each index it uniformly at random during the inner loop.These operations appear in the epoch pseudocode and support the stochastic update structure.
Sparse representation: Instead of storing wt and w′t explicitly, the implementation stores an auxiliary vector g, the full-data vector ũ, and scalars α, β, γ, δ, ζ.This representation allows sparse updates without reconstructing the dense iterates at every iteration.
Sparse representation: The stored state represents w′t as αg + βũ, while γ + 2δ + ζ tracks its squared norm for normalization.The scalar invariants preserve the norm needed by the normalization step.
Sparse representation: Initialization and normalization take O(d) time, while the update maintains the representation using sparse operations in O(ds) time.The construction assumes d ≤ dsn after dropping all-zero columns when necessary.