---
title: 'Certificates: reading the numerical evidence'
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2.0
    css: albers.css
    includes:
      in_header: albers-header.html
params:
  family: lapis
  preset: homage
resource_files:
- albers.css
- albers.js
- albers-header.html
- fonts

vignette: |
  %\VignetteIndexEntry{Certificates: reading the numerical evidence}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-opts, include = FALSE}
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
)
```

```{r albers-classes, echo=FALSE, results='asis'}
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
))
```

```{r helpers, include = FALSE}
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.

```{r setup}
library(eigencore)
```

## Anatomy of a certificate

```{r anatomy}
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:

```{r anatomy-bars, echo = FALSE, fig.cap = "Per-pair backward error for the five returned eigenpairs. The `max_backward_error` field is simply the height of the tallest stem; `passed` is TRUE because every stem clears the tolerance line.", fig.alt = "Stem plot of backward error for five eigenpairs on a log scale, all falling well below the dashed tolerance line."}
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:

```{r pass-vs-fail-fits}
# 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)
```

```{r pass-vs-fail, echo = FALSE, fig.width = 9, fig.height = 4, fig.cap = "Same matrix, same k, two targets. Left: the ten largest eigenpairs all clear the tolerance (blue, passed). Right: the ten smallest stall above it under a tight iteration budget (red, failed).", fig.alt = "Two side-by-side stem plots of per-pair backward error. Left panel shows ten blue points below the tolerance line; right panel shows ten red points above the tolerance line."}
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

```{r clean}
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:

```{r failed}
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.

```{r convergence-fits}
fit_ok <- eig_partial(A, k = 10, target = smallest(), maxit = 40)
fit_ok$certificate$passed
```

```{r convergence, echo = FALSE, fig.cap = "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.", fig.alt = "Line plot on a log scale of backward error versus restart number. The red maxit-15 curve plateaus above the tolerance line; the blue maxit-40 curve descends below it."}
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`:

```{r withheld}
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.

```{r generalized}
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:

```{r rspectra}
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.