Subspace Power Method (SPM)

• SPM is an iterative algorithm family that extracts low-dimensional subspaces using updates like QR and thin SVD, enabling efficient optimization of high-dimensional algebraic objects.
• It is applied in subspace optimization, eigenvalue problems, trace-ratio optimization, and tensor decompositions, with convergence guaranteed under specific spectral or geometric gap conditions.
• SPM achieves low per-iteration complexity (O(nk²)) by leveraging matrix-free techniques and block-orthogonal updates, making it suitable for scalable data analysis and scientific computing.

The Subspace Power Method (SPM) is a broad family of iterative algorithms for efficiently optimizing or decomposing high-dimensional algebraic objects via subspace-based updates. SPM generalizes and unifies classical power methods in several research areas, including subspace optimization, eigenvalue problems, large-scale trace-ratio optimization, and tensor decompositions. Its hallmark is the representation and extraction of low-dimensional subspaces through iterative projection, block-orthogonalization, and power-like iterations, typically paired with efficient matrix factorizations instead of full spectral decompositions. SPM achieves favorable convergence rates and computational complexity under certain spectral or geometric gap conditions, and has become a foundational technique in scalable data analysis, scientific computing, and signal processing.

1. Fundamental Principles and Algorithmic Structure

The unifying theme in SPM is the iterative update of a low-dimensional subspace (or its orthogonal basis), aiming to optimize a convexpolynomial, Rayleigh-type, or multilinear objective by alternating between subspace expansion, projection/extraction, and (optionally) deflation. Almost all SPM variants operate on objects of the form $X\in\mathbb{R}^{n\times k}$ with $X^\top X = I_k$, representing an orthonormal basis for the target $k$-dimensional subspace. The per-iteration cost scales as $O(nk^2)$ for a QR factorization or SVD of a thin matrix, in contrast to $O(n^3)$ for full-rank eigendecompositions.

A prototypical SPM iteration for subspace optimization is

$$\widetilde X_{t+1} = X_t - \eta \nabla f(X_t X_t^\top) X_t, \ (X_{t+1}, R_{t+1}) = \mathrm{QR}(\widetilde X_{t+1}),$$

with $f$ convex and smooth, and $X_{t+1}$ remaining orthonormal (Garber et al., 2022).

For generalized eigenvalue or trace-ratio problems, SPM maintains an evolving search subspace, projects the objective, solves the reduced problem, computes the residual, expands the subspace, and restarts as needed to preserve convergence and monotonicity (Ferrandi et al., 2024).

In tensor decomposition, SPM alternates between subspace extraction via flattening/SVD, and a power-type inner iteration (often blockwise), repeatedly seeking rank-1 summands from the extracted subspace (Kileel et al., 2019, Wang et al., 21 Oct 2025).

2. Subspace Power Method for Subspace Optimization

SPM plays a central role in minimizing convex, smooth loss functions over the Grassmann manifold of $k$-dimensional subspaces:

$$\min_{X\in \mathbb{R}^{n\times k}} f(X X^\top), \quad X^\top X = I_k,$$

where $X^\top X = I_k$0 is $X^\top X = I_k$1-smooth. By iteratively taking a gradient step on the Stiefel manifold, followed by QR orthonormalization, SPM obviates the need for full singular value decompositions:

$X^\top X = I_k$2

A local linear convergence result is established under a strict complementarity (eigengap) condition, which ensures that the optimizer to the convex relaxation over the Fantope,

$X^\top X = I_k$3

is unique and lies in the desired rank-$X^\top X = I_k$4 matrix set $X^\top X = I_k$5. Provided the initialization is within a specified Frobenius norm radius of the unique optimizer $X^\top X = I_k$6—which depends on the spectral gap $X^\top X = I_k$7 and smoothness $X^\top X = I_k$8—the method achieves a local contraction rate

$X^\top X = I_k$9

and per iteration cost is $k$0 (Garber et al., 2022).

The approach extends to projected gradient methods and Frank-Wolfe variants, operating either in matrix space or directly over the convex relaxation, provided iterates remain in a sufficiently small neighborhood of $k$1. Initialization is commonly performed using a PCA subspace or empirical covariance eigenvectors for data-driven problems.

3. SPM in Eigenvalue, Trace-Ratio, and Generalized Eigenproblems

SPM is adaptable to large-scale spectral optimization, such as trace-ratio problems and generalized eigenvalue problems. For $k$2 symmetric and $k$3 SPD, the goal is to find

$k$4

SPM organizes the iteration by:

• Maintaining a search subspace $k$5;
• Projecting $k$6 and $k$7 onto $k$8 and extracting the best $k$9-dimensional solution in the reduced basis;
• Computing the residual $O(nk^2)$0;
• Expanding $O(nk^2)$1 via the leading singular vectors of $O(nk^2)$2 and re-orthogonalizing;
• Restarting the process upon exceeding subspace dimension thresholds to preserve monotonicity.

