---
title: "Algorithm and implementation notes"
output: rmarkdown::html_vignette
bibliography: thebib.bib
vignette: >
  %\VignetteIndexEntry{Algorithm and implementation notes}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", eval = FALSE)
```

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 @cory2022sparse.

## 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 [@d2007direct; @yuan2013truncated;
@bertsimas2022solving]. 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 [@zou2006sparse; @journee2010generalized; @lu2012augmented;
@vu2013fantope; @benidis2016orthogonal; @cory2022sparse]. 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, @yuan2013truncated], 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 [@berk2017; @behdin2021sparse]. 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** [@yuan2013truncated]
>
> **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 [@mackey2008deflation] 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 @cory2022sparse 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`
[@eddelbuettel2011rcpp] and `RcppEigen` [@bates2013fast] 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:

```{r}
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:

```{r eval = FALSE}
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