Skip to contents

Worked example: msPCA on mtcars

Source: vignettes/msPCA.Rmd
msPCA.Rmd

Overview

This vignette shows the basic workflow of msPCA on the built-in mtcars dataset. We compute sparse principal components, inspect the solution with the print() and summary() S3 methods, and compare the sparse result with dense PCA.

Install and load

Install the package directly from CRAN.

Fit two sparse PCs

We work with the correlation matrix of mtcars and ask for two 4-sparse principal components under the default orthogonality constraint.

Sigma <- cor(mtcars)

set.seed(42)
res <- mspca(Sigma, r = 2, ks = c(4, 4), verbose = FALSE)

print() shows the sparse loading matrix restricted to the union of all active variables, together with the percentage of variance explained and the number of non-zero loadings per component.

print(res)
#> 
#> 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.4994875  0.0000000
#> cyl  -0.4952715  0.0000000
#> disp -0.5096593  0.0000000
#> hp    0.0000000 -0.5180511
#> wt   -0.4954450  0.0000000
#> qsec  0.0000000  0.5056635
#> vs    0.0000000  0.4935976
#> carb  0.0000000 -0.4819636

summary() gives a fuller breakdown: a per-PC table of variance explained, sparsity, and each component’s largest violation against any other component, followed by the full pairwise violation matrix, which shows how well the constraint is satisfied for each pair of components. The header states which constraint definition is reported: the one used to fit.

summary(res)
#> 
#> msPCA summary: 2 sparse PC(s)
#> Input type   : Sigma 
#> Runtime (s)  : 0.006 
#> Constraint   : orthogonality (as fitted) 
#> 
#> Per-component statistics:
#>   PC nonzero variance       fve cumulative_fve max_violation
#>  PC1       4 3.570419 0.3245835      0.3245835             0
#>  PC2       4 3.077834 0.2798031      0.6043866             0
#> 
#> Pairwise orthogonality violations (upper triangle):
#>     PC1 PC2
#> PC1   .   0
#> PC2   .   .
#> Total: 0e+00 
#> Violations under the other definition are stored in `$nonredundancy$uncorrelatedness`.

Working from the raw data matrix

By default (type = "Sigma") the first argument is a covariance/correlation matrix. Set type = "X" to pass the raw data matrix instead (rows are observations, columns are variables). With type = "X", msPCA applies the algorithm to the data directly: each matrix–vector product 𝚺𝜷=𝑿⊤(𝑿𝜷)/(n−1)\boldsymbol{\Sigma}\boldsymbol{\beta} = \boldsymbol{X}^\top(\boldsymbol{X}\boldsymbol{\beta})/(n - 1) is computed without ever forming the p×pp \times p matrix. This is mathematically equivalent but more scalable when p≫np \gg n.

The preprocessing arguments control which matrix is implicitly used:

  • center (default TRUE) subtracts column means.
  • scale (default TRUE) divides by column standard deviations; with scale = TRUE the algorithm operates on the correlation matrix, with scale = FALSE on the covariance matrix.
  • divisor selects the normalization: "n-1" (default, matching cov/cor) or "n".

With scale = TRUE and divisor = "n-1", the raw-data call targets exactly the same problem as the correlation-matrix call above, and returns the same solution here.

X <- as.matrix(mtcars)

