inst/doc/generalized-eigenproblems.R

In eigencore: Certified Partial Eigenvalue and Singular Value Computation

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

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

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

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

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

## ----validate-dense-spd, include=FALSE----------------------------------------
stopifnot(
  isTRUE(certificate(fit)$passed),
  isTRUE(all.equal(sort(Re(values(fit))), c(2, 4, 6)))
)

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

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

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

## ----validate-dense-general, include=FALSE------------------------------------
stopifnot(
  isTRUE(certificate(pencil)$passed),
  all(alpha_beta(pencil)$classification == "finite")
)

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

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

## ----validate-singular-pencil, include=FALSE----------------------------------
stopifnot(identical(
  sort(alpha_beta(singular)$classification),
  sort(c("finite", "infinite", "undefined"))
))

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

## ----validate-qz, include=FALSE-----------------------------------------------
stopifnot(
  length(values(qz)) == 2L,
  all(alpha_beta(qz)$classification == "finite")
)

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

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

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