Benchmarking against other sparse PCA packages • msPCA

This vignette gives the benchmarking protocol and results behind the summary in the package description. It compares msPCA against eight competing implementations, drawn from seven packages, on four datasets, and records how each competitor was configured. That configuration matters more than usual here, because the packages parameterize sparsity in mutually incompatible ways.

Reproducing the comparison requires eight additional R packages, plus callr and a Python installation, and the memory measurements require macOS or Linux. The scripts are in the package repository under replication/benchmarking/.

Getting the data

All four datasets are obtainable from R without manual downloads or accounts, so the comparison is reproducible end to end:

## mtcars -- base R
X_mtcars     <- as.matrix(datasets::mtcars)
Sigma_mtcars <- cor(X_mtcars)

## Pitprops -- the 13 x 13 correlation matrix of Jeffers (1967), Table 1.
## It is transcribed literally at the top of
## replication/benchmarking/notebook_pitprops.R; source that file to obtain `pitprops`.
Sigma_pitprops <- pitprops
## Methods that require a data matrix instead use pseudo-data drawn from this
## correlation structure (see "How each method was configured" below).
set.seed(42)
X_pitprops <- MASS::mvrnorm(n = 500, mu = rep(0, ncol(pitprops)),
                            Sigma = pitprops)

## Breast cancer -- one-off download via PMA, then keep the 500 highest-
## variance genes, following the preprocessing of Witten et al. (2009).
X            <- t(PMA::download_breast_data(
  url = "https://tibshirani.su.domains/PMA/breastdata.rda")$rna)
X_breast     <- X[, order(apply(X, 2, var), decreasing = TRUE)[1:500]]
Sigma_breast <- cor(X_breast)          # n = 89, p = 500

## Riboflavin -- ships with the hdi package
data("riboflavin", package = "hdi")
X_ribo <- as.matrix(riboflavin$x)      # n = 71, p = 4088

Methods compared

We compare msPCA::mspca() against eight functions drawn from seven competing packages. Each entry below names the exact function called, not merely the package, since some packages expose more than one sparse PCA routine.

Function	Reference	Language	Approach
`elasticnet::spca()`	Zou et al. (2006)	R	elastic-net regularization, LARS-like
`PMA::SPC()`	Witten et al. (2009)	R	penalized matrix decomposition
`sparsepca::spca()`	Erichson et al. (2020)	R	variable projection with $\ell_1$ penalty
`mixOmics::spca()`	Rohart et al. (2017)	R	regularized SVD (Shen and Huang 2008)
`nsprcomp::nsprcomp()`	Sigg (2019)	R	expectation-maximization (Sigg and Buhmann 2008) + deflation
`nsprcomp::nscumcomp()`	Sigg (2019)	R	cumulative variant, all components jointly
`amanpg::spca.amanpg()`	Chen et al. (2020)	R	alternating manifold proximal gradient
`sklearn.decomposition.SparsePCA()`	Pedregosa et al. (2011)	Python	dictionary learning (Mairal et al. 2009)

scikit-learn is called from R via the reticulate package (Ushey et al. 2025). Together these cover the main methodological approaches to sparse PCA: elastic-net regularization, penalized matrix decomposition, variable-projection optimization, regularized SVD, deflation, Riemannian proximal gradient, and dictionary learning.

How each function obtains multiple components, and what its tuning parameter controls

The comparison turns on two things that are easy to conflate: how a method gets from one sparse component to $r$ of them, and what the parameter the user actually turns is doing. The second column below states the mechanism; the last two name the parameter we tune and classify it as controlling sparsity or as an orthogonality penalty.