This matrix-free method needs only $O(nk^2)$3- and $O(nk^2)$4-vector products, with per-iteration costs $O(nk^2)$5 and storage $O(nk^2)$6. Rigorous convergence analysis tracks the principal angles between the search subspace and the global maximizer; convergence is established under mild spectral separation (eigengap) conditions (Ferrandi et al., 2024).

Numerical experiments demonstrate that SPM variants ("TR-SP", "FDA-SP") yield comparable or superior accuracy to partial-Ritz and other eigen-solvers, with significantly reduced matrix-vector products and runtime.

4. Subspace Power Methods for Tensor Decomposition

SPM is extended to decomposing low-rank tensors in both symmetric and partially symmetric settings. For a real symmetric tensor $O(nk^2)$7 (even $O(nk^2)$8), SPM alternates between:

A. Flattening $O(nk^2)$9 into a symmetric matrix $O(n^3)$0 representing the Gram structure on order-$O(n^3)$1 symmetric products; B. Performing a thin eigendecomposition, yielding a subspace $O(n^3)$2; C. Applying a shifted symmetric higher-order power method (SS-HOPM), which iteratively maximizes $O(n^3)$3, alternating with deflation to isolate each rank-1 component in turn (Kileel et al., 2019).

A fundamental property of SPM's objective is that its maximum value is known in advance and achieved precisely at the rank-1 factors $O(n^3)$4, with no spurious maxima under genericity and low-rank conditions (Kileel et al., 2021). Deflation steps are performed via rank-1 updates of the eigenspace for efficiency.

The method generalizes further to non-symmetric and partially symmetric tensors ("Multi-Subspace Power Method") via subspace extraction to orthonormal-slice form (using SVD of matricizations), followed by block-coordinate power updates to compute partially symmetric singular vector tuples (pSVTs) corresponding bijectively to rank-1 summands, with global and local linear convergence guarantees under suitable shift parameters (Wang et al., 21 Oct 2025).

5. Energy-Norm Error Analysis and Multigrid Integration

SPM underpins block inverse iteration and multilevel solvers for high-dimensional eigenproblems, as established in energy-norm error theory (He et al., 2017). Given an SPD matrix $O(n^3)$5, its eigenproblem is projected onto successively enriched trial spaces. The SPM update cycles between enrichment by current Ritz vectors and Galerkin projection (Rayleigh–Ritz), followed by inverse iteration or multigrid-accelerated linear solves:

• Error in the $O(n^3)$6-energy norm is quantified using subspace projectors and spectral gap parameters, yielding explicit convergence factors.
• Block versions accelerate convergence for eigenvalue clusters, and the multigrid context reduces the dominant cost of large linear solves.
• Optimality in both energy and $O(n^3)$7 norms is achieved up to a factor controlled by the trial space quality and eigenvalue separation.

SPM thus unifies subspace projection, power iteration, and multigrid methods under a single analytic framework with rigorous convergence rates.

6. Practical Implementation and Numerical Performance

Across applications, SPM is characterized by:

• Low per-iteration complexity: Only thin QR or SVD factorization on $O(n^3)$8 matrices, or blockwise subspace updates, are required;
• Matrix-free variants: Especially in large-scale trace-ratio and eigenvalue problems, only matrix-vector (or tensor-vector) products are needed;
• Initialization dependence: For local linear convergence, initialization must be within a basin of attraction determined by spectral or curvature gaps; standard initializations (e.g., leading eigenvectors, PCA directions) are commonly employed;
• Robust deflation mechanisms: In tensor decomposition, efficient rank-1 removal is performed by updating the subspace representation in compressed coordinates;
• Strong numerical stability and competitive runtime: Empirical studies document significant acceleration and robustness to noise compared to classical ALS, NLS, or full eigen-decomposition approaches in tensor settings, and greatly reduced matrix-vector products in spectral methods (Kileel et al., 2019, Ferrandi et al., 2024, Wang et al., 21 Oct 2025).

7. Theoretical Guarantees and Extensions

SPM enjoys rigorous analytical support across settings:

• Local linear convergence under eigengap or strict complementarity assumptions in manifold-constrained problems (Garber et al., 2022);
• Exact recovery or near-global minima in symmetric and partially symmetric tensor settings under genericity and rank conditions, with explicit thresholds separating true and spurious critical points (Kileel et al., 2021);
• Explicit energy- and $O(n^3)$9-norm error bounds in eigenvalue estimation (He et al., 2017);
• Extensions to generalized power methods, Frank–Wolfe variants, and multigrid cycles.

A plausible implication is that SPM acts as a bridge between structure-preserving geometry (e.g., Grassmannians, symmetric spaces) and scalable iterative optimization, supporting advanced statistical, scientific, and engineering applications with rigorous guarantees and practical efficiency.