Reuse an eigenspace across a parameter sweep • eigencore

When you solve a sequence of nearby Hermitian eigenproblems, the eigenvectors from one solve are often good starting directions for the next. Passing those directions through initial_subspace can reduce operator work while preserving the ordinary result contract: eigencore still recomputes the Ritz values, residuals, orthogonality checks, and certificate for the current operator.

This vignette follows the smallest eigenpairs of A - rho B as rho changes. At the end you will have a certified result at every step, provenance showing whether the supplied start was used, and operator-column counts for comparing the warm and cold solves.

library(eigencore)

What problem are we solving?

Here B is a perturbation term in the standard eigenproblem, not the optional generalized-eigenproblem metric argument of eig_partial(). Both A and B are symmetric, so every member of the family A - rho B is Hermitian.

n <- 180L
k <- 5L
A <- diag(2, n)
A[row(A) == col(A) + 1L] <- -1
A[row(A) + 1L == col(A)] <- -1
x <- seq_len(n) / (n + 1)
B <- diag(0.5 + x + 0.2 * sin(4 * pi * x))

Use an explicit Lanczos method for a warm start. This prevents auto() from routing a structured or nearest-target problem to shift-invert, which does not consume initial_subspace.

method <- lanczos(
  block = k,
  max_subspace = 8L * k,
  max_restarts = 200L
)
first <- eig_partial(
  A, k = k, target = smallest(), method = method,
  tol = 1e-8, seed = 2026
)
certificate(first)$passed
#> [1] TRUE

How do you reuse the previous eigenspace?

Pass the returned eigenvectors to the next solve. The values alone are not enough: initial_subspace expects directions with one row per operator dimension.

rho <- 0.02
warm <- eig_partial(
  A - rho * B, k = k, target = smallest(), method = method,
  tol = 1e-8, seed = 2026, initial_subspace = vectors(first)
)
certificate(warm)$passed
#> [1] TRUE

The same pattern extends to a sweep. This example also runs a cold solve at each step so the work and answers can be compared on equal terms.

rhos <- c(0.02, 0.04, 0.06)
start <- vectors(first)
rows <- vector("list", length(rhos))

for (i in seq_along(rhos)) {
  op <- A - rhos[[i]] * B
  cold <- eig_partial(
    op, k, smallest(), method = method, tol = 1e-8, seed = 2026
  )
  warm <- eig_partial(
    op, k, smallest(), method = method, tol = 1e-8, seed = 2026,
    initial_subspace = start
  )
  rows[[i]] <- data.frame(
    rho = rhos[[i]],
    cold_columns = work(cold)$operator_columns,
    warm_columns = work(warm)$operator_columns,
    max_value_difference = max(abs(sort(values(cold)) - sort(values(warm))))
  )
  start <- vectors(warm)
}
comparison <- do.call(rbind, rows)

rho	cold operator columns	warm operator columns	max \|cold value - warm value\|
0.02	1160	405	1.75e-13
0.04	1040	365	1.37e-13
0.06	960	565	1.70e-13

On this reproducible family, every continuation step uses fewer operator columns than its cold counterpart and agrees on the requested values well within the solve tolerance. work() separates solve columns from current- certificate columns, adjoint work, metric work, and preconditioner work; operator_block_calls is a different metric because one block call may apply the operator to several columns. The legacy matvecs field remains unchanged and should not be compared across solver families.

Do not assume the same speedup for every sweep. When the target eigenspace changes abruptly or the supplied directions have little overlap with it, a warm solve can cost about as much as a cold one. The reproducible benchmark inst/benchmarks/bench-warm-start-continuation.R includes both high-overlap continuation and a deliberate overlap-loss jump.

When should you use a restart state?

initial_subspace is the direct basis-hint interface. Use a restart state when the basis should be a versioned, integrity-checked artifact with explicit identity, invalidation, transition, serialization, and retained-memory records. Restart-state execution goes through an executable plan.

state_plan <- plan_solver(
  eigen_problem(A - 0.02 * B, target = smallest()),
  k = k,
  method = method,
  tol = 1e-8
)
state_first <- solve(state_plan, retain_state = "same_operator")
state <- restart_state(state_first, retention = "same_operator")

same_operator <- solve(
  state_plan,
  restart_state = state,
  reuse = "same_operator"
)

