Certificates: reading the numerical evidence

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

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.
# Converged pairs are blue, failed pairs red.
plot_backward_error <- function(cert, main = NULL, tol = cert$tolerance,
                                ylim = NULL) {
  be <- pmax(cert$backward_error, .Machine$double.eps)
  k  <- length(be)
  col_pt <- ifelse(cert$converged, ec_blue, ec_red)
  if (is.null(ylim)) {
    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)
}

Every result returned by eigencore carries a certificate — a small object that says how trustworthy the numbers are, not just what they are. Reading the certificate takes 30 seconds and replaces the usual ritual of "recompute the residual yourself to be sure."

This vignette walks through the certificate's fields, the three things they measure, and what to do when one of them fails.

For the V2 CRAN release, a certificate is not just display metadata. It is the evidence boundary that keeps promoted native paths, reference fallbacks, and V3 deferrals honest.

library(eigencore)

Anatomy of a certificate

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

fit  <- eig_partial(A, k = 5, target = largest())
cert <- fit$certificate
cert

The fields, in order:

• passed — overall verdict. TRUE means every requested pair satisfies every check (residual, backward error, orthogonality) and the norm bound used to build the scale is exact, not estimated.
• tolerance — the user-requested tolerance. Defaults to 1e-8.
• certificate_type — provenance: how the certificate was constructed. "residual_backward_error" means residuals are based on direct recomputation or native solver output, scaled by the shared backward-error denominator. "uncomputed" means vectors were not requested, so no certificate is possible.
• norm_bound_type — what was used as the norm in the scale. Common values:
• frobenius_exact — the Frobenius norm of an explicit dense matrix.
• frobenius_metadata — exact Frobenius norm derived from sparse/diagonal metadata.
• identity_exact — the standard problem B = I.
• frobenius_hutchinson_estimate — a stochastic Hutchinson estimate, used for matrix-free operators that do not expose a norm.
• scale_is_estimate — TRUE whenever the norm bound is stochastic. When this is TRUE, eigencore withholds passed even if every residual is below tolerance, because the denominator of the backward-error ratio is itself uncertain.
• max_residual — the worst absolute residual ||A v - lambda v|| (or ||A v - lambda B v|| for generalized problems) across the returned basis.
• max_backward_error — the worst residual divided by the labelled scale. This is the number that has to be below tolerance for a pair to be considered converged.
• max_orthogonality_loss — the worst entry of |V* V - I| (or |V* B V - I| in B-inner-product problems).
• failed_indices — which pairs failed. Empty when passed = TRUE.
• notes — free-form provenance text, used to explain unusual states.

The max_* fields are summaries. The certificate also keeps the underlying per-pair vectors (cert$backward_error, cert$converged), so you can see at a glance whether the whole basis cleared the bar or just barely missed:

plot_backward_error(cert)

Three things a certificate measures

1. Residual

For a Hermitian eigenproblem A v = lambda v, the absolute residual is ||A v_i - lambda_i v_i||. For an SVD, both the left and right relations matter, so the combined residual is

sqrt( ||A v - sigma u||^2 + ||A^T u - sigma v||^2 ).

For a generalized SPD problem A v = lambda B v, the residual is computed in the original coordinates: ||A v_i - lambda_i B v_i||. Eigencore never reports a residual computed in a transformed problem space without saying so in certificate_type.

2. Backward error

Absolute residuals can be misleading when the operator's norm is large or small. Backward error is the residual divided by a scale that captures "how big could a perturbation of A (and B) be that exactly makes (lambda, v) an eigenpair?":

eta_i = ||r_i|| / ( ||A|| + |lambda_i| ||B|| ).

A pair "converges" when eta_i <= tol. This is the criterion eigencore uses internally; the tolerance you pass to eig_partial(tol = ...) is the backward-error tolerance, not the residual tolerance.

3. Orthogonality

Iterative methods drift. After enough restarts the returned basis can lose orthogonality even when each per-pair residual is small. The certificate records max_abs(V* V - I) (or max_abs(V* B V - I)); its acceptance tolerance is sqrt(eps) by default. Even one passing residual should be treated with suspicion if orthogonality has collapsed: clustered or repeated eigenvalues will silently get returned as duplicates.

Reading common certificate states

The single most useful habit is to picture the per-pair backward error against the tolerance line. A passing certificate is "all stems below the line"; a failing one is "at least one stem above it." Here are the same ten eigenpairs of A, computed two ways:

# Largest eigenvalues are well separated -> easy, converges fast.
fit_pass <- eig_partial(A, k = 10, target = largest())

# Smallest eigenvalues are densely clustered near 1 -> a tight maxit
# budget leaves them short of tolerance.
fit_fail <- eig_partial(A, k = 10, target = smallest(), maxit = 15)

c(largest_passed  = fit_pass$certificate$passed,
  smallest_passed = fit_fail$certificate$passed)

op <- par(mfrow = c(1, 2))
ylim <- c(1e-16, 1e-1)
plot_backward_error(fit_pass$certificate, main = "largest(): passed", ylim = ylim)
plot_backward_error(fit_fail$certificate, main = "smallest(), maxit = 15: failed", ylim = ylim)
par(op)

Clean: every box ticked

fit_pass$certificate$passed
fit_pass$certificate$norm_bound_type
fit_pass$certificate$scale_is_estimate

This is the easy case. Every residual is below tolerance, orthogonality is near machine precision, and the norm bound used to scale the backward error is the exact Frobenius norm — no stochastic component.

Failed: residual too large

