Get started with eigencore

A matrix too large to fully diagonalize still has a handful of eigenvalues or singular values you actually need — the top few directions of variance, the leading modes of a graph, the dominant singular triplets behind a low-rank approximation. eigencore computes exactly that partial answer, and every result carries a certificate: residuals, a backward-error bound, and a pass/fail verdict, so you don’t have to recompute the check yourself to trust it.

library(eigencore)

Compute a certified partial eigendecomposition

Build a symmetric matrix and ask for its five largest eigenvalues:

set.seed(1)
n <- 200
A <- crossprod(matrix(rnorm(n * n), n, n)) / n + diag(n)
fit <- eig_partial(A, k = 5, target = largest())
fit
#> Partial eigen decomposition
#>   requested: 5 
#>   converged: 5 
#>   method: native scalar thick-restart Hermitian Lanczos 
#>   target: largest 
#>   restart:thick_restart(in_native_loop)
#>   locked: 5 
#>   max residual: 1.67046e-07 
#>   max backward error: 4.546939e-09 
#>   max orthogonality loss: 3.996803e-15 
#>   norm bound: frobenius_exact+identity_exact 
#>   scale estimated: FALSE 
#>   certificate: passed

The certificate is not a side note — fit$certificate$passed is the answer to “can I trust these five numbers?”

fit$certificate$passed
#> [1] TRUE
Scatter plot of all 200 eigenvalues sorted from largest to smallest in grey, with the five largest highlighted in blue at the top-left.

The five largest eigenvalues (blue) located within the full spectrum of A (grey). eigencore computes only the requested slice, then certifies it.

With n = 200 this computed 5 of 200 pairs. In a production problem with n = 1e6, storing all eigenvectors and computing the full spectrum is usually impractical; a certified partial result is often the useful result you can afford.

Generalized SPD eigenproblem (A v = lambda B v)

Pass a metric B to eig_partial() when the problem is A v = lambda B v rather than the standard A v = lambda v.

B <- diag(seq(1, 5, length.out = n))
fit_gen <- eig_partial(A, k = 5, target = largest(), B = B,
                       method = lobpcg(maxit = 200))
fit_gen
#> Partial eigen decomposition
#>   requested: 5 
#>   converged: 5 
#>   method: native generalized SPD LOBPCG (B-orthogonal, residual certified) 
#>   target: largest 
#>   restart: lobpcg 
#>   locked: 5 
#>   max residual: 1.695109e-08 
#>   max backward error: 1.253481e-10 
#>   max orthogonality loss: 1.332268e-15 
#>   norm bound: frobenius_exact+frobenius_exact 
#>   scale estimated: FALSE 
#>   certificate: passed

The absolute residual is ||A v - lambda B v||. The certificate divides that quantity by ||A|| + |lambda| ||B|| (and the vector norm, when needed) to form the backward error used for convergence. Orthogonality is measured in the B-inner product where appropriate. See vignette("generalized-eigenproblems") for dense pencils, singular B, and the QZ decomposition.

Partial SVD

For rectangular problems use svd_partial():

M <- matrix(rnorm(400 * 50), 400, 50)
svd_fit <- svd_partial(M, rank = 5, target = largest())
svd_fit
#> Partial SVD
#>   requested rank: 5 
#>   converged rank: 5 
#>   method: native certified Gram SVD special case 
#>   target: largest 
#>   max residual: 1.245159e-15 
#>   max backward error: 8.720147e-18 
#>   max orthogonality loss: 5.551115e-16 
#>   norm bound: frobenius_exact 
#>   scale estimated: FALSE 
#>   certificate: passed

The same mental model applies to singular values: you compute the leading few and leave the tail untouched.

Scatter plot of all 50 singular values of M sorted descending in grey, with the top five highlighted in blue.

The five leading singular values (blue) computed by eigencore, shown against the full singular-value spectrum of M (grey).

The reported method identifies the path — for very small or near-square problems eigencore may use a dense LAPACK SVD fallback rather than running its iterative Golub-Kahan kernel. Either way the certificate covers both ||A v - sigma u|| and ||A^T u - sigma v||.