Function	Multiple-component mechanism	Parameter tuned	What it controls
`msPCA::mspca()`	Iterative deflation heuristic; all $r$ components are refined together, with the non-redundancy constraint enforced across every pair and its penalty tightened over the iterations	`ks`	Sparsity — exact cardinality, one budget per component
`elasticnet::spca()`	All $r$ jointly: alternates an orthonormal Procrustes update of an auxiliary factor with per-component elastic-net regressions. Orthonormality is imposed on the auxiliary factor, not on the returned sparse loadings	`para` (with `sparse = "varnum"`)	Sparsity — exact cardinality per component
`PMA::SPC()`	Sequential rank-1 penalized matrix decompositions with deflation. `orth = TRUE` selects the variant of Section 3.2 of Witten et al. (2009), which requires the left factors $\boldsymbol{u}_k$ to be mutually orthogonal — not the sparse loadings $\boldsymbol{v}_k$ that the function returns and that we score	`sumabsv`	Sparsity — an $\ell_1$ bound on each loading vector, not a cardinality
`sparsepca::spca()`	All $r$ jointly, by variable projection on the SPCA objective	`alpha`	Sparsity — magnitude of the $\ell_1$ penalty, with no per-component control
`mixOmics::spca()`	Sequential regularized rank-1 SVDs with deflation	`keepX`	Sparsity — exact number of variables retained per component
`nsprcomp::nsprcomp()`	Sequential cardinality-constrained components via an EM-style iteration, combined with the deflation of Mackey (2008)	`k`	Sparsity — exact cardinality per component
`nsprcomp::nscumcomp()`	All $r$ jointly, under one cumulative budget	`k`, `gamma`	`k`: sparsity, as a total budget $\sum_t k_t$ ; `gamma`: orthogonality penalty, which we tune (see below) rather than leaving at its default of 0
`amanpg::spca.amanpg()`	All $r$ jointly: alternates a proximal-gradient step on the sparse loadings with a manifold update of an orthonormal factor on the Stiefel manifold	`lambda1`	Sparsity — a vector of per-component $\ell_1$ penalties
`sklearn...SparsePCA()`	All `n_components` jointly, by dictionary learning	`alpha`	Sparsity — magnitude of the $\ell_1$ penalty on the loadings

Compared with msPCA::mspca():

Among the competitors, only nsprcomp::nscumcomp() offers any control over non-redundancy (orthogonality), through the soft penalty gamma, whose default of 0 applies no penalty at all. For a fair comparison, we tune gamma along a grid: among the values that bring the orthogonality violation within 1e-4 (our default tolerance in msPCA), we keep the one with the highest FVE.

Not all methods control sparsity in the same way. Four functions — msPCA::mspca(), elasticnet::spca(), mixOmics::spca() and nsprcomp::nsprcomp() — accept an exact cardinality per component. The others expose either a penalty parameter that controls sparsity only indirectly (sparsepca::spca(), amanpg::spca.amanpg(), sklearn...SparsePCA()), an $\ell_1$ bound (PMA::SPC()), or a single cumulative sparsity budget covering all components together (nsprcomp::nscumcomp()). For these methods, we tune the hyperparameter so as to (i) return the requested number of non-empty components, $r$ , and (ii) come as close as possible to the requested total sparsity, $\sum_t k_t$ .

Environment

Experiments were conducted on a MacBook Air with an Apple M2 chip and 24 GB of RAM, with R version 4.6.0 (2026-04-24), msPCA v0.5.1, elasticnet v1.3, PMA v1.2.4, sparsepca v0.1.2, mixOmics v6.36.0, nsprcomp v0.5.1.2, amanpg v0.3.4, reticulate v1.46.0, callr v3.8.0, and scikit-learn v1.6.1 (Python 3.11). Full session information is recorded in replication/sessionInfo.txt.

Configuration

The table above gives each function’s mechanism and tuning parameter. What remains to be fixed is which input each one receives, since some accept a covariance or correlation matrix, some require the data matrix, and some take either.

