Instantly share code, notes, and snippets.

# jlmelville/spectral.R

Last active March 14, 2020 18:06
Star You must be signed in to star a gist
Graph Laplacians (in R)
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters. Learn more about bidirectional Unicode characters
 # 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))) }

### jlmelville commented Feb 23, 2020 • edited

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