GPCA at Scale and Special Cases

In genpca: Generalized Principal Component Analysis

if (requireNamespace("ragg", quietly = TRUE)) knitr::opts_chunk$set(dev = "ragg_png")
if (requireNamespace("systemfonts", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) albersdown::albers_register_fonts()
if (requireNamespace("ggplot2", quietly = TRUE) && requireNamespace("albersdown", quietly = TRUE)) ggplot2::theme_set(albersdown::theme_albers(family = params$family, preset = params$preset))
knitr::opts_chunk$set(
  collapse   = TRUE,
  comment    = "#>",
  message    = FALSE,
  warning    = TRUE,
  fig.width  = 6,
  fig.height = 4,
  out.width  = "85%"
)
library(genpca)
library(Matrix)

cat(sprintf(
  paste0(
    '<script>document.addEventListener("DOMContentLoaded",function(){',
    'document.body.classList.remove("palette-red","palette-lapis","palette-ochre","palette-teal","palette-green","palette-violet","preset-homage","preset-interaction","preset-study","preset-structural","preset-adobe","preset-midnight");',
    'document.body.classList.add("palette-%s","preset-%s");',
    '});</script>'
  ),
  params$family,
  params$preset
))

This vignette walks through the choices that matter once your data outgrow the defaults: which backend to pick, when to switch to a covariance-only fit, and how to project out-of-sample observations.

Backend selection

| Method | Best for | Pros | Cons | |:--|:--|:--|:--| | eigen | Small / medium dense problems | Robust reference behaviour | Can be expensive at scale; maxeig guards dense eigendecomposition of a singular general metric; it never truncates the metric | | spectra | Few components with factorizable metrics | Usually applies a whitened operator | Dense data copy, factorization costs, and dense fallbacks | | randomized | Wide (p >> n) low-rank workloads | Fast block GEMM / SpMM path | Approximation error depends on tuning | | deflation | Few components, tight memory | Low memory footprint | Can converge slowly; monitor iteration warnings | | auto | Automatic dispatch | Chooses a backend, including deflation when a singular metric exceeds the dense guard | Heuristics may not be optimal for every regime |

The default is "eigen"; pass method = "auto" to let the heuristics pick a backend for you on larger problems.

Backends on the same problem

Compare the dense reference with the randomized approximation on a full-rank noise matrix. Its slowly decaying spectrum makes approximation error visible. These single-run timings illustrate the calls; they are not a benchmark.

set.seed(11)
n <- 150; p <- 60
X <- matrix(rnorm(n * p), n, p)

t_eig <- system.time(
  fit_eig <- genpca(X, ncomp = 8, method = "eigen",
                    preproc = multivarious::center())
)
t_rnd <- system.time(
  fit_rnd <- genpca(X, ncomp = 8, method = "randomized",
                    preproc = multivarious::center())
)
data.frame(method = c("eigen", "randomized"),
           elapsed = c(t_eig["elapsed"], t_rnd["elapsed"]),
           top_sv  = c(fit_eig$sdev[1], fit_rnd$sdev[1]),
           max_relative_error = c(0, max(abs(fit_rnd$sdev / fit_eig$sdev - 1))))

plot(fit_eig$sdev, type = "b", pch = 19, col = "grey30",
     ylim = range(c(fit_eig$sdev, fit_rnd$sdev)),
     xlab = "Component", ylab = "Singular value",
     main = "Backend comparison")
lines(fit_rnd$sdev, type = "b", pch = 21, col = "steelblue")
legend("topright", legend = c("eigen", "randomized"),
       col = c("grey30", "steelblue"), pch = c(19, 21),
       bty = "n", cex = 0.85)

The maximum relative difference here is r sprintf("%.2f%%", 100 * max(abs(fit_rnd$sdev / fit_eig$sdev - 1))). Increase oversample, n_power, or n_polish when you need a more accurate approximation, then check the accuracy and time on a representative problem.

Sparse workflow (spectra)

The spectra backend factors each metric once (a sparse Cholesky here) and runs eigencore's iterative partial SVD on the whitened operator; this is useful when few components are needed and the data copy and metric factors fit in memory:

set.seed(42)
n <- 300; p <- 200
X_sparse <- rsparsematrix(n, p, density = 0.01)

# Sparse tridiagonal row/column metrics (mild AR(1)-style coupling)
M_sp <- bandSparse(n, k = c(-1, 0, 1),
                   diagonals = list(rep(0.1, n - 1), rep(1, n), rep(0.1, n - 1)))
A_sp <- bandSparse(p, k = c(-1, 0, 1),
                   diagonals = list(rep(0.1, p - 1), rep(1, p), rep(0.1, p - 1)))

fit_sp <- genpca(X_sparse, M = M_sp, A = A_sp, ncomp = 5, method = "spectra",
                 preproc = multivarious::pass())
fit_sp$sdev

What stays sparse

There are three separate storage costs: the data, the metrics or their factors, and the matrices used by the solver.

• "deflation" can retain sparse X, M, and A and apply the residual implicitly. Use preprocessing that preserves sparsity, such as pass() here: ordinary centering generally fills implicit zeros.
• "spectra" and "randomized" make a dense copy of X. Sparse input alone therefore does not bound their data storage by its nonzero count.
• The eigen and spectra factorization paths handle diagonal metrics elementwise, dense positive definite metrics by dense Cholesky, and sparse positive definite metrics by sparse Cholesky. Sparse Cholesky can add many nonzeros: fill-in depends on graph structure and ordering.
• A singular general metric on the smaller side requires dense eigendecomposition, refused above maxeig (default 5000). It is never truncated to meet that limit. A singular large-side metric is used in products without being factored. method = "auto" can route an oversized singular small-side case to deflation.
• Spectra usually applies the whitened operator without forming it, but can materialize it for a dense fallback. Its singular large-side route forms a smaller Gram matrix. The randomized method instead works with projected blocks and metric products; maxeig is not its workspace guard.

Metric validation can itself require a sparse Cholesky probe. Banded metrics such as those above have favourable fill-in; an arbitrary spatial graph need not. Budget for the factors and possible dense workspaces as well as the original sparse inputs.

Covariance-only GPCA

When you already have the cross-product C = X' M X, genpca_cov() avoids touching the full data matrix:

set.seed(123)
n <- 100; p <- 15
X <- matrix(rnorm(n * p), n, p)
M <- diag(runif(n, 0.8, 1.2))
A <- diag(runif(p, 0.7, 1.3))
C <- t(X) %*% M %*% X
fit_cov <- genpca_cov(C, R = A, ncomp = 5, method = "gmd")
fit_cov$d

barplot(fit_cov$d, names.arg = paste0("PC", seq_along(fit_cov$d)),
        col = "grey60", border = NA, ylab = "Singular value")

Out-of-sample projection

Fit on training rows, then project held-out observations into the same component space:

set.seed(7)
X <- matrix(rnorm(200 * 30), 200, 30)
fit <- genpca(X[1:150, ], ncomp = 4,
              preproc = multivarious::center())
scores_test <- multivarious::project(fit, X[151:200, ])
head(scores_test, 4)

S_train <- multivarious::scores(fit)
plot(rbind(S_train, scores_test)[, 1:2], type = "n",
     xlab = "PC1", ylab = "PC2",
     main = "Training vs out-of-sample")
points(S_train[, 1], S_train[, 2], pch = 19, col = "grey60")
points(scores_test[, 1], scores_test[, 2], pch = 19, col = "steelblue")
legend("topright", legend = c("Train", "OOS"),
       col = c("grey60", "steelblue"), pch = 19,
       bty = "n", cex = 0.85)

Performance tips

Choose preprocessing for the analysis first, then budget its storage: a centered sparse matrix can become dense. If a metric needs repair, use repair_metric() once and inspect its report before fitting. Limit ncomp to the components you intend to use, and consider the covariance route when n is large but p is moderate.

Where next

See GPCA Metrics for building metrics, and Getting Started for a getting-started walkthrough.

genpca documentation built on Sept. 17, 2026, 1:09 a.m.