msPCA::mspca() — run with both type = "Sigma" and type = "X" on all four datasets. On riboflavin ( $n \ll p$ ) the $p \times p$ covariance matrix is rank-deficient, making type = "X" the natural choice; we report both modes to quantify the runtime and memory gain.
elasticnet::spca() — takes a list of cardinality budgets (one per PC) via sparse = "varnum", and applies interchangeably to $\boldsymbol{\Sigma}$ (type = "Gram") or $\boldsymbol{X}$ (type = "predictor"). We report both input modes on mtcars, Pitprops and breast cancer.
PMA::SPC() — applied to the data matrix on mtcars, Pitprops and breast cancer, as documented (x is “a data matrix of dimension $n \times p$ ”, and SPC() centers its columns itself). On Pitprops that means the pseudo-data, since no real data matrix exists. sumabsv bounds the $\ell_1$ norm of each loading vector (subject to unit $\ell_2$ norm) and must lie between 1 and $\sqrt{p}$ ; we tune it per dataset to bring per-component sparsity as close to the target as the parameterization allows. We set orth = TRUE, which selects the multiple-component variant of Section 3.2 of Witten et al. (2009).
sparsepca::spca() — applied to $\boldsymbol{\Sigma}$ . It exposes only the magnitude of the $\ell_1$ penalty alpha, which we tune per dataset to achieve a total sparsity $\sum_t k_t$ close to the target. This method tends to return PCs with very unequal sparsity levels, the first typically denser than the rest.
mixOmics::spca() — takes the data matrix and the list of per-component budgets keepX, so it hits the cardinality target exactly.
nsprcomp::nsprcomp() and nsprcomp::nscumcomp() — two functions from the same package, both based on the EM algorithm of Sigg and Buhmann (2008). nsprcomp() computes the leading sparse PC then applies the deflation procedure of Mackey (2008) for subsequent PCs; it takes the data matrix and a vector of per-component cardinalities. nscumcomp() solves for all $r$ PCs simultaneously under a single cumulative budget $\sum_t k_t$ , which it distributes as it sees fit, plus the penalty gamma on divergence from orthonormality.
amanpg::spca.amanpg() — accepts either input and is told which via type: type = 0 declares z to be a data matrix, any non-zero value declares it a covariance matrix. We use type = 1 with the correlation matrix on mtcars and Pitprops, where the correlation matrix is the natural object, and type = 0 with the data matrix on breast cancer. It computes all $r$ components jointly, alternating a proximal-gradient step on the sparse loadings with a manifold update of an orthonormal factor on the Stiefel manifold; lambda1 is a vector of per-component $\ell_1$ penalties, which we tune per dataset, and we set lambda2 = Inf.
sklearn.decomposition.SparsePCA() — called via reticulate. The class computes n_components components simultaneously using dictionary learning with an $\ell_1$ penalty alpha on the loadings, applied to the data matrix. We tune alpha per dataset and fix random_state to the repetition’s seed.

For Pitprops only the correlation matrix is available, so we generate pseudo-data via MASS::mvrnorm() with the Pitprops matrix as the population covariance ( $n = 500$ ) for the methods that require $\boldsymbol{X}$ . FVE and the orthogonality violation are always evaluated against the published Pitprops correlation matrix, never against the covariance of the pseudo-data.

Every tuning parameter is selected by the scripts rather than fixed by hand: each is swept over a grid and the value whose realized sparsity comes closest to the target is kept, with ties broken on FVE and candidates that leave a component empty discarded. The selected values appear beneath each table. Tuning is excluded from all reported runtimes and memory figures.

How runtime and memory are measured

Every repetition of every method runs in its own fresh R process, driven by replication/benchmarking/bench_utils.R. Isolation is necessary because peak memory is read as a process-level high-water mark: several methods in one session, or even several repetitions of one method, contaminate each other.

Memory is obtained from getrusage(RUSAGE_SELF).ru_maxrss. R-level metrics such as gc() or bench::mark(mem_alloc) are deliberately not used: msPCA does its work in RcppEigen buffers and scikit-learn in the embedded CPython heap, neither of which R’s garbage collector observes, so an R-level metric would under-report memory usage for those two methods. The Memory column below reports the working set, defined as the size of the matrix handed to the method plus the additional peak RSS the method allocates on top of it.

Each method is repeated five times with a different seed. All four measured quantities are reported as medians. For FVE, runtime, and memory, we also report a bracketed min–max range wherever the spread is visible at the printed precision; for the orthogonality violation, we report the worst (i.e., largest) value.

