Algorithm and implementation notes

This vignette documents the algorithm behind mspca(), the implementation choices that make it fast, and how to set its parameters. Readers who want the theoretical analysis should consult Cory-Wright and Pauphilet (2026).

The problem

The goal of sparse PCA is to identify $r$ loading vectors $\boldsymbol{u}_1, \ldots, \boldsymbol{u}_r \in \mathbb{R}^p$ that collectively explain a large share of the variance in the data, while each vector involves only a small number of the $p$ original features. In the single-component case ($r = 1$) this is

$$\max_{\boldsymbol{u} \in \mathbb{R}^p} \ \boldsymbol{u}^\top \boldsymbol{\Sigma} \boldsymbol{u} \quad \text{s.t.} \quad \|\boldsymbol{u}\|_2 = 1, \ \|\boldsymbol{u}\|_0 \leq k,$$

for which many efficient algorithms exist (d’Aspremont et al. 2007; Yuan and Zhang 2013; Bertsimas et al. 2022). The challenge in the $r$-component case is coordinating the components so that they are non-redundant.

In standard dense PCA, non-redundancy is ensured by requiring the $r$ leading eigenvectors of $\boldsymbol{\Sigma}$ to be mutually orthogonal and their projections to be uncorrelated. These two properties hold simultaneously for eigenvectors but generally cannot both be enforced in the sparse setting. msPCA therefore supports either type of coupling constraint, selected with feasibilityConstraintType:

