I was reading some notes from QUT about using spectral graph theory to partition a graph using eigenvalues of unnormalised Laplacian. Their Matlab code looks like this:

% Form W = Gaussian distribution based on distances W = zeros(100, 100); D = zeros(100, 100); sigma = 2; N = 100; for i = 1:N for j = 1:N if (j ~= i) % Calculate the distance between two points dist = norm([A(i, 1) - A(j, 1); A(i, 2) - A(j, 2)]); expp = exp(-dist^2 / (2 * sigma^2)); adjacency(i, j) = 1; % Add the weights to the matrix W(i, j) = expp; % Add up the row sum as we go D(i, i) = D(i, i) + expp; end end end L = D - W; Find the eigenpairs of L [vec, val] = eig(L, 'vector'); % Ensure the eigenvalues and eigenvectors are sorted in ascending order [val,ind] = sort(val); vec = vec(:, ind); % Plot with clusters identified by marker % v2 figure; hold on; for i = 1:N plot(i, vec(i, 2), ['k' plot_syms{i}]); end

Personally, this gives me nightmares from undergrad numerical methods classes. So here’s how to do it in Haskell. Full source (including cabal config) for this post is here.

> module Notes where

To create the W matrix we use buildMatrix:

> buildW :: Double -> Matrix Double -> Matrix Double > buildW sigma a = buildMatrix n n $ (i,j) -> f sigma a i j > where > n = rows a > > f sigma a i j = if j /= i > then expp sigma a i j > else 0.0 > > dist :: Matrix Double -> Int -> Int -> Double > dist m i j = norm2 $ fromList [ m!i!0 - m!j!0, m!i!1 - m!j!1 ] > > expp :: Double -> Matrix Double -> Int -> Int -> Double > expp sigma m i j = exp $ (-(dist m i j)**2)/(2*sigma**2)

The D matrix is a diagonal matrix with each element being the sum of a row of W, which is nicely expressible by composing a few functions:

> buildD :: Matrix Double -> Matrix Double > buildD w = diag > . fromList > . map (foldVector (+) 0) > . toRows > $ w > where > n = rows w

The L matrix is real and symmetric so we use eigSH which provides the eigenvalues in descending order.

> lapEigs :: Double -> Matrix Double -> (Vector Double, Matrix Double) > lapEigs sigma m = eigSH l > where > w = buildW sigma m > d = buildD w > l = d - w

To finish up, the Fiedler eigenvector corresponds to the second smallest eigenvalue:

> fiedler :: Double -> Matrix Double -> (Double, Vector Double) > fiedler sigma m = (val ! (n-2), vector $ concat $ toLists $ vec ¿ [n-2]) > where > (val, vec) = lapEigs sigma m > n = rows m

To plot the eigenvalues and eigenvector we use the Chart library which uses the cairo backend.

> doPlot "eigenvalues.png" "Eigenvalues" "eigenvalue" $ zip [0..] (reverse $ toList val) > > doPlot "fiedler.png" > "Second eigenvalue of unnormalised Laplacian" > "fiedler eigenvector" > (zip [0..] $ toList algConnecEigVec)