# Representative call pattern; see replication/benchmarking/ for the full scripts.
set.seed(42)
res <- mspca(Sigma, r = 3, ks = rep(4, 3), verbose = FALSE)

fraction_variance_explained(Sigma, res$x_best)
feasibility_violation_off(Sigma, res$x_best, 0)
res$runtime

Results

mtcars

datasets::mtcars from base R ( $p = 11$ variables, $n = 32$ observations) serves as a simple reproducible baseline. We extract $r = 3$ sparse PCs targeting $k = 4$ per component.

Method	$(k_1,k_2,k_3)$	FVE	Orth. violation	Runtime (s)	Memory (MB)
`msPCA::mspca` $(\boldsymbol{\Sigma})$	(4,4,4)	0.829	9.4e-05	0.189	0.4 [0.4, 0.7]
`msPCA::mspca` $(\boldsymbol{X})$	(4,4,4)	0.829	9.4e-05	0.245	0.2
`elasticnet::spca` $(\boldsymbol{\Sigma})$	(4,4,4)	0.753	0.117	0.093	21
`elasticnet::spca` $(\boldsymbol{X})$	(4,4,4)	0.754	0.117	0.112	23 [23, 28]
`PMA::SPC`	(4,4,5)	0.551	0.221	0.009	0.6 [0.6, 3.6]
`sparsepca::spca`	(6,4,1)	0.581	0.133	0.032 [0.031, 0.036]	17
`mixOmics::spca`	(4,4,4)	0.633	0.133	0.003	0.5
`nsprcomp::nsprcomp`	(4,4,4)	0.832	0.078	0.004 [0.003, 0.005]	0.1
`nsprcomp::nscumcomp`	(5,5,2)	0.824 [0.798, 0.825]	1.1e-05 (max 3.3e-05)	0.241	30
`amanpg::spca.amanpg`	(4,4,1)	0.544	0.047	0.007	1.5
`sklearn...SparsePCA`	(4,3,4)	0.827	< 1e-15	0.012 [0.010, 0.015]	0.0 [0.0, 0.1]
PCA (dense)	(11,11,11)	0.899	< 1e-15	0.001	0.5

Tuning: sumabsv = 1.85, alpha = 0.0045 (sparsepca), lambda1 = c(4.5, 5.4, 0.1), alpha = 3.2 (sklearn), gamma = 1e6.

nscumcomp does well here: it is within 0.01 FVE of the best sparse method and comfortably inside the feasibility tolerance. However, it returns a sparsity pattern of (5,5,2) rather than the requested (4,4,4), and it is the only method whose FVE varies across seeds, spanning 0.798–0.825.

Because $p = 11$ , every method’s numerical working set is a few kilobytes, so the memory column measures fixed overheads: most methods sit at or below 1.5 MB, while elasticnet, sparsepca and nscumcomp carry a floor of 17–30 MB that is already present at this size.

Pitprops

The Pitprops dataset (Jeffers 1967) ( $p = 13$ physical measurements on $n = 180$ timber specimens) is a classic sparse PCA benchmark (e.g., Zou et al. 2006). We extract $r = 6$ sparse PCs targeting $k = 4$ nonzero loadings per component, as in Zou et al. (2006).