• Orthogonality (feasibilityConstraintType = 0, the default): the loading vectors are required to be mutually orthogonal, $\boldsymbol{u}_t^\top \boldsymbol{u}_{t'} = 0$ for all $t \neq t'$. This is the direct geometric analogue of standard PCA.
• Zero pairwise correlation (feasibilityConstraintType = 1): the projected components are required to be uncorrelated in the data, $\boldsymbol{u}_t^\top \boldsymbol{\Sigma} \boldsymbol{u}_{t'} = 0$ for all $t \neq t'$. This ensures each component captures statistically distinct information.

Writing $\boldsymbol{C} \in \{\mathbb{I}, \boldsymbol{\Sigma}\}$ to encode the constraint type, the $r$-component problem solved by msPCA is

$$\max_{\boldsymbol{U} \in \mathbb{R}^{p \times r}} \ \sum_{t=1}^r \boldsymbol{u}_t^\top \boldsymbol{\Sigma} \boldsymbol{u}_t \quad \text{s.t.} \quad \boldsymbol{u}_t^\top \boldsymbol{C} \boldsymbol{u}_{t'} = 0 \ \ \forall t \neq t', \quad \|\boldsymbol{u}_t\|_2 = 1, \ \|\boldsymbol{u}_t\|_0 \leq k_t \ \ \forall t \in [r]. \tag{1}$$

Orthogonality corresponds to $\boldsymbol{C} = \mathbb{I}$ and zero correlation to $\boldsymbol{C} = \boldsymbol{\Sigma}$. In the zero-correlation case the implementation uses $\boldsymbol{\Sigma}$ divided by the total variance $\mathrm{tr}(\boldsymbol{\Sigma})$. This leaves the feasible set unchanged and makes the violation measure invariant to a rescaling of the data.

The objective is the sum of per-component variances. Most approaches for sparse PCA with multiple PCs use this objective (Zou et al. 2006; Journée et al. 2010; Lu and Zhang 2012; Vu et al. 2013; Benidis et al. 2016; Cory-Wright and Pauphilet 2026). It corresponds to the variance of the orthogonal projection onto the span of $\boldsymbol{U}$ only when the loading vectors are orthogonal; in general it is the sum of the marginal variances of the sparse components.

Evaluation metrics

Variance explained. The cumulative fraction of total variance explained by $\boldsymbol{U} = [\boldsymbol{u}_1, \ldots, \boldsymbol{u}_r]$ is

$$\mathrm{FVE}(\boldsymbol{U}) = \frac{1}{\mathrm{tr}(\boldsymbol{\Sigma})} \sum_{t=1}^r \boldsymbol{u}_t^\top \boldsymbol{\Sigma} \boldsymbol{u}_t,$$

computed by fraction_variance_explained(). Because loading vectors may not be orthogonal, interpret this as a cumulative component-variance score rather than the variance of the orthogonal projection onto the span of $\boldsymbol{U}$. Per-component contributions are returned by fraction_variance_explained_perPC() and, unnormalized, by variance_explained_perPC().

Feasibility. The constraint violation measures how far the returned solution is from satisfying the coupling constraints. Under orthogonality it is

$$\mathrm{viol}_{\mathrm{off}}(\boldsymbol{U}) = \sum_{t > t'} \left| \boldsymbol{u}_t^\top \boldsymbol{u}_{t'} \right|,$$

and under zero pairwise correlation

$$\mathrm{viol}_{\mathrm{off}}(\boldsymbol{C}, \boldsymbol{U}) = \frac{1}{\mathrm{tr}(\boldsymbol{C})} \sum_{t > t'} \left| \boldsymbol{u}_t^\top \boldsymbol{C} \boldsymbol{u}_{t'} \right|,$$

both computed by feasibility_violation_off(). The second is normalized by the total variance $\mathrm{tr}(\boldsymbol{C})$: the loading vectors are unit-norm, so $|\boldsymbol{u}_t^\top \boldsymbol{C} \boldsymbol{u}_{t'}|$ is homogeneous of degree one in $\boldsymbol{C}$ and the unnormalized sum would depend on the units of the data. After normalization each pairwise term reads as a fraction of the total variance, hence is scale-invariant. The same convention is used inside the solver, and for the nonredundancy matrices stored on the fitted object.

Lagrangian alternating maximization

The key algorithmic idea is to handle the coupling constraints in (1) via a quadratic penalty in the objective. Introducing non-negative penalty parameters $\lambda_{t,t'}$ for each pair $t \neq t'$ (with $\lambda_{t,t'} = \lambda_{t',t}$) gives the penalized objective

$$\max_{\boldsymbol{U}} \ \sum_{t=1}^r \boldsymbol{u}_t^\top \boldsymbol{\Sigma} \boldsymbol{u}_t - \sum_{t \neq t'} \lambda_{t,t'} \left( \boldsymbol{u}_t^\top \boldsymbol{C} \boldsymbol{u}_{t'} \right)^2 \quad \text{s.t.} \quad \|\boldsymbol{u}_t\|_2 = 1, \ \|\boldsymbol{u}_t\|_0 \leq k_t \ \ \forall t. \tag{2}$$

For fixed components $\boldsymbol{u}_{t'}$, $t' \neq t$, and fixed penalties, the subproblem for $\boldsymbol{u}_t$ reduces to a non-convex single-component sparse PCA problem against the perturbed covariance matrix

$$\tilde{\boldsymbol{\Sigma}}_t = \boldsymbol{\Sigma} - \sum_{t' \neq t} \lambda_{t,t'} \boldsymbol{C} \boldsymbol{u}_{t'} \boldsymbol{u}_{t'}^\top \boldsymbol{C}. \tag{3}$$

This decomposition holds for both constraint types: with $\boldsymbol{C} = \mathbb{I}$ the perturbation is $\lambda_{t,t'} \boldsymbol{u}_{t'} \boldsymbol{u}_{t'}^\top$, and with $\boldsymbol{C} = \boldsymbol{\Sigma}$ it is $\lambda_{t,t'} \boldsymbol{\Sigma} \boldsymbol{u}_{t'} \boldsymbol{u}_{t'}^\top \boldsymbol{\Sigma}$.

Most methods for computing the leading sparse eigenvector require the matrix to be positive semidefinite. If $\tilde{\boldsymbol{\Sigma}}_t$ is not, we add a diagonal shift $\lambda_0 \mathbb{I}$, which does not change the optimal solution because all feasible vectors have unit norm. The shift used is $\lambda_0 = \sum_{t' \neq t} \lambda_{t,t'} \|\boldsymbol{C} \boldsymbol{u}_{t'}\|_2^2$, which bounds the deflation term by $\lambda_0 \|\boldsymbol{\beta}\|_2^2$ for every $\boldsymbol{\beta}$ and so guarantees $\tilde{\boldsymbol{\Sigma}}_t \succeq \boldsymbol{0}$. It is the smallest shift of this form, which preserves the eigengap the power method relies on, and it is recomputed at each inner step from the current components, so no estimate of the spectrum of $\boldsymbol{\Sigma}$ is needed. Under orthogonality $\|\boldsymbol{u}_{t'}\|_2 = 1$ and it reduces to $\sum_{t' \neq t} \lambda_{t,t'}$.

Iterating over $t = 1, \ldots, r$ and progressively increasing the penalties to drive constraint violations toward zero yields the following scheme.

Algorithm 1: Lagrangian alternating maximization for problem (1)

Require: covariance matrix $\boldsymbol{\Sigma}$, number of components $r$, sparsity budgets $k_1, \ldots, k_r$, constraint matrix $\boldsymbol{C} \in \{\mathbb{I}, \boldsymbol{\Sigma}\}$, iterations $L$, feasibility tolerance $\eta$

1. Initialize $\boldsymbol{u}_t^{(0)} \leftarrow \boldsymbol{0}$ for all $t \in [r]$; set $\lambda_{t,t'} \leftarrow 0$ for all $t \neq t'$
2. for $\ell = 1, \ldots, L$do
3.    for $t = 1, \ldots, r$do
4.      Compute $\tilde{\boldsymbol{\Sigma}}_t \leftarrow \boldsymbol{\Sigma} - \sum_{t' \neq t} \lambda_{t,t'} \boldsymbol{C} \boldsymbol{u}_{t'}^{(\ell-1)} \boldsymbol{u}_{t'}^{(\ell-1)\top} \boldsymbol{C} + \lambda_0 \mathbb{I}$
5.      Compute $\boldsymbol{u}_t^{(\ell)}$ via Algorithm 2 applied to $(\tilde{\boldsymbol{\Sigma}}_t, k_t)$
6.    end for
7.    if $\sum_t \left| \|\boldsymbol{u}_t^{(\ell)}\|^2 - 1 \right| + \sum_{t > t'} \left| \boldsymbol{u}_t^{(\ell)\top} \boldsymbol{C} \boldsymbol{u}_{t'}^{(\ell)} \right| \leq \eta$then
8.      Record $\{\boldsymbol{u}_t^{(\ell)}\}$ as feasible; update best solution if the objective improves
9.    end if
10.    Increase the $\lambda_{t,t'}$ values (see below)
11. end for
12. return best feasible solution found, or the last iterate if none was found

Each single-component subproblem is solved via the truncated power method (TPM, Yuan and Zhang 2013), which alternates between a power step (multiplying by $\tilde{\boldsymbol{\Sigma}}_t$) and a truncation step (retaining only the $k_t$ largest-magnitude entries). In practice TPM often finds high-quality solutions around two orders of magnitude faster than certifiably optimal methods (Berk and Bertsimas 2019; Behdin and Mazumder 2026). Each call starts from the current iterate and then draws random restarts, which guard against poor local optima.

Algorithm 2: truncated power method with random restarts (Yuan and Zhang 2013)

Require: matrix $\tilde{\boldsymbol{\Sigma}}$, sparsity budget $k$, iteration limit $L_{\mathrm{TPM}}$, time limit $T$

1. $\boldsymbol{u}_{\mathrm{best}} \leftarrow \boldsymbol{0}$
2. repeat
3.    Draw $\boldsymbol{u} \sim \mathcal{N}(\boldsymbol{0}, \mathbb{I})$
4.    repeat
5.      $\boldsymbol{u} \leftarrow \tilde{\boldsymbol{\Sigma}} \boldsymbol{u} / \|\tilde{\boldsymbol{\Sigma}} \boldsymbol{u}\|_2$   (power step)
6.      Zero out all but the $k$ entries of $\boldsymbol{u}$ largest in absolute value   (truncation step)
7.      $\boldsymbol{u} \leftarrow \boldsymbol{u} / \|\boldsymbol{u}\|_2$
8.    until $\boldsymbol{u}$ converges
9.    if $\boldsymbol{u}^\top \tilde{\boldsymbol{\Sigma}} \boldsymbol{u} > \boldsymbol{u}_{\mathrm{best}}^\top \tilde{\boldsymbol{\Sigma}} \boldsymbol{u}_{\mathrm{best}}$then $\boldsymbol{u}_{\mathrm{best}} \leftarrow \boldsymbol{u}$; reset iteration count
10. until time limit $T$ exceeded or no improvement after $L_{\mathrm{TPM}}$ iterations
11. return $\boldsymbol{u}_{\mathrm{best}}$

The scheme resembles an iterative deflation procedure (Mackey 2008) in which a single-component sparse PCA problem is solved against a surrogate matrix at each iteration. The key difference is that the deflated matrix is induced by an explicit penalty on the non-redundancy constraints, which is progressively increased throughout the algorithm.

Implementation details

Penalty update. The penalty parameters $\lambda_{t,t'}$ are initialized to zero and increased progressively across outer iterations, letting the algorithm explore freely at first and gradually tightening the feasibility requirement. During the first 15% of iterations the increment is proportional to the total constraint violation $\sum_{t > t'} |\boldsymbol{u}_t^\top \boldsymbol{C} \boldsymbol{u}_{t'}|$; for the remaining iterations we switch to a ratio-based update proportional to the ratio of the current objective to the current constraint violation, which produces larger and more decisive increases; during the last 25% of iterations the step-size coefficient is further increased by a factor of 5 to accelerate final convergence to feasibility. See Cory-Wright and Pauphilet (2026) for a full description and theoretical justification of the update rule.

Penalty weights. We write $\lambda_{t,t'} = \lambda \, w_{t'}$, with a single scalar $\lambda$ carrying the schedule above and a per-component weight $w_{t'}$ fixed at the first iteration. The weight is

$$w_{t} = \frac{\boldsymbol{u}_t^\top \boldsymbol{\Sigma} \boldsymbol{u}_t} {\|\boldsymbol{C} \boldsymbol{u}_t\|_2^2},$$

the same expression for both constraint types. Its effect is that the penalty a component can contribute is bounded by $\lambda \, w_{t'} \|\boldsymbol{C} \boldsymbol{u}_{t'}\|_2^2 = \lambda \, \boldsymbol{u}_{t'}^\top \boldsymbol{\Sigma} \boldsymbol{u}_{t'}$, i.e. $\lambda$ times the variance that component explains, whichever $\boldsymbol{C}$ is in force. The scalar $\lambda$ is thus a dimensionless penalty-to-objective ratio and the schedule behaves the same way under both constraints. Under orthogonality $\boldsymbol{C} = \mathbb{I}$ and $\|\boldsymbol{u}_t\|_2 = 1$, so $w_t$ is simply the variance explained by component $t$; under zero correlation $\boldsymbol{C} = \boldsymbol{\Sigma}/\mathrm{tr}(\boldsymbol{\Sigma})$ and the $1/\|\boldsymbol{C} \boldsymbol{u}_t\|_2^2$ factor offsets the shrinkage that normalizing $\boldsymbol{\Sigma}$ by its trace would otherwise apply to the penalty. Numerator and denominator are homogeneous of the same degree in $\boldsymbol{\Sigma}$, so the weights are invariant to a rescaling of the data. With this weight the PSD shift of the previous section also takes the common form $\lambda_0 = \lambda \sum_{t' \neq t} \boldsymbol{u}_{t'}^\top \boldsymbol{\Sigma} \boldsymbol{u}_{t'}$.

Termination. Algorithm 1 stops when the number of outer iterations reaches maxIter (default 200), or earlier at any iteration where the current solution is feasible and the change in objective value since the previous iteration falls below stallingTolerance (default 1e-8).

Feasibility tracking. At each iteration Algorithm 1 checks whether the current solution satisfies the coupling constraint up to feasibilityTolerance (default 1e-4). The best feasible solution encountered across all iterations is returned. If no feasible solution is found within the iteration budget, the algorithm returns the solution with the smallest observed constraint violation.

Software implementation. All computations are carried out in C++ via the Rcpp (Eddelbuettel and François 2011) and RcppEigen (Bates and Eddelbuettel 2013) interfaces, with a lightweight R wrapper providing the user-facing API. To avoid materializing the $p \times p$ perturbed matrix $\tilde{\boldsymbol{\Sigma}}_t$ at each inner-loop step, the C++ back-end represents it implicitly: each product $\tilde{\boldsymbol{\Sigma}}_t \boldsymbol{\beta}$ is evaluated as

$$\boldsymbol{\Sigma}\boldsymbol{\beta} - \boldsymbol{W}(\boldsymbol{d} \odot \boldsymbol{W}^\top \boldsymbol{\beta}) + \lambda_0 \boldsymbol{\beta},$$

where $\boldsymbol{W}$ collects the previously computed components ($\boldsymbol{u}_{t'}$, $t' \neq t$ for orthogonal loadings; $\boldsymbol{C} \boldsymbol{u}_{t'} = \boldsymbol{\Sigma} \boldsymbol{u}_{t'} / \mathrm{tr}(\boldsymbol{\Sigma})$, $t' \neq t$, for uncorrelated PCs) and $\boldsymbol{d}$ contains the corresponding scaled penalty coefficients. This eliminates the $O(r p^2)$ matrix-build cost per component update while keeping the per-step cost at $O(p^2)$.

When the raw data matrix $\boldsymbol{X}$ is provided (type = "X"), the product $\boldsymbol{\Sigma}\boldsymbol{\beta}$ is replaced by the two-pass evaluation $\boldsymbol{X}^\top(\boldsymbol{X}\boldsymbol{\beta})/(n-1)$ at cost $O(np)$ instead of $O(p^2)$, which is substantially more scalable when $n \ll p$ and avoids forming the $p \times p$ covariance matrix entirely. After the first outer iteration the previous iterate serves as a warm start for Algorithm 2, substantially reducing the number of random restarts required.

Computational complexity

The dominant cost of Algorithm 1 per outer iteration is $r$ calls to Algorithm 2. With the implicit matrix-vector representation above, applying $\tilde{\boldsymbol{\Sigma}}_t$ to a vector costs $O(p^2 + rp)$ with type = "Sigma" and $O(np + rp)$ with type = "X". In both cases the $O(rp)$ deflation term is negligible for moderate $r$. Each call performs at most $L_{\mathrm{TPM}}$ such products, giving a worst-case per-outer-iteration cost of $O(r \min(n,p) p \cdot L_{\mathrm{TPM}})$. In practice, warm-start initialization and early convergence detection reduce the effective number of TPM iterations substantially, so the empirical cost is much closer to $O(r \min(n,p) p)$ per outer iteration.

Guidance on parameter choices

Choosing the sparsity budgets ks

The budgets $k_1, \ldots, k_r$ are the primary tuning parameters. A practical approach is to run mspca() over a range of values and plot the trade-off between FVE and sparsity:

library("msPCA")
Sigma <- cor(datasets::mtcars)

ks_grid <- seq(2, 10, by = 1)
trade_off <- sapply(ks_grid, function(k) {
  set.seed(42)
  res <- mspca(Sigma, r = 3, ks = rep(k, 3), verbose = FALSE)
  fraction_variance_explained(Sigma, res$x_best)
})
plot(ks_grid, trade_off, type = "b",
     xlab = "sparsity budget k", ylab = "fraction of variance explained")

Domain knowledge often provides a natural guide: if each PC is expected to represent a distinct thematic cluster of features, setting $k_t$ to the anticipated cluster size is a good starting point.

Choosing the constraint type

Orthogonality (feasibilityConstraintType = 0) is appropriate when the loading vectors are to be used as a projection basis, or when the geometric structure of the components matters. Zero pairwise correlation (feasibilityConstraintType = 1) is preferable when the primary goal is statistical decorrelation of the projected data. In our experience the two options yield similar results when $\boldsymbol{\Sigma}$ is close to the identity, but can differ noticeably for strongly correlated datasets. See vignette("case-study-snp500", package = "msPCA") for a worked comparison.

Both sets of pairwise violations are computed at fit time and stored in nonredundancy, so a solution can be scored under the definition that was not enforced without a refit:

res <- mspca(Sigma, r = 3, ks = rep(5, 3), feasibilityConstraintType = 1, verbose = FALSE)
summary(res)                            # zero-correlation violations, as fitted
res$nonredundancy$orthogonality         # how far the same solution is from orthogonal

Iteration and restart budgets

maxIter (default 200) caps the number of outer iterations. Lowering it speeds up large problems at some risk of returning a less-refined solution; the case study uses maxIter = 100 on a 423-variable problem without noticeable loss. maxRestartTPM and minRestartTPM control the number of random restarts in the inner TPM call at the first and subsequent outer iterations respectively; the defaults (30 and 20) are conservative and can be reduced when runtime matters more than guarding against poor local optima.

References

Bates, Douglas, and Dirk Eddelbuettel. 2013. “Fast and Elegant Numerical Linear Algebra Using the RcppEigen Package.” Journal of Statistical Software 52 (5): 1–24.

Behdin, Kayhan, and Rahul Mazumder. 2026. “Sparse PCA: A New Scalable Estimator Based on Integer Programming.” Annals of Statistics.

Benidis, Konstantinos, Ying Sun, Prabhu Babu, and Daniel P. Palomar. 2016. “Orthogonal Sparse PCA and Covariance Estimation via Procrustes Reformulation.” IEEE Transactions on Signal Processing 64 (23): 6211–26.

Berk, Lauren, and Dimitris Bertsimas. 2019. “Certifiably Optimal Sparse Principal Component Analysis.” Mathematical Programming Computation 11 (3): 381–420.

Bertsimas, Dimitris, Ryan Cory-Wright, and Jean Pauphilet. 2022. “Solving Large-Scale Sparse PCA to Certifiable (Near) Optimality.” Journal of Machine Learning Research 23 (13): 1–35.

Cory-Wright, Ryan, and Jean Pauphilet. 2026. “Sparse PCA with Multiple Components.” Operations Research. https://doi.org/10.1287/opre.2023.0598.

d’Aspremont, Alexandre, Laurent El Ghaoui, Michael I. Jordan, and Gert R. G. Lanckriet. 2007. “A Direct Formulation for Sparse PCA Using Semidefinite Programming.” SIAM Review 49 (3): 434–48.

Eddelbuettel, Dirk, and Romain François. 2011. “Rcpp: Seamless R and C++ Integration.” Journal of Statistical Software 40 (8): 1–18.

Journée, Michel, Yurii Nesterov, Peter Richtárik, and Rodolphe Sepulchre. 2010. “Generalized Power Method for Sparse Principal Component Analysis.” Journal of Machine Learning Research 11 (2): 517–53.

Lu, Zhaosong, and Yong Zhang. 2012. “An Augmented Lagrangian Approach for Sparse Principal Component Analysis.” Mathematical Programming 135: 149–93.

Mackey, Lester. 2008. “Deflation Methods for Sparse PCA.” Advances in Neural Information Processing Systems 21.

Vu, Vincent Q., Juhee Cho, Jing Lei, and Karl Rohe. 2013. “Fantope Projection and Selection: A Near-Optimal Convex Relaxation of Sparse PCA.” Advances in Neural Information Processing Systems 26.

Yuan, Xiao-Tong, and Tong Zhang. 2013. “Truncated Power Method for Sparse Eigenvalue Problems.” Journal of Machine Learning Research 14 (4): 899–925.

Zou, Hui, Trevor Hastie, and Robert Tibshirani. 2006. “Sparse Principal Component Analysis.” Journal of Computational and Graphical Statistics 15 (2): 265–86.