The right-hand panel above is a genuine failure. The smallest eigenvalues of A sit in a dense cluster just above 1, so they need many more iterations than the well-separated largest ones. With maxit = 15 the solver runs out of budget before any pair converges, and the verdict flips. The failed_indices slot tells you which Ritz pairs missed:

fit_fail$certificate$passed
fit_fail$certificate$failed_indices
fit_fail$certificate$max_backward_error

You do not have to guess how far off it was, or whether it was inching toward convergence. The solver records a convergence_history; plotting the worst backward error per restart shows the failed run plateauing above the line while a generous budget drives it underneath.

fit_ok <- eig_partial(A, k = 10, target = smallest(), maxit = 40)
fit_ok$certificate$passed

h_fail <- diagnostics(fit_fail)$convergence_history
h_ok   <- diagnostics(fit_ok)$convergence_history
tol <- fit_fail$certificate$tolerance
yl  <- range(c(h_fail$max_backward_error, h_ok$max_backward_error, tol))
op <- par(mar = c(4, 4.6, 1, 1))
tryCatch({
  plot(h_fail$restart, h_fail$max_backward_error, log = "y", type = "l",
       lwd = 2, col = ec_red, bty = "n", ylim = yl,
       xlab = "restart", ylab = "max backward error (log scale)")
  lines(h_ok$restart, h_ok$max_backward_error, lwd = 2, col = ec_blue)
  abline(h = tol, col = ec_tol, lty = 2, lwd = 2)
  legend("right", c("maxit = 15 (failed)", "maxit = 40 (passed)", "tolerance"),
         col = c(ec_red, ec_blue, ec_tol), lwd = 2, lty = c(1, 1, 2),
         bty = "n", cex = 0.9, inset = 0.02)
}, finally = par(op))

What to do: increase maxit, raise tol, or — when you suspect a clustered spectrum — request a larger k (so the cluster is fully covered by the returned basis) and slice afterwards. Here, lifting maxit from 15 to 40 is enough.

Withheld: stochastic norm in the denominator

A matrix-free operator (one wrapped via linear_operator() with no exact norm metadata) forces eigencore onto a Hutchinson stochastic estimate of ||A||. Because the denominator of the backward-error ratio is then a sample, not a deterministic upper bound, eigencore withholds passed even if every sampled residual ratio is below tol:

set.seed(2)
op <- linear_operator(
  dim = c(n, n),
  apply = function(X, alpha = 1, beta = 0, Y = NULL) {
    Z <- alpha * (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 * (A %*% X)
    if (is.null(Y) || beta == 0) Z else Z + beta * Y
  },
  structure = hermitian(),
  name = "matrix-free Hermitian wrapper"
)
fit_mf <- eig_partial(op, k = 5, target = largest())
fit_mf$certificate$norm_bound_type
fit_mf$certificate$scale_is_estimate
fit_mf$certificate$passed
fit_mf$certificate$notes

eigencore does not have a clean Frobenius bound for this operator without paying for a full second pass over the matrix, so it falls back to a Hutchinson stochastic estimate. The honest output is "all evidence is consistent with convergence, but the certificate was not produced by an exact bound."

What to do: if you need a hard passed = TRUE, switch to a problem class where eigencore can carry exact norm metadata (built-in dense / dgCMatrix / ddiMatrix operators), or refine with a deterministic verification pass.

Generalized SPD: B-orthogonality matters

For A v = lambda B v, the certificate's orthogonality field is in the B-inner product: max_abs(V* B V - I). A passing certificate guarantees both the residual contract and B-orthogonality of the returned vectors — which is what downstream linear-algebra code typically needs.

set.seed(4)
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$certificate$norm_bound_type
fit_gen$certificate$max_orthogonality_loss
fit_gen$certificate$passed

If a B-orthogonality value comes back near machine precision, the B-inner-product Cholesky-QR refinement inside the solver did its job. If it comes back loose (say, 1e-4), increase maxit or lower tol — orthogonality loss is usually the first thing to surface in ill-conditioned-B problems.

Comparing eigencore certificates to RSpectra diagnostics

RSpectra::eigs_sym() returns nconv and niter but does not return residuals, backward errors, or an orthogonality measure. The eigencore shim (eigencore::eigs_sym()) returns the same RSpectra-shaped list with two added fields:

res <- eigs_sym(A, k = 5, which = "LA")
names(res)
res$certificate

Code already written against RSpectra::eigs_sym() ignores certificate and diagnostics silently; new code can opt in to certified results without changing call sites.

Cheat sheet

| You see | What it means | What to do | |-------------------------------------|--------------------------------------------------------------------------|-----------------------------------------------------------| | passed = TRUE, scale_is_estimate = FALSE | Trust the result. | Use the values/vectors. | | passed = FALSE, scale_is_estimate = FALSE, failed_indices non-empty | Some pairs hit maxit before converging. | Increase maxit, raise tol, or request a wider k. | | passed withheld, scale_is_estimate = TRUE | Norm bound is stochastic; evidence is consistent with convergence. | Use a problem class with deterministic norm metadata, or run a verification pass. | | max_orthogonality_loss near sqrt(eps) but residuals tiny | Iterative drift; clustered eigenvalues at risk of duplicates. | Increase maxit; check whether there are repeated eigenvalues. | | failed_indices is the first few pairs | Convergence has not started yet. | Almost certainly a maxit problem. | | failed_indices is the last few pairs | Convergence stalled near the requested cluster. | Increase k, request a wider window, or use lobpcg() with a preconditioner. |

The certificate is the bridge between "I called a solver" and "I have a trustworthy partial spectrum." Read it.

eigencore documentation built on July 26, 2026, 1:06 a.m.