Method	$(k_1, \ldots, k_6)$	FVE	Orth. violation	Runtime (s)	Memory (MB)
`msPCA::mspca` $(\boldsymbol{\Sigma})$	(4,4,4,4,4,4)	0.791 [0.786, 0.822]	9.7e-05 (max 9.8e-05)	0.647	0.7 [0.4, 0.8]
`msPCA::mspca` $(\boldsymbol{X})$	(4,4,4,4,4,4)	0.794 [0.784, 0.801]	9.4e-05 (max 9.9e-05)	3.958	0.3 [0.3, 0.9]
`elasticnet::spca` $(\boldsymbol{\Sigma})$	(4,4,4,4,4,4)	0.792	0.241	0.070	21
`elasticnet::spca` $(\boldsymbol{X})$	(4,4,4,4,4,4)	0.790	0.322	0.158	73 [70, 83]
`PMA::SPC`	(3,4,4,4,3,4)	0.714	0.491	0.021	15
`sparsepca::spca`	(7,3,4,2,4,4)	0.756	0.353	0.039	28
`mixOmics::spca`	(4,4,4,4,4,4)	0.788	0.215	0.010 [0.009, 0.011]	0.6 [0.5, 0.6]
`nsprcomp::nsprcomp`	(4,4,4,4,4,4)	0.816 [0.813, 0.833]	0.563 (max 1.105)	0.025 [0.023, 0.028]	34 [28, 42]
`nsprcomp::nscumcomp`	(5,6,5,2,4,2)	0.819 [0.814, 0.830]	5.1e-05 (max 1.0e-04)	1.098	68 [64, 129]
`amanpg::spca.amanpg`	(3,3,4,3,4,6)	0.863	1.654	0.011 [0.011, 0.013]	4.5
`sklearn...SparsePCA`	(7,3,4,3,5,4)	0.907	0.815	0.061	0.1
PCA (dense)	(13, …, 13)	0.877	< 1e-15	0.001	0.4

Tuning: sumabsv = 1.55, alpha = 0.005 (sparsepca), lambda1 = c(4.80, 4.40, 1.95, 5.50, 0.20, 0.05), alpha = 2.6 (sklearn), gamma = 1e7.

This dataset is the one where the feasibility differences are starkest. With $p = 13$ and 24 nonzero coordinates in total, disjoint-support solutions do not exist, so methods that do not explicitly enforce orthogonality — the deflation-based ones in particular — struggle: at the requested cardinality elasticnet, mixOmics and nsprcomp return violations of 0.215–0.563 against 9.7e-05 for msPCA. nscumcomp also reaches the tolerance (with $\gamma = 10^7$ against a default of 0), but at a fortyfold slowdown relative to nsprcomp (1.098 s against 0.025 s) and without achieving the requested sparsity.

Pitprops is also the one case where type = "X" is the wrong choice: with $n = 500$ pseudo-observations and $p = 13$ the data matrix is larger than the correlation matrix, and runtime is roughly six times worse. Across seeds msPCA ranges over 0.786–0.822 in FVE, the widest spread it shows on any of the four datasets; the orthogonality violation stays below 1e-04 throughout.

Breast cancer gene expression

Gene expression data from Chin et al. (2006) ( $n = 89$ tumor samples, 19,672 genes). Following the preprocessing of Witten et al. (2009) we retain the $p = 500$ genes with the highest marginal variance, yielding a high-dimensional ( $p > n$ ) benchmark. We set $r = 3$ and target $k = 20$ nonzero loadings per component.

Method	$(k_1, k_2, k_3)$	FVE	Orth. violation	Runtime (s)	Memory (MB)
`msPCA::mspca` $(\boldsymbol{\Sigma})$	(20,20,20)	0.093 [0.088, 0.095]	7.5e-05 (max 9.8e-05)	27.244 [14.119, 27.602]	9.4 [6.6, 11.3]
`msPCA::mspca` $(\boldsymbol{X})$	(20,20,20)	0.095 [0.088, 0.095]	7.5e-05 (max 9.8e-05)	13.268 [6.790, 13.399]	3.2 [2.8, 4.4]
`elasticnet::spca` $(\boldsymbol{\Sigma})$	(20,20,20)	0.043	0.040	16.104	192 [149, 220]
`elasticnet::spca` $(\boldsymbol{X})$	(20,20,20)	0.043	0.021	0.816	79 [79, 127]
`PMA::SPC`	(20,22,13)	0.042	0.008	0.041	33
`sparsepca::spca`	(43,14,22)	0.082	< 1e-15	4.841	142 [46, 176]
`mixOmics::spca`	(20,20,20)	0.077	< 1e-15	0.032	7.0 [6.9, 27.6]
`nsprcomp::nsprcomp`	(20,20,20)	0.083 [0.070, 0.084]	< 1e-15	0.032 [0.030, 0.035]	48
`nsprcomp::nscumcomp`	(29,23,8)	0.081 [0.080, 0.085]	< 1e-15	0.229 [0.219, 0.261]	55
`amanpg::spca.amanpg`	(68,1,2)	0.045	< 1e-15	0.071	50 [49, 65]
`sklearn...SparsePCA`	(22,3,1)	0.020 [0.020, 0.032]	< 1e-15	0.079	0.6 [0.5, 1.3]
PCA (dense)	(500,500,500)	0.352	< 1e-15	0.010	2.6 [2.5, 3.6]

