Skip to contents

Multiple Sparse PCA

Source: R/main.R
mspca.Rd

Returns multiple sparse principal components of a dataset using an iterative deflation heuristic. As in the elasticnet package, the data is passed as a single argument M whose interpretation is set by type: "Sigma" (the default) treats M as a covariance/correlation matrix (p x p) and "X" treats M as a raw data matrix (n observations x p variables). With type = "X" the algorithm operates on the data directly via the products \(X^\top(X\beta)\) and never forms the p x p matrix, which is substantially more scalable when \(n \ll p\).

Usage

mspca(
  M,
  r,
  ks,
  type = c("Sigma", "X"),
  feasibilityConstraintType = 0,
  verbose = TRUE,
  maxIter = 200,
  feasibilityTolerance = 1e-04,
  stallingTolerance = 1e-08,
  timeLimitTPM = 20,
  maxRestartTPM = 30,
  minRestartTPM = 20,
  center = TRUE,
  scale = TRUE,
  divisor = c("n-1", "n"),
  checkPSD = TRUE,
  symTolerance = 1e-08,
  psdTolerance = 1e-08
)

Arguments

M

A matrix. The data, interpreted according to type: a covariance/ correlation matrix (p x p) when type = "Sigma", or a raw data matrix (n x p) when type = "X".

r

An integer. Number of principal components (PCs) to be computed. Must be a single whole number greater than or equal to 1.

ks

An integer vector. Target sparsity of each PC. Every element must be a whole number greater than or equal to 1 (a sparsity level of zero or less does not define a component and is rejected). length(ks) must not exceed r; if it is shorter than r, a warning is issued and the algorithm is run for length(ks) PCs.

type

(optional) Either "Sigma" (default; M is a covariance/correlation matrix) or "X" (M is a raw data matrix).

feasibilityConstraintType

(optional) An integer. Type of feasibility constraints to be enforced. 0: orthogonality constraints; 1: uncorrelatedness constraints. Must be exactly 0 or 1; fractional values are rejected rather than rounded. Default 0.

verbose

(optional) A Boolean. Controls console output. TRUE/FALSE or the numbers 0/1; any other value is an error. Default TRUE.

maxIter

(optional) An integer. Maximum number of iterations of the algorithm. Must be a whole number greater than or equal to 1. Default 200.

feasibilityTolerance

(optional) A float. Tolerance for constraint violation (orthogonality/uncorrelatedness, according to feasibilityConstraintType). Under uncorrelatedness the violation is normalized by the total variance tr(Sigma). Must be non-negative; Inf is allowed and accepts any solution as feasible. Default 1e-4.

stallingTolerance

(optional) A float. Controls the objective improvement below which the algorithm is considered to have stalled. Must be non-negative; Inf is allowed. Default 1e-8.

timeLimitTPM

(optional) An integer. Maximum time in seconds for the truncated power method (inner iteration). Must be a finite, non-negative whole number; Inf is not accepted, so use a large finite value to effectively disable the limit. Default 20.

maxRestartTPM

(optional) An integer. Number of random restarts of the truncated power method (inner iteration) for the first outer iteration. Must be a whole number greater than or equal to 0; zero means no random restarts. Default 30.

minRestartTPM

(optional) An integer. Number of random restarts of the truncated power method (inner iteration) for outer iterations >= 2. Must be a whole number greater than or equal to 0. Default 20.

center

(optional, type = "X") A Boolean. Center the columns of M before computing the covariance. TRUE/FALSE or the numbers 0/1. Default TRUE.

scale

(optional, type = "X") A Boolean. Scale the columns of M to unit variance, i.e. operate on the correlation matrix. TRUE/FALSE or the numbers 0/1. Default TRUE.

divisor

(optional, type = "X") Either "n-1" (default, sample covariance, matches cov/cor) or "n" (population covariance). Default "n-1".

checkPSD

(optional, type = "Sigma") A Boolean. Verify that M is positive semidefinite. TRUE/FALSE or the numbers 0/1; any other value is an error. Default TRUE.