RSpectra-compatible workflow

If your existing code uses RSpectra::eigs_sym(), you can call eigencore in the same shape — the return list extends RSpectra’s by adding certificate and diagnostics:

res <- eigs_sym(A, k = 5, which = "LA")
str(res, max.level = 1)
#> List of 7
#>  $ values     : num [1:5] 5.01 4.77 4.7 4.57 4.5
#>  $ vectors    : num [1:200, 1:5] 0.0227 0.0103 0.0923 -0.1479 -0.0414 ...
#>  $ nconv      : int 5
#>  $ niter      : int 60
#>  $ nops       : int 62
#>  $ certificate:List of 18
#>   ..- attr(*, "class")= chr "eigencore_certificate"
#>  $ diagnostics:List of 15
res$certificate
#> eigencore certificate
#>   passed: TRUE 
#>   tolerance: 1e-08 
#>   type: residual_backward_error 
#>   norm bound: frobenius_exact+identity_exact 
#>   scale estimated: FALSE 
#>   max residual: 2.074433e-09 
#>   max backward error: 5.646539e-11 
#>   max orthogonality loss: 2.109424e-15 
#>   orthogonality tolerance: 1.490116e-08 
#>   orthogonality required: TRUE

eigs(), eigs_sym(), and svds() accept the same which codes as RSpectra — "LM", "SM", "LA", "SA", "LR", "SR", "LI", "SI", and "BE".

Under the hood: operators, problems, and plans

eig_partial() and svd_partial() cover the common case, but every call you’ve made in this vignette runs through the same four-stage architecture, and you can drop down to it directly when you want to inspect or reuse a piece of that pipeline.

1. Operators. Any matrix is wrapped in a block-native operator the moment you hand it to a problem constructor. You can do this explicitly:

Aop <- as_operator(A)
Aop
#> <eigencore operator>
#>   name: dense_matrix 
#>   dim: 200 x 200 
#>   dtype: double 
#>   structure: hermitian

The operator carries dimensions, structure tags, and a flag indicating whether the underlying storage has a native kernel (in this case dense double — yes). vignette("sparse-pca") builds operators that center and scale a sparse matrix without densifying it — the same object type, doing more work per call.

2. Problems. eigen_problem() and svd_problem() describe what to solve: the matrix, its structure (hermitian(), general(), …), and the spectral target (largest(), smallest(), nearest(), …).

P <- eigen_problem(A, structure = hermitian(), target = largest())

3. Plans. plan_solver() reports the primary method selected from the problem’s current structure, target, and policy, together with any allowed fallback. It is an inspection object, not an executable or frozen plan: certification can still trigger a recorded fallback while solving.

plan <- plan_solver(P, k = 5)
plan
#> eigencore solver plan
#>   problem: eigen 
#>   requested: 5 
#>   target: largest 
#>   method: native scalar thick-restart Hermitian Lanczos 
#>   reasons:
#>    - structure: hermitian 
#>    - target: largest 
#>    - standard eigenproblem 
#>    - built-in dense operator has native block apply 
#>   controls:
#>    - block : 1 
#>    - max_subspace : 35 
#>    - max_restarts : 100 
#>    - reorthogonalize : TRUE 
#>   fallback: dense oracle prototype if unsupported

4. Solve. solve() on a problem is what eig_partial() calls internally. Reuse the problem object when you want to solve it again. Inspect the returned method, warnings, and restart diagnostics to see the route that actually ran.

same_fit <- solve(P, k = 5)
isTRUE(all.equal(same_fit$values, fit$values))
#> [1] TRUE

Where to go next

  • vignette("sparse-pca") — the flagship workflow: centering and scaling a sparse matrix without densifying it, and building matrix-free operators by hand.
  • vignette("certificates") is the deep dive on reading the numerical evidence — what each field means and what to do when a check fails.
  • vignette("generalized-eigenproblems") covers dense pencils, singular B, and the QZ decomposition.
  • Run help(package = "eigencore") to browse the installed help index.
  • ?certificate documents the certificate fields in detail.
  • ?plan_solver explains how operator structure, target, and method combine to choose a kernel.