Graph Laplacians (in R)

spectral.R

# Create the degree matrix D from an affinity/adjacency matrix
degmat <- function(X) {
  diag(colSums(X))
}

# A random affinity matrix
randw <- function(n = 3) {
  X <- matrix(rnorm(n * n), nrow = n)
  X <- X * X
  X <- t(X) + X
  diag(X) <- 0
  
  list(W = X, D = degmat(X))
}

# https://stackoverflow.com/questions/16584948/how-to-create-weighted-adjacency-list-matrix-from-edge-list
# g1 graph from https://dominikschmidt.xyz/spectral-clustering-exp/
# edgelist <- list(c(1, 2), c(1, 3), c(1, 6), c(1, 9), c(1, 10), c(2, 3), 
#                  c(9, 10), c(6, 4), c(6, 5), c(6, 7), c(6, 8), c(4, 5), 
#                  c(7, 8))
# G <- edge2adj(edgelist)
edge2adj <- function(edgelist, n = max(sapply(edgelist, max))) {
  G <- matrix(0, n, n)
  edges <- cbind(a = sapply(edgelist, `[[`, 1), 
                 b = sapply(edgelist, `[[`, 2))
  G[edges] <- 1
  G + t(G)
}

# The equivalent of randw, but using an edgelist
edge2graph <- function(edgelist, n = max(sapply(edgelist, max))) {
  G <- edge2adj(edgelist, n = n)
  list(W = G, D = degmat(G))
}

# Generate various Laplacians
lapm <- function(WD) {
  W <- WD$W
  D <- WD$D

  L <- D - W
  
  # Commented out code more closely follows the mathematical definitions, but
  # pointlessly stores and inverts the full diagonal matrix as well as carrying
  # out matrix multiplications
  # Dinv <- solve(D)
  Dinv <- 1 / diag(D)
  
  # Lsym <- Dinvs %*% L %*% Dinvs
  Dinvs <- sqrt(Dinv)
  Lsym <- Dinvs * sweep(L, 2, Dinvs, '*')

  #P <- Dinv %*% W
  P <- Dinv * W
  
  # I <- diag(nrow = nrow(D), ncol = ncol(D))
  # Lrw <- I - P
  Lrw <- -P
  diag(Lrw) <- 1 - diag(Lrw)
  
  list(L = L, Lsym = Lsym, Lrw = Lrw, P = P)
}

# Calculate eigenvectors/values
eig <- function(X, norm = FALSE, val1 = FALSE) {
  res <- eigen(X)
  sorteig(res, norm = norm, val1 = val1)
}

# Calculate generalized eigenvectors/values
geig <- function(A, B, norm = FALSE, val1 = FALSE) {
  res <- geigen::geigen(A, B)
  sorteig(res, norm = norm, val1 = val1)
}

# Calculate k smallest eigenvectors/values
reig <- function(X, k, norm = FALSE, val1 = FALSE) {
  res <- RSpectra::eigs(X, k = k, which = "SM")
  res$vectors <- Re(res$vectors)
  res$values <- Re(res$value)
  sorteig(res, norm = norm, val1 = val1)
}

# Sort eigenvectors by value
sorteig <- function(X, norm = FALSE, val1 = FALSE) {
  vectors <- X$vectors
  values <- X$values

  if (val1) {
    values <- 1 - values
  }
  
  if ((is.logical(norm) && norm) || is.numeric(norm) || (is.character(norm) && norm == "n")) {
    if (is.logical(norm)) {
      m <- 1
    }
    else if (is.numeric(norm)) {
      m <- norm
    }
    else {
      m <- sqrt(nrow(vectors)) # make smallest eigenvector all 1s
    }
    sqrtcsums <- sqrt(colSums(vectors * vectors))
    vectors <- m * sweep(vectors, 2, sqrtcsums, "/")
  }

  vectors <- vectors[, order(values)]
  values <- sort(values)

  list(vectors = vectors, values = values, lengths = sqrt(colSums(vectors ^ 2)))
}

Code for exploring basic Graph Laplacian properties, as discussed at https://jlmelville.github.io/smallvis/spectral.html. Python code is at https://gist.github.com/jlmelville/8b7655fb4803ce49e4f560d316b04a46.

• randw: generate a list containing: W, a typical random affinity matrix: positive semi definite, symmetric, with zeros on the diagonal, and D, the degree matrix.

• edge2graph: take a list of edges (as two-element vectors with vertices numbered from 1) and generate the equivalent of randw, where W is the adjacency matrix. Edges all have a weight of one.

• lapm: takes the list of matrices returned by randw (or edge2graph) and returns its own list containing the various Laplacian matrices: L, the un-normalized Laplacian; Lsym, the symmetrized normalized Laplacian; Lrw, the random walk Laplacian; and P the random walk transition matrix.

• eig: calculates the eigenvectors and eigenvalues of an input matrix X, and returns them in increasing order. If you set norm = TRUE, it will return the eigenvectors after normalizing their length to 1. If you set norm to a numeric value, it will return the eigenvectors scaled such that their length is that value. To get the lowest eigenvector to be all 1s, you should set norm to the square root of the number of rows in X. As a short-cut you can set norm = "n". If val1 = TRUE then the eigenvalues are subtracted from 1 before returning. This is useful for the comparison of the eigenvalues of P and Lrw (and hence the connection between Laplacian Eigenmaps and Diffusion Maps). The return value is the vectors as a matrix with eigenvectors in each column, and values, with the corresponding eigenvalues. Also, lengths, which gives the length of each vector (after any scaling by norm).

• geig: calculates the generalized eigenvectors and eigenvalues for A v = lambda B v. Has the same parameters and return value structure as eig. You need to install geigen for this to work.

• reig: uses the RSpectra package to only calculate the first k eigenvectors, so you must provide k. If you set k to be the full matrix it will just use eigen anyway.

Examples:

WD <- randw(5)
lap <- lapm(WD)

# Laplacian Eigenmaps: these two should give the same results
# use norm = "n", because otherwise the eigenvectors can have different lengths
geig(lap$L, WD$D, norm = "n")
eig(lap$Lrw, norm = "n")

# g1 graph from https://dominikschmidt.xyz/spectral-clustering-exp/
edgelist <- list(c(1, 2), c(1, 3), c(1, 6), c(1, 9), c(1, 10), c(2, 3), 
                 c(9, 10), c(6, 4), c(6, 5), c(6, 7), c(6, 8), c(4, 5), 
                 c(7, 8))
G <- edge2graph(edgelist)
lapG <- lapm(G)

# compare with results at https://dominikschmidt.xyz/spectral-clustering-exp/
eig(lapG$L)