## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7,
                      fig.height = 4.2, dpi = 96, out.width = "100%")

## -----------------------------------------------------------------------------
library(matsketch)

## ----data---------------------------------------------------------------------
set.seed(1)
X <- rbind(
  matrix(rnorm(2 * 1300, sd = 0.5), ncol = 2),     # large cluster
  matrix(rnorm(2 * 180, sd = 0.2), ncol = 2) + 4,  # small cluster
  matrix(runif(2 * 20, -6, 10), ncol = 2)          # outliers
)
K <- kernel_matrix(X, "gaussian", bandwidth = 0.5)
K

## ----fits---------------------------------------------------------------------
rules <- c(rpchol = "accelerated", greedy = "greedy", uniform = "uniform")
fits <- lapply(rules, function(r) rpchol(K, k = 100, method = r))
ev <- eigen(as.matrix(K), symmetric = TRUE, only.values = TRUE)$values
best <- 1 - cumsum(ev[1:100]) / sum(ev)

ranks <- c(20, 40, 100)
tab <- cbind(sapply(fits, function(f) f$trace_path[ranks]), best = best[ranks])
rownames(tab) <- paste("rank", ranks)
signif(tab, 2)

## ----error-curves-------------------------------------------------------------
cols <- c(rpchol = "#D55E00", greedy = "#0072B2", uniform = "#E69F00")
plot(best, log = "y", type = "l", lwd = 2, lty = 2, col = "grey30",
     ylim = range(best, unlist(lapply(fits, `[[`, "trace_path"))),
     xlab = "rank", ylab = "relative trace error")
for (r in names(fits)) lines(fits[[r]]$trace_path, lwd = 2, col = cols[[r]])
legend("topright", c(names(fits), "best possible"), col = c(cols, "grey30"),
       lty = c(1, 1, 1, 2), lwd = 2, bty = "n")

## ----pivots, fig.height = 3---------------------------------------------------
op <- par(mfrow = c(1, 3), mar = c(0.5, 0.5, 2, 0.5))
for (r in names(fits)) {
  plot(X, pch = 16, cex = 0.35, col = "grey75", asp = 1, axes = FALSE,
       xlab = "", ylab = "", main = r)
  points(X[fits[[r]]$pivots[1:40], ], pch = 16, cex = 0.7, col = "#D55E00")
  box(col = "grey85")
}
par(op)

## ----entries------------------------------------------------------------------
sapply(fits, function(f) f$entries / K$n^2)

## ----trace-setup--------------------------------------------------------------
set.seed(2)
n <- 500
U <- qr.Q(qr(matrix(rnorm(n * n), n)))
A <- U %*% ((1:n)^-2 * t(U))
trA <- sum(diag(A))
trA
trace_est(A, m = 60)

## ----trace-curves-------------------------------------------------------------
budgets <- c(12, 24, 48, 96, 192)
methods <- c(hutchinson = "#999999", hutchpp = "#E69F00",
             xtrace = "#D55E00", xnystrace = "#0072B2")
err <- sapply(names(methods), function(meth) {
  sapply(budgets, function(m) {
    median(replicate(10, abs(trace_est(A, m, meth)$estimate - trA) / trA))
  })
})
matplot(budgets, err, log = "xy", type = "b", pch = 19, lty = 1, lwd = 2,
        col = methods, xlab = "matrix-vector products",
        ylab = "relative error")
legend("bottomleft", names(methods), col = methods, lwd = 2, pch = 19,
       bty = "n")

## ----diag---------------------------------------------------------------------
cor(diag_est(A, m = 100), diag(A))

## ----pcg----------------------------------------------------------------------
set.seed(3)
Z <- matrix(rnorm(3 * 1000), ncol = 3)
K2 <- as.matrix(kernel_matrix(Z, "gaussian", bandwidth = 2))
b <- rnorm(1000)
mu <- 1e-4

ed <- effective_dim(rpchol(K2, k = 200), mu)
ed

## ----pcg-runs-----------------------------------------------------------------
precond_of_rank <- function(k) nystrom_precond(rpchol(K2, k), mu)
runs <- list(
  plain = pcg(K2, b, mu = mu, maxit = 5000),
  `rank 50` = pcg(K2, b, mu = mu, precond = precond_of_rank(50)),
  recommended = pcg(K2, b, mu = mu,
                    precond = precond_of_rank(ed$recommended_rank))
)
sapply(runs, `[[`, "iterations")

## ----pcg-plot-----------------------------------------------------------------
cols <- c("#999999", "#E69F00", "#D55E00")
plot(runs$plain$residuals, log = "y", type = "l", lwd = 2, col = cols[1],
     xlab = "iteration", ylab = "relative residual")
for (i in 2:3) {
  lines(runs[[i]]$residuals, type = "o", pch = 19, cex = 0.4, lwd = 2,
        col = cols[i])
}
legend("topright", names(runs), col = cols, lwd = 2, bty = "n")

## ----reml---------------------------------------------------------------------
set.seed(4)
dat <- sim_genomic(n = 1000, p = 2000, h2 = 0.5, pops = 4)
G <- grm_matrix(dat$M)
G
fit <- reml_sketch(dat$y, G)
fit
exact <- reml_exact(dat$y, as.matrix(G))
rbind(sketched = c(fit$sigma2, h2 = fit$h2),
      exact = c(exact$sigma2, h2 = exact$h2))