GPCA at Scale and Special Cases

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))))
#>       method elapsed   top_sv max_relative_error
#> 1      eigen   0.002 19.48896         0.00000000
#> 2 randomized   0.005 19.19057         0.02819281
The randomized approximation underestimates the reference singular values on this full-rank example. The table reports the largest relative difference.

The randomized approximation underestimates the reference singular values on this full-rank example. The table reports the largest relative difference.

The maximum relative difference here is 2.82%. 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
#> [1] 5.153523 4.533004 4.258609 4.174038 3.991699

What stays sparse

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

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
#> [1] 13.80217 12.42550 11.92054 11.15895 10.96560
Singular values from the covariance-only fit.

Singular values from the covariance-only fit.

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)
#>            PC1        PC2        PC3        PC4
#> [1,] 1.9323426  0.4080526  0.1924407  0.6673048
#> [2,] 0.3498745 -0.6490485 -0.2579339 -0.9996621
#> [3,] 0.9590113 -0.9305126  1.4603670  1.1496702
#> [4,] 0.1125973 -1.0845757  0.2493420  1.2509707
Training scores (grey) and out-of-sample scores (blue) projected into the same component space.

Training scores (grey) and out-of-sample scores (blue) projected into the same component space.

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.