PCA on sparse data without densifying

Principal component analysis starts with centering: subtract each column’s mean before finding directions of maximum variance. For a sparse matrix, that one subtraction is the problem. A - 1 %*% t(col_means) fills in every zero with a small nonzero number, so the matrix that started sparse becomes dense before you’ve computed a single eigenvector.

eigencore solves the centered (and optionally scaled) problem as an operator, so the subtraction never happens as a stored matrix. This vignette builds one sparse dataset, centers and scales it, and computes a certified partial SVD without ever forming the dense copy.

library(eigencore)
library(Matrix)

Build a sparse dataset with a low-rank signal

Simulate the kind of matrix behind collaborative filtering and topic models: sparse interaction counts between 2,000 users and 500 items, built from five latent factors so a handful of singular vectors should recover most of the structure.

set.seed(1)
m <- 2000L; n <- 500L; k_true <- 5L
U <- qr.Q(qr(matrix(rnorm(m * k_true), m, k_true)))
V <- qr.Q(qr(matrix(rnorm(n * k_true), n, k_true)))
signal <- U %*% diag(seq(8, 4, length.out = k_true)) %*% t(V)

Only 4% of user-item pairs are observed, and the observations are noisy:

mask  <- rsparsematrix(m, n, density = 0.04)
noise <- rsparsematrix(m, n, density = 0.01, rand.x = function(k) rnorm(k, sd = 0.05))
A <- as((signal * (mask != 0)) + noise, "dgCMatrix")
dim(A)
#> [1] 2000  500

At 2,000 x 500 this is a small matrix, but the storage gap is already visible: a dense copy would need 7.6 MB, while the sparse representation needs 0.57 MB — about 13 times less.

Center without densifying

center() wraps A in an operator that subtracts column means on the fly, computing the means directly from the sparse matrix rather than requiring you to supply them.

Ac <- center(A, columns = TRUE)
Ac
#> <eigencore operator>
#>   name: center(sparse_csc_matrix) 
#>   dim: 2000 x 500 
#>   dtype: double 
#>   structure: general

Pass the centered operator straight to svd_partial(), exactly as you would the original matrix.

fit <- svd_partial(Ac, rank = 5, target = largest())
fit
#> Partial SVD
#>   requested rank: 5 
#>   converged rank: 5 
#>   method: native matrix-free Golub-Kahan callback cycle + native Ritz extraction (callback boundary) 
#>   target: largest 
#>   max residual: 5.968384e-11 
#>   max backward error: 1.040313e-11 
#>   max orthogonality loss: 2.331468e-15 
#>   norm bound: frobenius_hutchinson_estimate 
#>   scale estimated: TRUE 
#>   certificate: failed
Scatter plot of all 500 singular values sorted descending in grey, with the top five highlighted in blue.

The five largest singular values of the centered matrix (blue) against the full singular spectrum (grey).

The five requested singular values occupy the leading edge of the spectrum, although the fifth-to-sixth gap is modest after masking and noise. This example checks the non-densifying computation rather than claiming reliable recovery of the five simulated factors.

Read the certificate on a composed operator

fit$certificate$passed is FALSE here, but look at why before treating that as a problem:

fit$certificate$max_backward_error
#> [1] 1.040313e-11
fit$certificate$norm_bound_type
#> [1] "frobenius_hutchinson_estimate"
fit$certificate$scale_is_estimate
#> [1] TRUE

The backward error is tiny — far below the 1e-8 tolerance. What’s withheld is the scale: the current centered sparse operator does not propagate an exact Frobenius norm, so eigencore reports a stochastic (Hutchinson) estimate and refuses to mark the certificate passed while that estimate is the only evidence. This is the certificate bookkeeping vignette("certificates") covers in detail. It is a current metadata limitation, not a mathematical requirement of centering; supplying exact norm metadata gives a deterministic certificate scale.

Does it match what you’d get by densifying?

The whole point of center() is that it should be numerically indistinguishable from centering the dense matrix directly. Check it:

max(abs(sort(fit$d, decreasing = TRUE) - sort(all_sv[1:5], decreasing = TRUE)))
#> [1] 1.276756e-15

The two agree to machine precision. dense_centered above was built only to draw the comparison plot and run this check — the eigencore computation itself never materialized it.

