Generalized eigenvalue problems with eigencore
In eigencore: Certified Partial Eigenvalue and Singular Value Computation

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
library(eigencore)
library(Matrix)

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

A generalized eigenvalue problem asks for nonzero vectors $x$ and scalars $\lambda$ satisfying

$$ A x = \lambda B x. $$

The second matrix $B$ may define a positive-definite metric, or it may be part of a more general matrix pencil. eigencore supports both cases, along with partial solves when you need only a few eigenpairs and QZ decompositions when you need the full structure of a dense pencil.

Start with a positive-definite metric

When A is symmetric (or Hermitian) and B is symmetric positive definite, pass B directly to eig_full(). The returned vectors are normalized in the $B$ metric.

A <- diag(c(2, 8, 18))
B <- diag(c(1, 2, 3))

fit <- eig_full(A, B = B)
values(fit)
certificate(fit)$passed

Here the generalized eigenvalues are 2, 4, and 6: each diagonal entry of A is divided by the corresponding entry of B.

stopifnot(
  isTRUE(certificate(fit)$passed),
  isTRUE(all.equal(sort(Re(values(fit))), c(2, 4, 6)))
)

Compute only part of the spectrum

Use eig_partial() when you need only a few eigenpairs. A target such as smallest() states which part of the spectrum to compute.

part <- eig_partial(A, B = B, k = 2, target = smallest())
values(part)
part$method
certificate(part)$passed

stopifnot(
  isTRUE(certificate(part)$passed),
  isTRUE(all.equal(sort(Re(values(part))), c(2, 4), tolerance = 1e-7))
)

The main choice is the structure of the problem and how much of its spectrum you need:

| Problem | Function | | --- | --- | | Symmetric/Hermitian A with positive-definite B, full spectrum | eig_full(A, B = B) | | Symmetric/Hermitian A with positive-definite B, partial spectrum | eig_partial(A, B = B, k = ...) | | Dense general pencil | eig_full(A, B = B, structure = general()) | | Dense generalized Schur decomposition | generalized_schur(A, B) |

Handle a dense general pencil

Use structure = general() when B is indefinite, singular, nonsymmetric, or when the pair does not satisfy the symmetric/Hermitian positive-definite contract. This path represents eigenvalues by homogeneous coordinates $(\alpha, \beta)$, where a finite eigenvalue is $\alpha / \beta$.

A_general <- matrix(c(1, 4, 2, 3), 2, 2)
B_general <- matrix(c(2, 1, 0, -1), 2, 2)

pencil <- eig_full(A_general, B = B_general, structure = general())
values(pencil)
alpha_beta(pencil)$classification
certificate(pencil)$passed

stopifnot(
  isTRUE(certificate(pencil)$passed),
  all(alpha_beta(pencil)$classification == "finite")
)

Homogeneous coordinates are especially useful when B is singular. eigencore classifies finite, infinite, and undefined eigenvalues instead of forcing every pair into an ordinary numeric ratio.

singular <- eig_full(
  diag(c(2, 3, 0)),
  B = diag(c(1, 0, 0)),
  structure = general()
)

alpha_beta(singular)$classification
certificate(singular)$failed_indices

stopifnot(identical(
  sort(alpha_beta(singular)$classification),
  sort(c("finite", "infinite", "undefined"))
))

Use QZ when you need the Schur form

generalized_schur() computes a dense generalized Schur, or QZ, decomposition. Use values() for the finite ratios and alpha_beta() when you need the homogeneous coordinates or finite/infinite classification.

qz <- generalized_schur(A_general, B_general)
values(qz)
alpha_beta(qz)$classification
qz$method

stopifnot(
  length(values(qz)) == 2L,
  all(alpha_beta(qz)$classification == "finite")
)

For pencils with singular B, sort = "finite" or sort = "infinite" moves the requested class to the leading part of the decomposition.

qz_singular <- generalized_schur(
  diag(c(2, 3, 0)),
  diag(c(1, 0, 0)),
  sort = "infinite"
)
alpha_beta(qz_singular)$classification

Keep sparse partial problems sparse

Sparse symmetric/Hermitian problems with a positive-definite metric use eig_partial(). Set allow_dense_fallback = "never" when preserving sparsity is a hard requirement.

A_sparse <- Diagonal(x = c(1, 4, 9, 16, 25, 36))
B_sparse <- Diagonal(x = c(1, 2, 3, 4, 5, 6))

sparse_fit <- eig_partial(
  A_sparse,
  B = B_sparse,
  k = 3,
  target = smallest(),
  method = lanczos(max_subspace = 6),
  allow_dense_fallback = "never"
)

values(sparse_fit)
sparse_fit$method
certificate(sparse_fit)$passed

stopifnot(
  isTRUE(certificate(sparse_fit)$passed),
  isTRUE(all.equal(
    sort(Re(values(sparse_fit))), c(1, 2, 3), tolerance = 1e-7
  ))
)

Full-spectrum functions accept base dense matrices and reject sparse or operator inputs rather than silently converting them to dense storage.

A note for older code

The retired geigen package exposed related operations, but eigencore is not a namespace-compatible replacement. When maintaining older code, translate the mathematical operation: use eig_full() for a full generalized eigensolve, generalized_schur() for QZ, and alpha_beta() to inspect homogeneous eigenvalue coordinates.

For more on result validation, see vignette("certificates").