Benchmark smoke tests

In eigencore: Certified Partial Eigenvalue and Singular Value Computation

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

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

`%||%` <- 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.

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:

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

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

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.

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.")
)

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

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.

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

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:

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:

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.

eigencore documentation built on July 26, 2026, 1:06 a.m.