next_plan <- plan_solver(
  eigen_problem(A - 0.04 * B, target = smallest()),
  k = k,
  method = method,
  tol = 1e-8
)
changed_operator <- solve(
  next_plan,
  restart_state = state,
  reuse = "auto"
)

data.frame(
  solve = c("same operator", "changed operator"),
  relation = c(
    same_operator$state_transition$relation,
    changed_operator$state_transition$relation
  ),
  basis_used = c(
    same_operator$state_transition$basis_used,
    changed_operator$state_transition$basis_used
  ),
  method_state_used = c(
    same_operator$state_transition$method_state_used,
    changed_operator$state_transition$method_state_used
  )
)
#>              solve         relation basis_used method_state_used
#> 1    same operator    same_operator       TRUE              TRUE
#> 2 changed operator changed_operator       TRUE             FALSE

The exact-revision solve may consume the retained fitted start block. The changed matrix has a new built-in operator identity, so auto keeps only the coordinate-compatible public basis and invalidates every operator-dependent claim. Both results have a fresh current-operator certificate. The state owns ordinary R data only in version 1, so retained_bytes(state) is complete and reports zero native bytes.

Version 1 deliberately has no receiving adapter for SVD, generalized, shift-invert, Arnoldi, LOBPCG, structured, or dense-fallback routes. Supplying a state to one of those routes is an error; it is never treated as a silent cold start.

How can you tell what happened to the start?

diagnostics() exposes the start source and the boundary-processing counts.

d <- diagnostics(warm)
data.frame(
  start_source = d$start_source,
  supplied = d$initial_subspace$supplied,
  accepted = d$initial_subspace$accepted,
  rejected = d$initial_subspace$rejected,
  augmented = d$initial_subspace$augmented,
  rank = d$initial_subspace$rank
)
#>    start_source supplied accepted rejected augmented rank
#> 1 user_supplied        5        5        0         0    5

The supplied matrix must be numeric, finite, and have n rows. Eigencore orthonormalizes its columns and detects numerical rank at the solver boundary. If too few independent directions survive, the start block is augmented; if more directions survive than the method’s block width, a seeded rotation lets all accepted directions contribute to the fitted block. The seed therefore makes augmentation and compression reproducible as well as controlling a cold random start.

Passing initial_subspace = NULL is exactly the cold-start contract. It leaves the previous random-start sequence and result unchanged.

Which solver plans accept a warm start?

The supported surface is deliberately narrow:

Problem and plan	Warm-start status
Standard real Hermitian Lanczos, explicit dense double matrix	Native and supported
Standard real Hermitian Lanczos, `dgCMatrix`	Native and supported
Standard real Hermitian Lanczos, matrix-free operator, `block > 1`	Native callback path and supported
Standard real Hermitian Lanczos, matrix-free operator, `block = 1`	Reference-labelled path and supported
Generalized problem using the `B` argument	Rejected
Shift-invert or `nearest()` plan	Rejected
Dense fallback, nonsymmetric, or non-Lanczos plan	Rejected

Use method = lanczos() when the start is required. With method = auto(), inspect plan_solver() first: sparse algebraic-edge targets and nearest() may select a factorized shift-invert plan. Supplying a start to an unsupported plan raises an error instead of silently ignoring it.

For native block paths, lanczos(check_stride = N) with a positive integer checks convergence every N block iterations and can stop a strong warm start before a full cold-sized sweep completes. The default check_stride = 0L preserves full-sweep behavior. Mid-sweep checks do not add operator applications; enable them when workload-specific measurements show a benefit.

What does the fresh certificate prove?

A starting subspace is only a hint. Each solve certifies the returned eigenpairs against the current A - rho B; it never reuses the previous certificate.

There is one important distinction: a small residual proves that a returned pair is an eigenpair, but not that it belongs to the requested end of the spectrum. A fully supplied subspace that is already invariant at the requested tolerance is therefore discarded in favor of a cold start, and start_source records the invariant-guard decision. For reliable continuation, supply directions that overlap the requested target eigenspace and always check certificate(fit)$passed.

Where should you go next?

vignette("reproducible-solver-plans") explains frozen plans, operator identity, RDS replay, and typed execution failures.
?eig_partial documents the complete initial_subspace contract.
vignette("certificates") explains residuals, backward error, and orthogonality checks.
diagnostics(fit) reports start provenance and separates operator block calls, operator columns, and certification columns.
Run Rscript inst/benchmarks/bench-warm-start-continuation.R --strict for the non-quick continuation evidence on your own machine.