Tuning: sumabsv = 3.3, alpha = 0.0035 (sparsepca), lambda1 = c(4.5, 1.8, 3.3), alpha = 15 (sklearn), gamma = 316.

This is where the penalty-parameterized methods break down most clearly in terms of achieved sparsity. Against a request of (20,20,20), amanpg returns (68,1,2), sklearn (22,3,1) and sparsepca (43,14,22); no value on their grids does better. Only the four functions taking a per-component cardinality budget return the desired sparsity levels. Feasibility, by contrast, is uninformative here: with 60 nonzeros over $p = 500$ disjoint supports are easy to find, and most methods return solutions that are trivially orthogonal. msPCA is the exception, at 7.5e-05: it stops once the violation is inside its 1e-4 tolerance rather than driving it to zero.

This is the smallest dataset at which the memory column measures the algorithms rather than fixed overheads, and msPCA on the data matrix is the lightest of the methods that meet the budget: 3.2 MB, against 7.0 MB for mixOmics, 9.4 MB for msPCA on the correlation matrix, 48 MB for nsprcomp, and 79 and 192 MB for elasticnet. (sparsepca at 142 MB and nscumcomp at 55 MB do not meet the budget.) The type = "X" saving is modest in absolute terms because $\boldsymbol{\Sigma}$ is only about 2 MB at this size.

Riboflavin

The riboflavin (vitamin B2) production dataset (Bühlmann et al. 2014) comprises log-transformed expression levels of $p = 4{,}088$ genes in $n = 71$ samples of Bacillus subtilis. This is a challenging $p \gg n$ benchmark: the empirical correlation matrix is rank-deficient (rank $\leq 70$ ). We set $r = 2$ and target $k = 20$ nonzero loadings per component. To avoid time-consuming parameter tuning, we consider only methods that accept exact cardinality budgets.

Method	$(k_1, k_2)$	FVE	Orth. violation	Runtime (s)	Memory (MB)
`msPCA::mspca` $(\boldsymbol{\Sigma})$	(20,20)	0.008 [0.008, 0.009]	< 1e-15	57.226 [41.791, 87.700]	640
`msPCA::mspca` $(\boldsymbol{X})$	(20,20)	0.008 [0.008, 0.009]	< 1e-15	1.557 [0.777, 2.177]	44
`elasticnet::spca` $(\boldsymbol{X})$	(20,20)	0.004	< 1e-15	6.303	512 [421, 577]
`mixOmics::spca`	(20,20)	0.006	< 1e-15	0.335	24 [24, 50]
`nsprcomp::nsprcomp`	(20,20)	0.009	< 1e-15	0.370 [0.265, 0.391]	61 [57, 124]
PCA (dense)	(4088,4088)	0.478	< 1e-15	0.046	39

msPCA is run with maxIter = 100 here, and every method other than msPCA $(\boldsymbol{\Sigma})$ takes the data matrix.

This dataset isolates the benefit of type = "X". Passing the raw data matrix replaces the $O(p^2)$ product $\boldsymbol{\Sigma}\boldsymbol{\beta}$ with the $O(np)$ two-pass product $\boldsymbol{X}^\top(\boldsymbol{X}\boldsymbol{\beta})/(n-1)$ and never forms the covariance matrix: the working set falls from 640 MB to 44 MB and median runtime from 57.2 seconds to 1.6 seconds, at identical FVE and feasibility. The correlation matrix accounts for 128 MB of the difference, the rest being the Eigen copy of it and the dense intermediates the $O(p^2)$ path requires. Both rows vary by a factor of two or more in runtime across seeds, so those ratios are good to one significant figure at most.