set.seed(42)
res_X <- mspca(X, r = 2, ks = c(4, 4), type = "X", scale = TRUE, 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.4994875  0.0000000
#> cyl  -0.4952714  0.0000000
#> disp -0.5096593  0.0000000
#> hp    0.0000000 -0.5180511
#> wt   -0.4954451  0.0000000
#> qsec  0.0000000  0.5056627
#> vs    0.0000000  0.4935983
#> carb  0.0000000 -0.4819636

The same dual interface is available for the single-component tpm().

Orthogonality versus zero correlation

Sparse loading vectors are not automatically non-redundant, so mspca() imposes an explicit constraint between components. The default (feasibilityConstraintType = 0) enforces orthogonality of the loading vectors. Setting feasibilityConstraintType = 1 instead enforces zero pairwise correlation between the resulting scores. The choice can lead to different solutions when the variables are strongly correlated.

set.seed(42)
res_corr <- mspca(Sigma, r = 2, ks = c(4, 4),
                  feasibilityConstraintType = 1, verbose = FALSE)
print(res_corr)
#> 
#> msPCA solution: 2 sparse PCs
#> Pct. variance explained: 24.59905 22.95676 
#> Non-zero loadings per PC: 4 4 
#> 
#> Sparse PCs
#>            [,1]       [,2]
#> hp   -0.3070017  0.0000000
#> drat  0.0000000  0.3396703
#> wt    0.0000000 -0.0926711
#> qsec  0.6767011  0.0000000
#> vs    0.2739785  0.0000000
#> am    0.0000000  0.6236279
#> gear  0.0000000  0.6979430
#> carb -0.6105418  0.0000000
summary(res_corr)
#> 
#> msPCA summary: 2 sparse PC(s)
#> Input type   : Sigma 
#> Runtime (s)  : 0.125 
#> Constraint   : uncorrelatedness (as fitted) 
#> 
#> Per-component statistics:
#>   PC nonzero variance       fve cumulative_fve max_violation
#>  PC1       4 2.705895 0.2459905      0.2459905  9.559826e-05
#>  PC2       4 2.525243 0.2295676      0.4755581  9.559826e-05
#> 
#> Pairwise uncorrelatedness violations (upper triangle):
#>     PC1          PC2
#> PC1   . 9.559826e-05
#> PC2   .            .
#> Total: 9.559826e-05 
#> Violations under the other definition are stored in `$nonredundancy$orthogonality`.

Both sets of violations are computed at fit time and stored, so the solution can be inspected under the other definition without a refit:

res_corr$nonredundancy$uncorrelatedness   # the constraint that was enforced
#>     PC1          PC2
#> PC1  NA 9.559826e-05
#> PC2  NA           NA
res_corr$nonredundancy$orthogonality      # the one that was not
#>     PC1 PC2
#> PC1  NA   0
#> PC2  NA  NA

Diagnostics

The utility functions feasibility_violation_off() and fraction_variance_explained() can be called directly for custom reporting or for comparing solutions across methods — in particular for scoring loadings that did not come from mspca(), such as those of a competing package, which carry no stored diagnostics.

# Orthogonality and zero-correlation violations for the default solution
feasibility_violation_off(Sigma, res$x_best, feasibilityConstraintType = 0)
#> [1] 0
feasibility_violation_off(Sigma, res$x_best, feasibilityConstraintType = 1)
#> [1] 0.2123274

# The same two numbers, already stored on the fitted object
sum(res$nonredundancy$orthogonality,    na.rm = TRUE)
#> [1] 0
sum(res$nonredundancy$uncorrelatedness, na.rm = TRUE)
#> [1] 0.2123274

# Total and per-PC fraction of variance explained
fraction_variance_explained(Sigma, res$x_best)
#> [1] 0.6043866
fraction_variance_explained_perPC(Sigma, res$x_best)
#> [1] 0.3245835 0.2798031

The zero-correlation violation is normalized by the total variance tr(Σ)\mathrm{tr}(\Sigma), so that each pairwise term is a fraction of the total variance and the measure is comparable between a covariance matrix and the corresponding correlation matrix; the orthogonality violation needs no such normalization. Note that res$feasibility_violation, the quantity the solver compares against feasibilityTolerance, also includes the diagonal norm terms |∥𝒖t∥22−1|\bigl|\,\|\boldsymbol{u}_t\|_2^2 - 1\bigr|, and so is at least as large as the off-diagonal diagnostic above.

Comparison with dense PCA

The first two dense principal components explain more variance, but all variables receive non-zero loadings.

pca_res <- prcomp(mtcars, scale. = TRUE)
fraction_variance_explained(Sigma, pca_res$rotation[, 1:2])
#> [1] 0.8417153

Sparse PCA trades explained variance for a more interpretable loading pattern.

Where to go next