Instantly share code, notes, and snippets.

# jlmelville/spectral.py

Created March 14, 2020 18:00
Show Gist options
• Save jlmelville/8b7655fb4803ce49e4f560d316b04a46 to your computer and use it in GitHub Desktop.
Graph Laplacians (in Python)
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
 import collections import pprint import numpy as np import scipy as sp import scipy.linalg AffDeg = collections.namedtuple("AffDeg", ["W", "D"]) Lap = collections.namedtuple("Lap", ["L", "Lsym", "Lrw", "P"]) Eig = collections.namedtuple("Eig", ["vectors", "values", "lengths"]) def colsums(X): return np.sum(X, axis=0) def degmat(X): return np.diag(colsums(X)) def randw(n=3): X = np.random.randn(n, n) X = X * X X = X.transpose() + X np.fill_diagonal(X, 0) return AffDeg(X, degmat(X)) def lapm(aff_deg): W = aff_deg.W D = aff_deg.D L = D - W Dinv = np.diag(1 / np.diag(D)) Dinvs = np.sqrt(Dinv) Lsym = np.matmul(Dinvs, np.matmul(L, Dinvs)) P = np.matmul(Dinv, W) Lrw = -P np.fill_diagonal(Lrw, 1 - np.diag(Lrw)) return Lap(L, Lsym, Lrw, P) def eig(X, norm=False, val1=False): res = np.linalg.eig(X) return sorteig(res[1], res[0], norm=norm, val1=val1) def geig(A, B, norm=False, val1=False): res = sp.linalg.eigh(A, B, eigvals_only=False) return sorteig(res[1], res[0], norm=norm, val1=val1) def sorteig(vectors, values, norm=False, val1=False): if val1: values = 1 - values if ( (isinstance(norm, bool) and norm) or isinstance(norm, float) or (isinstance(norm, str) and norm == "n") ): if isinstance(norm, bool): m = 1 elif isinstance(norm, float): m = norm else: # norm = "n" make smallest eigenvector all 1s m = np.sqrt(vectors.shape[0]) sqrtcsums = np.sqrt(colsums(vectors * vectors)) vectors = m * vectors / sqrtcsums vectors = vectors.transpose()[np.argsort(values)].transpose() values = np.sort(values) return Eig(vectors, values, np.sqrt(colsums(vectors * vectors)))

### jlmelville commented Mar 14, 2020 • edited

Python code to accompany https://jlmelville.github.io/smallvis/spectral.html. The R equivalent is at https://gist.github.com/jlmelville/772060a26001d7d25d7453b0df4feff9.

```    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
print(geig(lap.L, WD.D, norm="n"))
print(eig(lap.Lrw, norm="n"))

# Using P also gives the same results, if the eigenvalues are transformed by 1 - lambda
print(eig(lap.P, norm="n", val1=True))

# eigenvectors are stored by row, i.e. eig(lap.Lrw, norm="n").vectors[:, 0] should be all 1s```

to join this conversation on GitHub. Already have an account? Sign in to comment