The memory figures are far steadier: msPCA $(\boldsymbol{\Sigma})$ spans 639–641 MB over the same repetitions, whereas elasticnet’s five repetitions span 421–577 MB and its median moves between runs of the identical script. msPCA allocates in Eigen, where the footprint follows from the algorithm; the pure-R methods allocate on R’s heap, where the peak depends on garbage-collector timing. Together with the Pitprops result the rule is just the relative size of the two inputs: type = "X" when $n \ll p$ , type = "Sigma" when $p \lesssim n$ .

Discussion

Sparsity control. PMA::SPC(), sparsepca::spca(), nsprcomp::nscumcomp(), amanpg::spca.amanpg() and sklearn.decomposition.SparsePCA() give no direct control over the sparsity of each loading vector. PMA::SPC() exposes only a global $\ell_1$ bound and sklearn only a global $\ell_1$ penalty; sparsepca::spca(), nscumcomp() and amanpg tend to return unbalanced components, with one PC carrying most of the budget. None can be calibrated to an exact per-component target, which is what makes an apples-to-apples comparison awkward. Where the discussion below needs methods that meet the budget exactly, those are msPCA::mspca(), elasticnet::spca(), mixOmics::spca() and nsprcomp::nsprcomp().

Feasibility. Two things have to be separated here, because the answer depends on whether disjoint supports exist at all. On breast cancer and riboflavin, where $p$ is large relative to the total budget, orthogonality is nearly free: most methods return exactly orthogonal loadings simply because their supports do not overlap. On mtcars ( $p = 11$ , 12 nonzeros) and Pitprops ( $p = 13$ , 24 nonzeros) no disjoint solution exists, and that is where the differences appear.

There, msPCA::mspca() reaches 1e-4 while returning the requested cardinality exactly, and does so at every repetition: the largest violation it produced anywhere in the benchmark is 9.9e-05. Two other methods reach the tolerance, both with qualifications. nsprcomp::nscumcomp() reaches it on both datasets, but only with gamma driven to $10^6$ – $10^7$ , returning (5,5,2) and (5,6,5,2,4,2) against budgets of (4,4,4) and (4,4,4,4,4,4), and its worst Pitprops repetition is 1.03e-04. sklearn...SparsePCA() reaches it on mtcars by returning 11 nonzeros in total rather than the requested 12, which at $p = 11$ admits a disjoint solution and so makes orthogonality free. Every other method, at the requested cardinality, returns violations between 0.08 and 0.56 on these two datasets.

The claim supported by this is narrower than “only msPCA can be feasible”, and more useful: on the two datasets where feasibility actually binds, it is the only method that delivers feasibility and the requested sparsity together, without a search over a penalty parameter to get there.

Variance explained. Among the methods that meet the requested cardinality, msPCA::mspca() and nsprcomp::nsprcomp() are generally the best performers, followed by mixOmics::spca() and elasticnet::spca(). The FVE cost of enforcing feasibility therefore looks small: the largest gap between msPCA and the best budget-respecting competitor is about 0.025 FVE, on Pitprops, where nsprcomp reaches 0.816 at an orthogonality violation of 0.563.

Several methods post higher FVE while missing the budget, and those figures should be read in that light: sklearn...SparsePCA() reaches 0.907 on Pitprops at (7,3,4,3,5,4) with a violation of 0.815, and amanpg::spca.amanpg() 0.863 at a violation of 1.654. The exception is nscumcomp at 0.824 on mtcars with (5,5,2), which is both high and feasible, and is the most substantive competition in the comparison.

Reproducibility. Repeating each method under several seeds separates the implementations that returned identical results at every seed (elasticnet, PMA, sparsepca, mixOmics, amanpg, dense PCA) from the rest. nsprcomp and nscumcomp use random initialization and move between runs; nsprcomp’s orthogonality violation on Pitprops spans 0.56–1.10. sklearn...SparsePCA(), seeded from the repetition’s seed, also moves (0.020–0.032 FVE on breast cancer). msPCA uses random restarts and is likewise not bit-reproducible across seeds, but its variation is confined to the objective: the feasibility target was met at every seed on every dataset, the largest violation observed anywhere being 9.9e-05.

