inst/doc/certificates.R

In eigencore: Certified Partial Eigenvalue and Singular Value Computation

params <-
list(family = "lapis", preset = "homage")

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

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

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

## ----setup--------------------------------------------------------------------
library(eigencore)

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

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

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

## ----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--------------------------------------------------------------------
fit_pass$certificate$passed
fit_pass$certificate$norm_bound_type
fit_pass$certificate$scale_is_estimate

## ----failed-------------------------------------------------------------------
fit_fail$certificate$passed
fit_fail$certificate$failed_indices
fit_fail$certificate$max_backward_error

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

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

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

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

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