Certificates: reading the numerical evidence

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.

Clean: every box ticked

fit_pass$certificate$passed
#> [1] TRUE
fit_pass$certificate$norm_bound_type
#> [1] "frobenius_exact+identity_exact"
fit_pass$certificate$scale_is_estimate
#> [1] FALSE

This is the easy case. Every backward error is below tolerance, orthogonality is near machine precision, and the norm used to scale the backward error is the exact Frobenius norm — no stochastic component. The absolute residuals need not themselves be below the dimensionless backward-error tolerance.

Failed: backward error 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
#> [1] FALSE
fit_fail$certificate$failed_indices
#>  [1]  1  2  3  4  5  6  7  8  9 10
fit_fail$certificate$max_backward_error
#> [1] 0.001240683

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
#> [1] TRUE

Worst backward error per restart for the ten smallest eigenpairs. A tight budget (red) stalls above the tolerance; a generous one (blue) drives the error under the line and the certificate passes.

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

norm_bound_type names what the certificate used as the scale in the backward-error ratio. The common values, from most to least trustworthy:

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

A matrix-free operator (one wrapped via linear_operator() with no exact norm metadata) forces eigencore onto that last, 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 — this is what scale_is_estimate flags:

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
#> [1] "frobenius_hutchinson_estimate+identity_exact"
fit_mf$certificate$scale_is_estimate
#> [1] TRUE
fit_mf$certificate$passed
#> [1] FALSE
fit_mf$certificate$notes
#> [1] "certificate scale uses a stochastic norm estimate; passed is withheld"

This callback operator does not expose exact norm metadata, so eigencore falls back to a Hutchinson stochastic estimate. Accordingly, the residual and orthogonality checks are consistent with convergence, but passed is withheld because the certificate scale is estimated rather than exact.

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
#> [1] "frobenius_exact+frobenius_exact"
fit_gen$certificate$max_orthogonality_loss
#> [1] 2.220446e-15
fit_gen$certificate$passed
#> [1] TRUE

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.