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

## ----setup-opts, include = FALSE----------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE,
  fig.width = 7,
  fig.height = 4.1,
  fig.align = "center",
  out.width = "92%",
  dpi = 100
)

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

## ----helpers, include = FALSE-------------------------------------------------
# Shared plotting helper. Kept out of the reader's view: the reader sees
# the API call and the picture, never the plotting boilerplate.
ec_blue <- "#2166ac"; ec_grey <- "grey78"

# A computed slice (blue) drawn on top of the full dense spectrum (grey).
plot_spectrum <- function(all_vals, computed, ylab = "value", main = NULL) {
  s <- sort(all_vals, decreasing = TRUE)
  k <- length(computed)
  op <- par(mar = c(4, 4.6, if (is.null(main)) 1 else 2.4, 1)); on.exit(par(op))
  plot(seq_along(s), s, pch = 19, cex = 0.4, col = ec_grey, bty = "n",
       xlab = "rank (largest to smallest)", ylab = ylab, main = main)
  points(seq_len(k), sort(computed, decreasing = TRUE),
         pch = 19, cex = 1.2, col = ec_blue)
  legend("topright", c("full spectrum (dense reference)", "computed by eigencore"),
         pch = 19, pt.cex = c(0.7, 1.2), col = c(ec_grey, ec_blue),
         bty = "n", cex = 0.9)
}

## ----setup--------------------------------------------------------------------
library(eigencore)
library(Matrix)

## ----simulate-signal----------------------------------------------------------
set.seed(1)
m <- 2000L; n <- 500L; k_true <- 5L
U <- qr.Q(qr(matrix(rnorm(m * k_true), m, k_true)))
V <- qr.Q(qr(matrix(rnorm(n * k_true), n, k_true)))
signal <- U %*% diag(seq(8, 4, length.out = k_true)) %*% t(V)

## ----simulate-sparsify--------------------------------------------------------
mask  <- rsparsematrix(m, n, density = 0.04)
noise <- rsparsematrix(m, n, density = 0.01, rand.x = function(k) rnorm(k, sd = 0.05))
A <- as((signal * (mask != 0)) + noise, "dgCMatrix")
dim(A)

## ----simulate-memory, echo = FALSE--------------------------------------------
dense_mb  <- as.numeric(object.size(matrix(0, m, n))) / 1024^2
sparse_mb <- as.numeric(object.size(A)) / 1024^2

## ----center-------------------------------------------------------------------
Ac <- center(A, columns = TRUE)
Ac

## ----first-svd----------------------------------------------------------------
fit <- svd_partial(Ac, rank = 5, target = largest())
fit

## ----first-scree, echo = FALSE, fig.cap = "The five largest singular values of the centered matrix (blue) against the full singular spectrum (grey).", fig.alt = "Scatter plot of all 500 singular values sorted descending in grey, with the top five highlighted in blue."----
dense_centered <- scale(as.matrix(A), center = TRUE, scale = FALSE)
all_sv <- svd(dense_centered, nu = 0, nv = 0)$d
plot_spectrum(all_sv, fit$d, ylab = "singular value")

## ----composed-certificate-----------------------------------------------------
fit$certificate$max_backward_error
fit$certificate$norm_bound_type
fit$certificate$scale_is_estimate

## ----verify-dense-------------------------------------------------------------
max(abs(sort(fit$d, decreasing = TRUE) - sort(all_sv[1:5], decreasing = TRUE)))

## ----scale-cols---------------------------------------------------------------
col_var <- (Matrix::colSums(A^2) - m * Matrix::colMeans(A)^2) / (m - 1)
w <- 1 / sqrt(col_var)
Acs <- scale_cols(Ac, w)
fit_scaled <- svd_partial(Acs, rank = 5, target = largest())
fit_scaled$d

## ----validate-scale-cols, include = FALSE-------------------------------------
stopifnot(
  all(is.finite(col_var)),
  all(col_var > 0),
  max(abs(col_var * w^2 - 1)) < 1e-12
)

## ----plan---------------------------------------------------------------------
plan <- plan_solver(svd_problem(Acs), rank = 5, target = largest())
plan$method
plan$reasons

## ----matrix-free--------------------------------------------------------------
op <- linear_operator(
  dim = dim(A),
  apply = function(X, alpha = 1, beta = 0, Y = NULL) {
    # Sparse %*% dense returns an S4 dgeMatrix; as.matrix() keeps the
    # callback's return type consistent with what the native solver expects.
    Z <- alpha * as.matrix(A %*% X)
    if (is.null(Y) || beta == 0) Z else Z + beta * Y
  },
  apply_adjoint = function(X, alpha = 1, beta = 0, Y = NULL) {
    Z <- alpha * as.matrix(crossprod(A, X))
    if (is.null(Y) || beta == 0) Z else Z + beta * Y
  },
  name = "matrix-free A"
)
fit_mf <- svd_partial(op, rank = 5, target = largest())
fit_mf$certificate$norm_bound_type

## ----matrix-free-hint---------------------------------------------------------
op_hinted <- linear_operator(
  dim = dim(A),
  apply = op$apply,
  apply_adjoint = op$apply_adjoint,
  name = "matrix-free A (exact norm)",
  metadata = list(frobenius_norm = norm(A, type = "F"))
)
fit_hinted <- svd_partial(op_hinted, rank = 5, target = largest())
fit_hinted$certificate$norm_bound_type
fit_hinted$certificate$passed