Memory. On the two high-dimensional datasets, where the memory column reflects the algorithms rather than fixed overheads, msPCA is among the lightest of the sparse methods that respect an exact cardinality budget. On breast cancer it uses 3.2 MB on the data matrix against 79 MB for elasticnet on the same input and 192 MB in Gram mode, 7.0 MB for mixOmics and 48 MB for nsprcomp (with sparsepca at 142 MB and nscumcomp at 55 MB, neither meeting the budget). On riboflavin it uses 44 MB against 512 MB for elasticnet and 61 MB for nsprcomp, though mixOmics is lighter still at 24 MB. Within msPCA, the $\boldsymbol{\Sigma}$ -versus- $\boldsymbol{X}$ choice matters more than the choice of package: the same function spans 640 MB to 44 MB on riboflavin depending on the interface.

Scalability. msPCA::mspca() scales to datasets with $p$ in the thousands in reasonable time, and the type = "X" interface yields substantial runtime and memory gains when $n \ll p$ : roughly 35x faster and 15x lighter on riboflavin. When $n > p$ , as in the Pitprops pseudo-data, the ordering reverses and type = "Sigma" is preferable.

References

Bühlmann, Peter, Markus Kalisch, and Lukas Meier. 2014. “High-Dimensional Statistics with a View Toward Applications in Biology.” Annual Review of Statistics and Its Application 1 (1): 255–78.

Chen, Shixiang, Shiqian Ma, Anthony Man-Cho So, and Tong Zhang. 2020. “Proximal Gradient Method for Nonsmooth Optimization over the Stiefel Manifold.” SIAM Journal on Optimization 30 (1): 210–39.

Chin, Koei, Sandy DeVries, Jane Fridlyand, et al. 2006. “Genomic and Transcriptional Aberrations Linked to Breast Cancer Pathophysiologies.” Cancer Cell 10 (6): 529–41.

Erichson, N. Benjamin, Peng Zheng, Krithika Manohar, Steven L. Brunton, J. Nathan Kutz, and Aleksandr Y. Aravkin. 2020. “Sparse Principal Component Analysis via Variable Projection.” SIAM Journal on Applied Mathematics 80 (2): 977–1002.

Jeffers, J. N. R. 1967. “Two Case Studies in the Application of Principal Component Analysis.” Applied Statistics 16 (3): 225–36.

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

Mairal, Julien, Francis Bach, Jean Ponce, and Guillermo Sapiro. 2009. “Online dictionary learning for sparse coding.” Proceedings of the 26th International Conference on Machine Learning, 689–96.

Pedregosa, Fabian, Gaël Varoquaux, Alexandre Gramfort, et al. 2011. “Scikit-learn: Machine Learning in Python.” Journal of Machine Learning Research 12: 2825–30.

Rohart, Florian, Benoît Gautier, Amrit Singh, and Kim-Anh Lê Cao. 2017. “mixOmics: An R Package for ’Omics Feature Selection and Multiple Data Integration.” PLOS Computational Biology 13 (11): e1005752.

Shen, Haipeng, and Jianhua Z. Huang. 2008. “Sparse Principal Component Analysis via Regularized Low Rank Matrix Approximation.” Journal of Multivariate Analysis 99 (6): 1015–34.

Sigg, Christian D. 2019. nsprcomp: Non-Negative and Sparse PCA. https://CRAN.R-project.org/package=nsprcomp.

Sigg, Christian D., and Joachim M. Buhmann. 2008. “Expectation-Maximization for Sparse and Non-Negative PCA.” Proceedings of the 25th International Conference on Machine Learning, 960–67.

Ushey, Kevin, J. J. Allaire, and Yuan Tang. 2025. reticulate: Interface to Python. https://CRAN.R-project.org/package=reticulate.

Witten, Daniela M., Robert Tibshirani, and Trevor Hastie. 2009. “A Penalized Matrix Decomposition, with Applications to Sparse Principal Components and Canonical Correlation Analysis.” Biostatistics 10 (3): 515–34.

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