Spectral and random walks

A **random walk** on a directed graph moves from the current vertex u to one of its out-neighbours chosen uniformly at random. The **transition matrix** P has P[u][v] = 1 / outdeg(u) for every edge u → v and 0 otherwise, so every row sums to 1. A vertex with no out-edges (a **dangling** vertex) is handled by letting it jump to any vertex with equal probability.

If dist is a row vector of probabilities over the vertices, one step of the walk gives dist · P. Applying this over and over (**power iteration**) moves the distribution toward a **stationary distribution** π that satisfies π = π · P: once there, one more step changes nothing.

Open the animation with the preset “Pure random walk converges”. The graph 0→1, 1→2, 2→0, 0→2 settles at 0.4, 0.2, 0.4 after a few dozen rounds. The preset “Pure walk on a 2-cycle never settles” shows why a plain walk is not enough: on 0↔1, 0↔2 the distribution alternates between uniform and 2/3, 1/6, 1/6 forever, because the graph is periodic. Damping removes this problem.

1// p[u][v] = probability that a random walker standing at u moves to v next (a vertex without out-edges jumps anywhere)
2static double[][] transitionMatrix(int n, List<Integer>[] out) {
3    double[][] p = new double[n][n];
4    for (int u = 0; u < n; u++) {
5        if (out[u].isEmpty()) { Arrays.fill(p[u], 1.0 / n); continue; }
6        for (int v : out[u]) p[u][v] += 1.0 / out[u].size();
7    }
8    return p;
9}
10static double[] step(double[] dist, double[][] p) {                   // the distribution after one step: dist * P
11    int n = dist.length;
12    double[] next = new double[n];
13    for (int u = 0; u < n; u++) for (int v = 0; v < n; v++) next[v] += dist[u] * p[u][v];
14    return next;
15}

1// out[v] = successors of v; d is the damping factor in (0, 1]; returns the rank vector after the given number of rounds (it sums to 1).
2// A vertex without out-edges is "dangling": its rank is shared equally by all vertices.
3static double[] pageRank(int n, List<Integer>[] out, double d, int rounds) {
4    double[] rank = new double[n];
5    Arrays.fill(rank, 1.0 / n);                                       // start from the uniform distribution
6    for (int round = 0; round < rounds; round++) {
7        double[] next = new double[n];
8        double dangling = 0;
9        for (int u = 0; u < n; u++) {
10            if (out[u].isEmpty()) { dangling += rank[u]; continue; }
11            for (int v : out[u]) next[v] += d * rank[u] / out[u].size();   // u hands its rank equally to its out-neighbours
12        }
13        for (int v = 0; v < n; v++) next[v] += d * dangling / n + (1 - d) / n;   // dangling mass and teleportation, spread evenly
14        rank = next;
15    }
16    return rank;
17}

**PageRank** models a surfer who, with probability d (the damping factor, classically 0.85), follows a random link from the current page and, with probability 1 − d, jumps to a page chosen uniformly at random. A page's rank is the long-run fraction of time the surfer spends on it: pages with many incoming links from important pages rank high.

The teleportation makes the chain strongly connected and aperiodic, so the stationary distribution exists, is unique and power iteration converges **geometrically at rate d per round**: after 60 rounds with d = 0.85 the error is about 0.85^60 ≈ 6 · 10^-5. Each round hands d · rank[u] / outdeg(u) to every out-neighbour and adds (1 − d) / n (plus the dangling mass) to everyone.

On the animation's default web 0→1, 0→2, 1→2, 2→0, 3→2 with d = 0.85 the ranks after 30 rounds are 0.3725, 0.1958, 0.3941, 0.0375. Page 2 is the best ranked: it has three in-links. Page 3 has no in-links at all, so its rank is exactly the teleportation share (1 − d) / n = 0.15 / 4 = 0.0375. With d = 0.5 page 3 gets 0.125. The dangling-page preset shows the leaked rank being redistributed.

1// out[v] = successors of v; d is the damping factor in (0, 1]; returns the rank vector after the given number of rounds (it sums to 1).
2// A vertex without out-edges is "dangling": its rank is shared equally by all vertices.
3static double[] pageRank(int n, List<Integer>[] out, double d, int rounds) {
4    double[] rank = new double[n];
5    Arrays.fill(rank, 1.0 / n);                                       // start from the uniform distribution
6    for (int round = 0; round < rounds; round++) {
7        double[] next = new double[n];
8        double dangling = 0;
9        for (int u = 0; u < n; u++) {
10            if (out[u].isEmpty()) { dangling += rank[u]; continue; }
11            for (int v : out[u]) next[v] += d * rank[u] / out[u].size();   // u hands its rank equally to its out-neighbours
12        }
13        for (int v = 0; v < n; v++) next[v] += d * dangling / n + (1 - d) / n;   // dangling mass and teleportation, spread evenly
14        rank = next;
15    }
16    return rank;
17}

1// out[v] = successors of v; d is the damping factor in (0, 1]; returns the rank vector after the given number of rounds (it sums to 1).
2// A vertex without out-edges is "dangling": its rank is shared equally by all vertices.
3static double[] pageRank(int n, List<Integer>[] out, double d, int rounds) {
4    double[] rank = new double[n];
5    Arrays.fill(rank, 1.0 / n);                                       // start from the uniform distribution
6    for (int round = 0; round < rounds; round++) {
7        double[] next = new double[n];
8        double dangling = 0;
9        for (int u = 0; u < n; u++) {
10            if (out[u].isEmpty()) { dangling += rank[u]; continue; }
11            for (int v : out[u]) next[v] += d * rank[u] / out[u].size();   // u hands its rank equally to its out-neighbours
12        }
13        for (int v = 0; v < n; v++) next[v] += d * dangling / n + (1 - d) / n;   // dangling mass and teleportation, spread evenly
14        rank = next;
15    }
16    return rank;
17}

