---
title: Benchmark smoke tests
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{Benchmark smoke tests}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup-opts, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 8,
  fig.height = 4.6,
  fig.align = "center",
  out.width = "96%",
  dpi = 100
)

benchmark_ready <- all(vapply(
  c("bench", "RSpectra", "irlba"),
  requireNamespace,
  logical(1),
  quietly = TRUE
))
benchmark_iterations <- 3L
options(knitr.kable.NA = "not run")
```

```{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}
`%||%` <- function(x, y) if (is.null(x)) y else x

fro_norm <- function(A) {
  if (inherits(A, "sparseMatrix")) {
    sqrt(sum(A@x^2))
  } else {
    sqrt(sum(abs(A)^2))
  }
}

relative_error <- function(x, truth) {
  x <- sort(Re(x), decreasing = TRUE)
  truth <- sort(Re(truth), decreasing = TRUE)
  max(abs(x - truth) / pmax(1, abs(truth)))
}

eigen_backward_error <- function(A, values, vectors) {
  A_dense <- as.matrix(A)
  values <- Re(values)
  vectors <- as.matrix(vectors)
  residual <- A_dense %*% vectors - vectors %*% diag(values, nrow = length(values))
  scale <- fro_norm(A_dense) + abs(values)
  max(sqrt(colSums(abs(residual)^2)) / pmax(scale, .Machine$double.eps))
}

svd_backward_error <- function(A, d, u, v) {
  A_dense <- as.matrix(A)
  d <- Re(d)
  u <- as.matrix(u)
  v <- as.matrix(v)
  left_residual <- A_dense %*% v - u %*% diag(d, nrow = length(d))
  right_residual <- t(A_dense) %*% u - v %*% diag(d, nrow = length(d))
  residual <- sqrt(colSums(abs(left_residual)^2) + colSums(abs(right_residual)^2))
  max(residual / pmax(fro_norm(A_dense) + d, .Machine$double.eps))
}

bench_eval <- function(expr, iterations = benchmark_iterations) {
  expr <- substitute(expr)
  env <- parent.frame()
  last_result <- NULL
  mark <- tryCatch(
    bench::mark(
      last_result <- eval(expr, env),
      iterations = iterations,
      check = FALSE,
      memory = isTRUE(capabilities("profmem")),
      filter_gc = FALSE
    ),
    error = function(e) e
  )
  if (inherits(mark, "error")) {
    return(list(
      result = NULL,
      median_ms = NA_real_,
      mem_mb = NA_real_,
      error = conditionMessage(mark)
    ))
  }
  list(
    result = last_result,
    median_ms = as.numeric(mark$median[[1L]]) * 1000,
    mem_mb = as.numeric(mark$mem_alloc[[1L]]) / 1024^2,
    error = NA_character_
  )
}

make_dense_hermitian <- function(n, seed = 1L) {
  set.seed(seed)
  X <- matrix(rnorm(n * n), n, n)
  crossprod(X) / n + diag(seq(1, 1.2, length.out = n))
}

make_dense_low_rank <- function(m, n, rank = 8L, noise = 1e-3, seed = 1L) {
  set.seed(seed)
  U <- qr.Q(qr(matrix(rnorm(m * rank), m, rank)))
  V <- qr.Q(qr(matrix(rnorm(n * rank), n, rank)))
  signal <- U %*% diag(seq(rank, 1, length.out = rank), nrow = rank) %*% t(V)
  signal + noise * matrix(rnorm(m * n), m, n)
}

path_laplacian <- function(n) {
  Matrix::bandSparse(
    n,
    k = c(-1L, 0L, 1L),
    diagonals = list(rep(-1, n - 1L), c(1, rep(2, n - 2L), 1), rep(-1, n - 1L))
  )
}

