# Worked example: msPCA on mtcars

## 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.

    install.packages("msPCA")

    library(msPCA)

## 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.009
    #> 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 \(\boldsymbol{\Sigma}\boldsymbol{\beta} = \boldsymbol{X}^\top(\boldsymbol{X}\boldsymbol{\beta})/(n - 1)\) is computed without ever forming the \(p \times p\) matrix. This is mathematically equivalent but more scalable when \(p \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.189
    #> 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 \(\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 \(\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

- `vignette("case-study-snp500", package = "msPCA")` applies the same workflow to a 423-stock correlation matrix shipped with the package, and shows how the two constraint types can yield noticeably different factor structures.
- `vignette("algorithm-and-implementation", package = "msPCA")` documents the optimization problem, both algorithms, the implementation, and how to set `ks`, `feasibilityConstraintType` and the iteration budgets.
- A comparison against eight competing functions, drawn from seven packages, on four datasets is available on the package website, under Benchmarking against other sparse PCA packages.