Scale columns too

A common next step in PCA preprocessing is scaling each column to unit sample variance, turning a covariance-matrix PCA into a correlation-matrix one. Compute the sample variances from the sparse matrix and pass the inverse standard deviations to scale_cols(), composed with the centering operator you already have:

col_var <- (Matrix::colSums(A^2) - m * Matrix::colMeans(A)^2) / (m - 1)
w <- 1 / sqrt(col_var)
Acs <- scale_cols(Ac, w)
fit_scaled <- svd_partial(Acs, rank = 5, target = largest())
fit_scaled$d
#> [1] 76.87178 72.70075 70.03808 68.13807 66.94059

center() and scale_cols() compose freely because both return the same eigencore_operator type that every solver accepts — there’s no separate “scaled matrix” object to keep track of.

What did the planner actually do?

Every solve is preceded by a plan, and you can inspect it directly instead of trusting a black box:

plan <- plan_solver(svd_problem(Acs), rank = 5, target = largest())
plan$method
#> [1] "native matrix-free Golub-Kahan callback cycle + native Ritz extraction (callback boundary)"
plan$reasons
#> [1] "target: largest"                                                                                                                                  
#> [2] "rectangular SVD problem"                                                                                                                          
#> [3] "adjoint is available"                                                                                                                             
#> [4] "matrix-free operator uses a native Golub-Kahan callback cycle with native Ritz extraction; sparse/matrix-free performance promotion remains gated"
#> [5] "operator uses R-level apply path in current prototype"

Centering and scaling remain non-densifying, but they are not currently one fully fused native operator. The centered sparse apply is native; the scaling layer composes through an R-level operator apply. The solver therefore drives the result through a matrix-free Golub-Kahan callback cycle rather than a fully block-native kernel — the reasons make that boundary explicit. plan_solver() takes the same problem/rank/target arguments as the solve, so you can inspect the intended primary route before paying for the computation; the result diagnostics record any fallback that actually occurs.

Going fully matrix-free

Sometimes there’s no matrix at all — just a function that knows how to apply A and its adjoint. linear_operator() wraps arbitrary callbacks the same way center() wraps a sparse matrix.

op <- linear_operator(
  dim = dim(A),
  apply = function(X, alpha = 1, beta = 0, Y = NULL) {
    # Sparse %*% dense returns an S4 dgeMatrix; as.matrix() keeps the
    # callback's return type consistent with what the native solver expects.
    Z <- alpha * as.matrix(A %*% X)
    if (is.null(Y) || beta == 0) Z else Z + beta * Y
  },
  apply_adjoint = function(X, alpha = 1, beta = 0, Y = NULL) {
    Z <- alpha * as.matrix(crossprod(A, X))
    if (is.null(Y) || beta == 0) Z else Z + beta * Y
  },
  name = "matrix-free A"
)
fit_mf <- svd_partial(op, rank = 5, target = largest())
fit_mf$certificate$norm_bound_type
#> [1] "frobenius_hutchinson_estimate"

Without any norm metadata, this is another stochastic-estimate certificate. Supply the exact Frobenius norm — cheap to compute once for a sparse matrix — and the certificate stops withholding passed:

op_hinted <- linear_operator(
  dim = dim(A),
  apply = op$apply,
  apply_adjoint = op$apply_adjoint,
  name = "matrix-free A (exact norm)",
  metadata = list(frobenius_norm = norm(A, type = "F"))
)
fit_hinted <- svd_partial(op_hinted, rank = 5, target = largest())
fit_hinted$certificate$norm_bound_type
#> [1] "frobenius_metadata"
fit_hinted$certificate$passed
#> [1] TRUE

Whether you reach for center()/scale_cols() or hand-write a linear_operator(), the certificate always tells you which kind of evidence backed the result.

Why this matters at scale

The 2,000 x 500 example above keeps this vignette fast to build, but the memory argument only gets sharper as the matrix grows. A production recommender with 5 million users and 200,000 items would need roughly 8 TB to store as a dense double matrix. At densities from 0.01% to 0.05%, the sparse value and row-index arrays alone require roughly 1.2 to 6 GB, before format overhead. Densifying to center it would defeat the point of using a sparse format.

Where to go next