Graph Laplacians (in Python)

spectral.py

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

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