eigen_rows <- function(name, A, k = 6L, tol = 1e-8, seed = 1L,
                       iterations = benchmark_iterations) {
  truth <- eigen(as.matrix(A), symmetric = TRUE, only.values = TRUE)$values[seq_len(k)]
  methods <- c("eigencore", "RSpectra", "base")
  rows <- lapply(methods, function(method) {
    timed <- switch(
      method,
      eigencore = bench_eval({
        set.seed(seed)
        eig_partial(A, k = k, target = largest(), tol = tol)
      }, iterations = iterations),
      RSpectra = bench_eval({
        RSpectra::eigs_sym(A, k = k, which = "LA", opts = list(tol = tol, maxitr = 1000L))
      }, iterations = iterations),
      base = bench_eval({
        eigen(as.matrix(A), symmetric = TRUE)
      }, iterations = iterations)
    )
    if (!is.na(timed$error)) {
      return(data.frame(
        regime = name,
        task = "eigen",
        method = method,
        median_ms = timed$median_ms,
        mem_mb = timed$mem_mb,
        rel_error = NA_real_,
        backward_error = NA_real_,
        residual_check = FALSE,
        eigencore_label = NA_character_,
        status = timed$error,
        stringsAsFactors = FALSE
      ))
    }
    result <- timed$result
    extracted <- switch(
      method,
      eigencore = list(values = values(result), vectors = vectors(result)),
      RSpectra = list(values = result$values, vectors = result$vectors),
      base = {
        ord <- order(result$values, decreasing = TRUE)[seq_len(k)]
        list(values = result$values[ord], vectors = result$vectors[, ord, drop = FALSE])
      }
    )
    backward <- eigen_backward_error(A, extracted$values, extracted$vectors)
    data.frame(
      regime = name,
      task = "eigen",
      method = method,
      median_ms = timed$median_ms,
      mem_mb = timed$mem_mb,
      rel_error = relative_error(extracted$values, truth),
      backward_error = backward,
      residual_check = isTRUE(backward <= tol),
      eigencore_label = if (method == "eigencore") result$method else NA_character_,
      status = "ok",
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

svd_rows <- function(name, A, rank = 6L, tol = 1e-8, seed = 1L,
                     iterations = benchmark_iterations) {
  truth <- svd(as.matrix(A), nu = 0, nv = 0)$d[seq_len(rank)]
  methods <- c("eigencore", "RSpectra", "irlba", "base")
  rows <- lapply(methods, function(method) {
    timed <- switch(
      method,
      eigencore = bench_eval({
        set.seed(seed)
        svd_partial(A, rank = rank, target = largest(), tol = tol)
      }, iterations = iterations),
      RSpectra = bench_eval({
        RSpectra::svds(A, k = rank, nu = rank, nv = rank, opts = list(tol = tol, maxitr = 1000L))
      }, iterations = iterations),
      irlba = bench_eval({
        set.seed(seed)
        irlba::irlba(A, nv = rank, nu = rank, tol = tol)
      }, iterations = iterations),
      base = bench_eval({
        svd(as.matrix(A), nu = rank, nv = rank)
      }, iterations = iterations)
    )
    if (!is.na(timed$error)) {
      return(data.frame(
        regime = name,
        task = "SVD",
        method = method,
        median_ms = timed$median_ms,
        mem_mb = timed$mem_mb,
        rel_error = NA_real_,
        backward_error = NA_real_,
        residual_check = FALSE,
        eigencore_label = NA_character_,
        status = timed$error,
        stringsAsFactors = FALSE
      ))
    }
    result <- timed$result
    extracted <- switch(
      method,
      eigencore = list(d = values(result), u = left_vectors(result), v = right_vectors(result)),
      RSpectra = list(d = result$d, u = result$u, v = result$v),
      irlba = list(d = result$d, u = result$u, v = result$v),
      base = list(d = result$d[seq_len(rank)], u = result$u, v = result$v)
    )
    backward <- svd_backward_error(A, extracted$d, extracted$u, extracted$v)
    data.frame(
      regime = name,
      task = "SVD",
      method = method,
      median_ms = timed$median_ms,
      mem_mb = timed$mem_mb,
      rel_error = relative_error(extracted$d, truth),
      backward_error = backward,
      residual_check = isTRUE(backward <= tol),
      eigencore_label = if (method == "eigencore") result$method else NA_character_,
      status = "ok",
      stringsAsFactors = FALSE
    )
  })
  do.call(rbind, rows)
}

metric_table <- function(rows, metric, digits = 4L) {
  regimes <- unique(rows$regime)
  methods <- c("eigencore", "RSpectra", "irlba", "base")
  out <- data.frame(regime = regimes, check.names = FALSE)

  for (method in methods) {
    out[[method]] <- vapply(regimes, function(regime) {
      keep <- rows$regime == regime & rows$method == method & rows$status == "ok"
      if (!any(keep)) NA_real_ else rows[[metric]][which(keep)[1L]]
    }, numeric(1))
    out[[method]] <- signif(out[[method]], digits)
  }

  out
}

timing_table <- function(rows) {
  out <- metric_table(rows, "median_ms")
  method_columns <- setdiff(names(out), "regime")
  out$lowest_median <- vapply(out$regime, function(regime) {
    current <- rows[
      rows$regime == regime & rows$status == "ok",
      c("method", "median_ms"),
      drop = FALSE
    ]
    if (!nrow(current)) {
      return(NA_character_)
    }
    current$method[[which.min(current$median_ms)]]
  }, character(1))
  out[, c("regime", "lowest_median", method_columns)]
}

quality_table <- function(rows) {
  regimes <- unique(rows$regime)
  do.call(rbind, lapply(regimes, function(regime) {
    current <- rows[rows$regime == regime, , drop = FALSE]
    ok <- current[current$status == "ok", , drop = FALSE]
    if (!nrow(ok)) {
      return(data.frame(
        regime = regime,
        successful_methods = sprintf("0/%d", nrow(current)),
        max_relative_error = NA_character_,
        max_backward_error = NA_character_,
        residual_checks = "not available",
        stringsAsFactors = FALSE
      ))
    }
    data.frame(
      regime = regime,
      successful_methods = sprintf("%d/%d", nrow(ok), nrow(current)),
      max_relative_error = formatC(max(ok$rel_error), format = "e", digits = 2),
      max_backward_error = formatC(max(ok$backward_error), format = "e", digits = 2),
      residual_checks = if (all(ok$residual_check)) "all pass" else "one or more fail",
      stringsAsFactors = FALSE
    )
  }))
}
```

This vignette runs five small eigenvalue and SVD cases with eigencore,
`RSpectra`, `irlba`, and base R. The matrices are intentionally small and each
method is timed three times so the article remains practical to rebuild.
Results describe one render on one machine; they are not evidence of package
rankings or scaling behavior.

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

## What is included?

Dense and sparse inputs exercise different code paths, and SVD has different
comparison methods than Hermitian eigenproblems. The smoke suite covers these
five cases:

```{r regimes}
benchmark_regimes <- data.frame(
  regime = c(
    "dense Hermitian",
    "sparse path Laplacian",
    "dense low-rank SVD",
    "tall sparse SVD",
    "wide sparse SVD"
  ),
  input = c("120 x 120 dense", "300 x 300 dgCMatrix", "180 x 70 dense",
            "320 x 60 dgCMatrix", "60 x 320 dgCMatrix"),
  compared_methods = c("eigencore, RSpectra, base",
                       "eigencore, RSpectra, base",
                       "eigencore, RSpectra, irlba, base",
                       "eigencore, RSpectra, irlba, base",
                       "eigencore, RSpectra, irlba, base")
)
knitr::kable(benchmark_regimes)
```

Base R supplies the dense reference values. For sparse matrices it converts
the input with `as.matrix()`, so its timing and memory results apply only to
these small cases. `irlba` appears only in the SVD cases.

## What is recorded?

Each row records four things:

- median elapsed time from `bench::mark()`;
- allocated memory reported by `bench`;
- relative error against a dense base truth for the requested values;
- an explicit residual/backward-error check on the returned vectors.

For eigencore, the planner label identifies the implementation used in each
case.

## Results from this render

```{r build-cases, eval = benchmark_ready, include = FALSE}
set.seed(1001)
bench_cases <- list(
  dense_hermitian = make_dense_hermitian(120L, seed = 1001L),
  sparse_laplacian = path_laplacian(300L),
  dense_low_rank_svd = make_dense_low_rank(180L, 70L, rank = 8L, seed = 1002L),
  tall_sparse_svd = Matrix::rsparsematrix(320L, 60L, density = 0.035),
  wide_sparse_svd = Matrix::rsparsematrix(60L, 320L, density = 0.035)
)
```

```{r run-benchmarks, eval = benchmark_ready, include = FALSE}
benchmark_rows <- rbind(
  eigen_rows("dense Hermitian", bench_cases$dense_hermitian,
             k = 6L, seed = 2001L),
  eigen_rows("sparse path Laplacian", bench_cases$sparse_laplacian,
             k = 6L, seed = 2002L),
  svd_rows("dense low-rank SVD", bench_cases$dense_low_rank_svd,
           rank = 6L, seed = 2003L),
  svd_rows("tall sparse SVD", bench_cases$tall_sparse_svd,
           rank = 6L, seed = 2004L),
  svd_rows("wide sparse SVD", bench_cases$wide_sparse_svd,
           rank = 6L, seed = 2005L)
)
stopifnot(all(is.finite(benchmark_rows$median_ms[benchmark_rows$status == "ok"])))
stopifnot(all(is.finite(benchmark_rows$rel_error[benchmark_rows$status == "ok"])))
stopifnot(all(is.finite(benchmark_rows$backward_error[benchmark_rows$status == "ok"])))
```

### Elapsed time

The entries are median solver-call milliseconds; lower values are better. An
eigencore solver call includes its native certificate. The external calls do
not, and the shared residual check used elsewhere in this vignette runs after
the timer. These smoke timings therefore describe call overhead, not a common
time-to-certified-answer. The installed benchmark scripts report both raw
solver time and a common eigencore-certified total. The `lowest_median` column
is selected from unrounded measurements and reports the smallest value in each
row.

```{r timing-table, eval = benchmark_ready, echo = FALSE}
timing_rows <- timing_table(benchmark_rows)
knitr::kable(
  timing_rows,
  col.names = c("case", "lowest median", "eigencore", "RSpectra", "irlba", "base R"),
  align = c("l", "l", "r", "r", "r", "r"),
  caption = paste("Median solver-call time in milliseconds from",
                  benchmark_iterations, "iterations per method.")
)
```

```{r timing-summary, eval = benchmark_ready, echo = FALSE, results = 'asis'}
eigencore_lowest <- sum(timing_rows$lowest_median == "eigencore")
cat(sprintf(
  paste0(
    "In this render, eigencore recorded the lowest median in **%d of %d** ",
    "cases. The result is mixed, and the sub-millisecond rows are especially ",
    "sensitive to setup overhead and run-to-run noise."
  ),
  eigencore_lowest,
  nrow(timing_rows)
))
```

### Why can larger cases look different?

The smoke cases above are small enough that fixed work matters: dispatch,
result construction, certification, and the dense eigensolve for a bounded
Gram matrix can be a substantial fraction of total time. They are useful for
checking behavior, but they are poor evidence about scaling.

For a highly rectangular sparse SVD, eigencore can form the smaller Gram
problem and certify the resulting triplets in the original coordinates. If the
smaller dimension stays bounded while the long dimension grows, the Gram
eigensolve retains the same dimension while iterative methods perform more
expensive matrix-vector products. A crossover is therefore plausible. The
installed large endpoint cases show that eigencore can win in this regime, but
those endpoints do not locate a crossover because their dimensions, density,
rank, and random fixture all change together. The controlled mode of
`bench-svd-gram-cutoff.R` fixes the small side, rank, density, and a nested
sparse fixture while varying only the long side, and reports raw solver and
common certified-total ratios separately.

In the 2026-07-15 installed five-iteration sweep with small side 90, the
certified-total ratio crossed above one by long side 5,000 in both
orientations. The raw-call ratio crossed there for the wide case and by 25,000
for the tall case. These are machine-specific release measurements, not a
universal crossover constant.

This does not mean that increasing either matrix dimension helps. The current
planner limits the smaller dimension to 512 for tall matrices and 1024 for
wide matrices, and also checks aspect ratio, requested-rank fraction, and a
memory budget. The repository's cutoff evidence is mixed outside the core
regime: increasing the smaller side can remove the advantage, particularly for
tall matrices. The current full endpoint cases are `100000 x 500` and its
transpose; see the
[benchmark manifest](https://github.com/bbuchsbaum/eigencore/blob/main/docs/v1-benchmark-manifest.md)
for the installed evidence and exact commands.

The larger benchmark scripts pass the requested tolerance to RSpectra and
irlba, then recompute eigencore certificates for all external results. They
retain both raw solver time and common time through certification. A row that
still fails the requested certificate is excluded from time-to-certified-answer
claims rather than being treated as a slow certified result.

### Allocated memory

These values are megabytes allocated during the timed expression. They do not
measure peak resident memory.

```{r memory-table, eval = benchmark_ready, echo = FALSE}
memory_rows <- metric_table(benchmark_rows, "mem_mb")
knitr::kable(
  memory_rows,
  col.names = c("case", "eigencore", "RSpectra", "irlba", "base R"),
  align = c("l", "r", "r", "r", "r"),
  caption = "Allocated memory in megabytes."
)
```

### Numerical checks

The summary below reports the largest error observed across the methods in
each case. Timing should be ignored for any case that has a failed method or
residual check.

```{r quality-table, eval = benchmark_ready, echo = FALSE}
quality_rows <- quality_table(benchmark_rows)
knitr::kable(
  quality_rows,
  col.names = c("case", "methods completed", "max relative error",
                "max backward error", "residual checks"),
  align = c("l", "c", "r", "r", "l"),
  caption = "Numerical checks across all methods in each case."
)
```

### Eigencore planner paths

```{r planner-table, eval = benchmark_ready, echo = FALSE}
planner_rows <- benchmark_rows[
  benchmark_rows$method == "eigencore",
  c("regime", "eigencore_label", "status"),
  drop = FALSE
]
rownames(planner_rows) <- NULL
knitr::kable(
  planner_rows,
  col.names = c("case", "planner label", "status"),
  align = c("l", "l", "l")
)
```

## Limits of this comparison

The rows above are deliberately small. They keep the vignette runnable, and
they catch obvious problems: wrong target, broken sparse handling, missing
vectors, certification drift, and large timing regressions. They do not prove
large-scale speed claims.

For larger sparse matrices, base R is a truth oracle only while `as.matrix(A)`
is affordable. Studying scaling and memory behavior requires installed-package
runs, larger fixtures, repeated measurements, and saved artifacts.

## Run the larger benchmark scripts

The scripts under `inst/benchmarks` run larger cases with stricter checks and
saved outputs. Run them from an installed package rather than a mutable source
tree:

```sh
R CMD INSTALL --library=/tmp/eigencore-bench-lib .
R_LIBS=/tmp/eigencore-bench-lib \
  Rscript inst/benchmarks/bench-native-hermitian-gate.R \
  --strict --save --cases=path_laplacian:1000
```

For the SVD surface, compare eigencore against `RSpectra`, `irlba`, and base
rows where those baselines are meaningful:

```sh
R_LIBS=/tmp/eigencore-bench-lib \
  Rscript inst/benchmarks/bench-svd-surface.R \
  --iterations=1 --h-candidate \
  --methods=eigencore,RSpectra,irlba,base \
  --cases=tall_sparse,wide_sparse \
  --subject=eigencore --strict --save
```

Those scripts add strict pass/fail gates, larger sparse fixtures, saved RDS
artifacts, memory ratios, and planner-provenance checks. Use
`docs/v1-benchmark-manifest.md` for the current release-gate inventory and
`docs/post-v1-benchmark-gates.md` for future-promotion gates.