---
title: 'Certificates: reading the numerical evidence'
output:
  rmarkdown::html_vignette:
    toc: yes
    toc_depth: 2
    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}
# Shared plotting helper. 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"

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

```{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 you act on, in order:

- **`passed`** — overall verdict. `TRUE` means every requested pair's backward
  error is at or below `tolerance`, the required orthogonality check passes,
  and the norm used to build the scale is exact rather than estimated. The
  absolute residual is recorded, but it is not compared directly with
  `tolerance`.
- **`tolerance`** — the user-requested tolerance. Defaults to `1e-8`.
- **`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`.

Two more fields explain *how* the verdict was reached rather than what it
is: `norm_bound_type` names the scale used to compute backward error, and
`scale_is_estimate` flags when that scale is itself a stochastic estimate
rather than an exact bound — the "Withheld" section below shows exactly
what that means for `passed`. `certificate_type` and `notes` carry
provenance text for unusual states. `?certificate` documents every field.

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 falls below 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 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:

```{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

`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:

```{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
```

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.

```{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 returned pairs exceed the backward-error tolerance.                    | Inspect convergence history, then consider more iterations, a different tolerance, or 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` contains particular positions | Those returned pairs missed the tolerance; index position alone does not identify the cause. | Inspect `convergence_history` and the requested target before changing solver controls. |

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

## Where to go next

- `vignette("sparse-pca")` shows a realistic withheld certificate in
  context — the current centered sparse operator does not propagate an exact
  norm, so certification uses a stochastic estimate unless exact metadata is
  supplied.
- `vignette("eigencore")` and `vignette("generalized-eigenproblems")` cover
  the workflows that produce the certificates read here.