The **graph Laplacian** of an undirected simple graph is L = D − A, where D is the diagonal matrix of degrees and A the adjacency matrix. Each row sums to 0, so the all-ones vector is an eigenvector with eigenvalue 0. For a path of three vertices L = [[1, −1, 0], [−1, 2, −1], [0, −1, 1]] with eigenvalues 0, 1 and 3.

The Laplacian measures smoothness: for any vector x, xᵀ L x = Σ over edges (x_u − x_v)². This makes L positive semidefinite, with all eigenvalues non-negative, and it explains the main structural fact: **the multiplicity of the eigenvalue 0 equals the number of connected components**. A vector that is constant on each component has zero energy, and nothing else does.

The code builds L and counts the zero eigenvalues as n − rank(L) by Gaussian elimination with partial pivoting. This is not faster than a BFS for counting components, but it shows the algebra at work, and the same matrix gives more: the second-smallest eigenvalue (the **algebraic connectivity**) is positive exactly when the graph is connected and is large for well-connected graphs.

1static double[][] laplacian(int n, int[][] edges) {                   // L = D - A of an undirected simple graph
2    double[][] l = new double[n][n];
3    for (int[] e : edges) {
4        l[e[0]][e[0]]++;
5        l[e[1]][e[1]]++;
6        l[e[0]][e[1]]--;
7        l[e[1]][e[0]]--;
8    }
9    return l;
10}
11static int componentsFromLaplacian(double[][] l) {                    // the number of zero eigenvalues is n - rank(L), and it equals the number of components
12    int n = l.length, rank = 0;
13    double[][] m = new double[n][];
14    for (int i = 0; i < n; i++) m[i] = l[i].clone();
15    for (int col = 0; col < n && rank < n; col++) {                   // Gaussian elimination with partial pivoting
16        int pivot = rank;
17        for (int r = rank + 1; r < n; r++) if (Math.abs(m[r][col]) > Math.abs(m[pivot][col])) pivot = r;
18        if (Math.abs(m[pivot][col]) < 1e-9) continue;
19        double[] t = m[pivot]; m[pivot] = m[rank]; m[rank] = t;
20        for (int r = rank + 1; r < n; r++) {
21            double f = m[r][col] / m[rank][col];
22            for (int c = col; c < n; c++) m[r][c] -= f * m[rank][c];
23        }
24        rank++;
25    }
26    return n - rank;
27}

**Spectral clustering** cuts a graph where it is weakly connected, using the eigenvector of the second-smallest Laplacian eigenvalue, the **Fiedler vector**. It minimises Σ (x_u − x_v)² over edges for vectors that are orthogonal to the all-ones vector and have unit length: vertices joined by many edges get similar values. The signs of the Fiedler vector therefore split the vertices into two sides with few edges between them.

To compute it without a linear-algebra library, run **power iteration** on shift · I − L, where shift is at least the largest eigenvalue. That matrix has the same eigenvectors as L with the order of the eigenvalues reversed, so the dominant eigenvector is that of the smallest eigenvalue of L: the all-ones vector. Projecting every iterate orthogonal to the all-ones vector removes it, and the iteration converges to the Fiedler vector. The code does exactly this and returns the side of each vertex.

The verification ran two cliques of 3 to 5 vertices joined by a single edge and checked that the sign split recovers the two cliques; it also checked that a path with an even number of vertices splits into its two halves. For more than two clusters, take the first k eigenvectors, treat each vertex as a point in k dimensions and run k-means. The method has no guarantee of the optimal cut, but the Cheeger inequality bounds how good the sign cut is in terms of the eigenvalue.

1// splits a connected graph in two by the sign of the Fiedler vector: the eigenvector of the second-smallest eigenvalue of L
2static boolean[] spectralBisection(int n, int[][] edges) {
3    double[][] l = laplacian(n, edges);
4    double shift = 0;
5    for (int i = 0; i < n; i++) shift = Math.max(shift, 2 * l[i][i]);  // shift >= the largest eigenvalue, so shift*I - L has the same eigenvectors, reversed in order
6    double[] x = new double[n];
7    Random random = new Random(1);
8    for (int i = 0; i < n; i++) x[i] = random.nextDouble() - 0.5;
9    for (int iter = 0; iter < 3000; iter++) {                         // power iteration on shift*I - L
10        double mean = 0;
11        for (double v : x) mean += v / n;
12        for (int i = 0; i < n; i++) x[i] -= mean;                     // stay orthogonal to the all-ones vector (eigenvalue 0 of L)
13        double[] y = new double[n];
14        for (int i = 0; i < n; i++) {
15            y[i] = shift * x[i];
16            for (int j = 0; j < n; j++) y[i] -= l[i][j] * x[j];
17        }
18        double norm = 0;
19        for (double v : y) norm += v * v;
20        norm = Math.sqrt(norm);
21        for (int i = 0; i < n; i++) x[i] = y[i] / norm;
22    }
23    boolean[] side = new boolean[n];
24    for (int i = 0; i < n; i++) side[i] = x[i] > 0;
25    return side;
26}