symTolerance

(optional, type = "Sigma") A float. Tolerance for the symmetry check on M. Must be non-negative. Default 1e-8.

psdTolerance

(optional, type = "Sigma") A float. Tolerance (on the smallest eigenvalue) for the PSD check on M. Must be non-negative. Default 1e-8.

Value

An object of class "mspca" (a list) with fields: x_best (p x r matrix of sparse PC loadings), objective_value, feasibility_violation, runtime, variance_explained (per-PC explained variance), total_variance (trace of the covariance matrix), feasibilityConstraintType (the value used to fit, reused as the default for all diagnostics), and nonredundancy. With type = "X" it additionally records inputType, center, scale, divisor, nObs, and p. Use print() to display the sparse loadings and summary() for a full per-PC breakdown.

nonredundancy is a list of two r x r matrices, orthogonality (\(|u_t^\top u_s|\)) and uncorrelatedness (\(|u_t^\top \Sigma u_s| / \mathrm{tr}(\Sigma)\)), computed once at fit time from the sorted loadings. Only the strict upper triangle is populated; the diagonal and lower triangle are NA. Both are stored regardless of which constraint was enforced, so a summary under either definition is available without the covariance matrix and without refitting.

The uncorrelatedness terms are normalized by the total variance \(\mathrm{tr}(\Sigma)\). The loadings being unit-norm, \(|u_t^\top \Sigma u_s|\) scales linearly with \(\Sigma\), so normalization makes the uncorrelatedness terms scale invariant.

Note that feasibility_violation (returned by the solver) is the quantity compared against feasibilityTolerance during the fit.

Examples

# From a covariance/correlation matrix (the default type):
TestMat <- cor(mtcars)
res <- mspca(TestMat, r = 2, ks = c(4, 4), verbose = FALSE)
print(res)
#> 
#> msPCA solution: 2 sparse PCs
#> Pct. variance explained: 29.83393 29.16617 
#> Non-zero loadings per PC: 4 4 
#> 
#> Sparse PCs
#>            [,1]       [,2]
#> mpg   0.5480282  0.0000000
#> cyl   0.0000000 -0.5061954
#> disp -0.3738294  0.0000000
#> hp    0.0000000 -0.5095722
#> drat  0.4980638  0.0000000
#> wt   -0.5584346  0.0000000
#> qsec  0.0000000  0.4714732
#> vs    0.0000000  0.5116791
summary(res)
#> 
#> msPCA summary: 2 sparse PC(s)
#> Input type   : Sigma 
#> Runtime (s)  : 0.014 
#> Constraint   : orthogonality (as fitted) 
#> 
#> Per-component statistics:
#>   PC nonzero variance       fve cumulative_fve max_violation
#>  PC1       4 3.281733 0.2983393      0.2983393             0
#>  PC2       4 3.208278 0.2916617      0.5900010             0
#> 
#> Pairwise orthogonality violations (upper triangle):
#>     PC1 PC2
#> PC1   .   0
#> PC2   .   .
#> Total: 0e+00 
#> Violations under the other definition are stored in `$nonredundancy$uncorrelatedness`.
# Equivalent call from the raw data matrix:
res_X <- mspca(as.matrix(mtcars), r = 2, ks = c(4, 4), type = "X", verbose = FALSE)
print(res_X)
#> 
#> msPCA solution: 2 sparse PCs
#> Pct. variance explained: 32.45835 27.98031 
#> Non-zero loadings per PC: 4 4 
#> 
#> Sparse PCs
#>            [,1]       [,2]
#> mpg   0.4994876  0.0000000
#> cyl  -0.4952709  0.0000000
#> disp -0.5096594  0.0000000
#> hp    0.0000000  0.5180509
#> wt   -0.4954454  0.0000000
#> qsec  0.0000000 -0.5056645
#> vs    0.0000000 -0.4935973
#> carb  0.0000000  0.4819631