Get started with eigencore
In eigencore: Certified Partial Eigenvalue and Singular Value Computation

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 4.1,
  fig.align = "center",
  out.width = "92%",
  dpi = 100
)

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
))

# Shared palette + two small base-graphics helpers used throughout the
# vignette. Kept out of the reader's view: the reader sees the API call and
# the picture, never the plotting boilerplate.
ec_blue <- "#2166ac"; ec_red <- "#b2182b"; ec_tol <- "#ef8a62"; ec_grey <- "grey78"

# Per-pair backward error as a stem/dot plot, against the tolerance line.
plot_backward_error <- function(cert, main = NULL, tol = cert$tolerance) {
  be <- pmax(cert$backward_error, .Machine$double.eps)
  k  <- length(be)
  col_pt <- ifelse(cert$converged, ec_blue, ec_red)
  rng  <- range(c(be, tol))
  ylim <- c(rng[1] / 6, rng[2] * 6)
  op <- par(mar = c(4, 4.6, if (is.null(main)) 1 else 2.4, 1)); on.exit(par(op))
  plot(seq_len(k), be, log = "y", type = "n", xlim = c(0.5, k + 0.5),
       ylim = ylim, xaxt = "n", bty = "n", main = main,
       xlab = "returned pair", ylab = "backward error (log scale)")
  axis(1, at = seq_len(k))
  segments(seq_len(k), ylim[1], seq_len(k), be, col = col_pt, lwd = 2)
  points(seq_len(k), be, pch = 19, cex = 1.3, col = col_pt)
  abline(h = tol, col = ec_tol, lty = 2, lwd = 2)
  text(k + 0.45, tol, sprintf("tol = %.0e", tol), col = ec_tol,
       pos = 3, cex = 0.85, xpd = NA)
}

# A computed slice (blue) drawn on top of the full dense spectrum (grey).
plot_spectrum <- function(all_vals, computed, ylab = "value", main = NULL) {
  s <- sort(all_vals, decreasing = TRUE)
  k <- length(computed)
  op <- par(mar = c(4, 4.6, if (is.null(main)) 1 else 2.4, 1)); on.exit(par(op))
  plot(seq_along(s), s, pch = 19, cex = 0.4, col = ec_grey, bty = "n",
       xlab = "rank (largest to smallest)", ylab = ylab, main = main)
  points(seq_len(k), sort(computed, decreasing = TRUE),
         pch = 19, cex = 1.2, col = ec_blue)
  legend("topright", c("full spectrum (dense reference)", "computed by eigencore"),
         pch = 19, pt.cex = c(0.7, 1.2), col = c(ec_grey, ec_blue),
         bty = "n", cex = 0.9)
}

eigencore turns spectral computation into a five-step workflow:

Build an operator describing the action of A.
Describe the problem (eigen or SVD; metric B; spectral target).
Plan a solver and inspect what it will run.
Solve.
Inspect the certificate — the numerical evidence that the result is trustworthy.

This vignette walks through the workflow on three problem classes: a standard Hermitian eigenproblem, a generalized SPD eigenproblem with a metric B, and a partial SVD. It also shows the RSpectra-compatible shim for users migrating from existing code.

For the V2 CRAN release, planner labels are part of the contract: promoted paths are native and benchmark-backed in their documented regimes, while reference, prototype, oracle, and diagnostic paths remain clearly labelled.

library(eigencore)

1. Standard Hermitian eigenproblem

Build a small symmetric matrix and ask for the five largest eigenvalues.

set.seed(1)
n <- 200
A <- crossprod(matrix(rnorm(n * n), n, n)) / n + diag(n)

eigencore wraps the matrix in a block-native operator the moment you hand it to a problem constructor. You can do this explicitly:

Aop <- as_operator(A)
Aop

The operator carries dimensions, structure tags, and a flag indicating whether the underlying storage has a native kernel (in this case dense double — yes).

Build the problem and inspect the plan

P    <- eigen_problem(A, structure = hermitian(), target = largest())
plan <- plan_solver(P, k = 5)
plan

The plan tells you exactly which kernel will be invoked. There is no silent dispatch.

Solve and read the certificate

fit <- solve(P, k = 5)
fit

fit$certificate

A passing certificate means: every requested pair has

residual ||A v - lambda v|| / s below tol, where s is the labeled norm bound (here frobenius_exact),
backward error below tol,
inter-vector orthogonality loss below the orthogonality tolerance.

The certificate is not one number — it is one verdict per returned pair. Plot the per-pair backward error and the tolerance becomes a line you can see every pair clearing.

plot_backward_error(fit$certificate)

If any pair were above the line, the certificate's failed_indices slot would name it. Here the worst pair clears the tolerance by orders of magnitude.

fit$values

These five values are a certified slice of a 200-dimensional spectrum. Plotting them against the full dense spectrum shows exactly which part of the problem you paid to compute.

all_vals <- eigen(A, symmetric = TRUE, only.values = TRUE)$values
plot_spectrum(all_vals, fit$values, ylab = "eigenvalue")

With n = 200 we asked for 5 of 200 pairs. In a production problem with n = 1e6, computing the full spectrum is impossible — the partial result is the only result, which is exactly why a certificate matters.

2. Generalized SPD eigenproblem (`A v = lambda B v`)

Pass a metric B to eigen_problem() (or to eig_partial() via the B argument).

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

The certificate's residual is ||A v - lambda B v|| / (||A|| + |lambda| ||B||), and orthogonality is measured in the B-inner product where appropriate. The norm_bound_type now reports a bound for both A and B.

3. 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

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

all_sv <- svd(M, nu = 0, nv = 0)$d
plot_spectrum(all_sv, svd_fit$d, ylab = "singular value")

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

4. 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)

res$certificate

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

Where to go next

vignette("certificates") is the deep dive on reading the numerical evidence — what each field means and what to do when a check fails.
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.
linear_operator() lets you wrap a matrix-free apply callback as a first-class operator; eigencore's planner will warn when an R-level callback is being driven from a block-native solver hot loop, so you can decide whether to invest in a native operator implementation.