A trick to compute the svd of the population regularized Laplacian with a single extra line of code

[alexpghayes]

library(fastRG)
#> Loading required package: Matrix

# construct any undirected factor model

k <- 5
n <- 1000
B <- matrix(stats::runif(k * k), nrow = k, ncol = k)

theta <- round(stats::rlnorm(n, 2))

pi <- c(1, 2, 4, 1, 1)

mod <- dcsbm(
  theta = theta,
  B = B,
  pi = pi,
  expected_degree = 50
)

# now suppose we want the population svds of E[A], E[L], and E[L_tau]

# we already have a method the gives us the svds of E[A] with high computational
# efficiency, using the low rank structure in E[A]
s_EA <- svds(mod)

# that is, we already
# have code that computes the svds of E[A] = XSX'. we know that
# L_tau = D_tau A D_tau, and thus we have (roughly)
# E[L_tau] = E[D_tau] E[A] E[D_tau] = E[D_tau] XSX' E[D_tau]
# so we're gonna construct a new object with X_new = E[D_tau] X and run
# the old svd code on it

# we don't actually need to form the pop_D_tau matrix explicitly
# because rowScale does this for us, but here's what that would look like

# pop_degs <- expected_degrees(mod)
# pop_D_tau <- Diagonal(n = mod$n, x = 1 / sqrt(pop_degs + tau))
# X_new <- rowScale(mod$X, 1 / sqrt(pop_degs + tau))

# the very concise version of this
tau <- 10
mod$X <- rowScale(mod$X, 1 / sqrt(expected_degrees(mod) + tau))

# getting population eigenvalues of E[L_tau] only takes a single additional line
# of code
s_EL <- svds(mod)

# NOTE! Do not call `sample_sparse(mod)` anymore, it won't sample from the model
# you expect

# quick check that the SVDs are different
s_EL$d - s_EA$d
#> [1] -114.024473  -27.336110  -10.149304   -5.945424   -4.218959

^{Created on 2023-07-21 with reprex v2.0.2}