The graph Laplacian: why its second eigenvector cuts a network in two (Python)

The question

Zachary's karate club is a network of 34 members and 78 friendships, recorded in the early 1970s while a dispute between the instructor, Mr. Hi, and the club's president split the club in two. Seventeen members went each way. NetworkX ships the friendships together with the side each member joined, and here they are in the spring layout of NetworkX from the ground up, where linked members sit close together.

Show code

import math

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from scipy.linalg import eigh

plt.rcParams.update({
    "figure.figsize": (7, 3.6), "figure.dpi": 110,
    "axes.spines.top": False, "axes.spines.right": False,
    "axes.grid": True, "grid.alpha": 0.25,
    "font.size": 11, "lines.linewidth": 1.8,
})
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"
LAYOUT_SEED = 3

G = nx.karate_club_graph()
n = G.number_of_nodes()
side = {k: G.nodes[k]["club"] for k in G}
pos = nx.spring_layout(G, seed=LAYOUT_SEED, weight=None)   # the layout of py-networkx
SHAPES = [("Mr. Hi", "o"), ("Officer", "s")]

fig, ax = plt.subplots(figsize=(7, 5.5))
nx.draw_networkx_edges(G, pos, ax=ax, edge_color=MUTED, alpha=0.5, width=1)
for club, shape in SHAPES:
    nx.draw_networkx_nodes(G, pos, nodelist=[k for k in G if side[k] == club], node_shape=shape,
                           node_color=INK, node_size=80, edgecolors="white", ax=ax)
nx.draw_networkx_labels(G, {k: pos[k] + (0, 0.08) for k in (0, 33)}, labels={0: 0, 33: 33},
                        ax=ax, font_color=INK, bbox=dict(fc="white", ec="none", pad=1))
handles = [plt.Line2D([], [], ls="", marker=m, ms=8, color=INK, label=f"joined {name}")
           for name, m in [("Mr. Hi", "o"), ("the officer", "s")]]
ax.legend(handles=handles, frameon=False, loc="lower left")
ax.set_axis_off()
plt.show()
print(f"{n} members, {G.number_of_edges()} friendships; "
      f"node 0 has {G.degree(0)} friends, node 33 has {G.degree(33)}")

Zachary's karate club in a spring layout. Circles joined Mr. Hi, squares the officer. Two dense groups form around node 0 and node 33, with a few members in between that have friends in both.

34 members, 78 friendships; node 0 has 16 friends, node 33 has 17

Two things are obvious. The friendships crowd around two hubs, node 0 (Mr. Hi, 16 friends) and node 33 (the officer, 17), and a few members sit in between with friends on both sides. NetworkX counts members from 0 and Zachary from 1, so node 8 here is his member 9. In the NetworkX tutorial, Louvain community detection found the two halves from the friendships alone and put two of the 34 on the wrong side.

Less obvious is this. Write the friendships into a 34 × 34 matrix, one small change away from the adjacency matrix, take its eigenvectors with eigh, and look at the signs of the second one. The positive entries fall on Mr. Hi's side and the negative ones on the officer's, again with two of 34 wrong. An eigenvector is a statement about how a matrix stretches vectors, and the matrix knows nothing about karate. Why should it know where a network is thin?

The idea: the network as masses on springs

Make every member a unit mass and every friendship a unit spring, and let each mass move only up and down, by \(x_i\). The spring to a friend \(j\) pulls mass \(i\) with the force \(x_j - x_i\), so the total force on \(i\) is \(\sum_j A_{ij}(x_j - x_i) = -\bigl(d_i x_i - \sum_j A_{ij}x_j\bigr)\), with \(A_{ij} = 1\) for a friendship and 0 otherwise, and \(d_i\) the number of friends of \(i\), its degree. Newton's law for all 34 masses at once is

\[\ddot{\mathbf x} = -L\,\mathbf x ,\]

where \(L\) holds the degrees on its diagonal and \(-1\) at every friendship off it. These are the gliders of Eigenvalues with numpy.linalg, \(m\ddot{\mathbf x} = -K\mathbf x\) with \(m = 1\), only with the springs where the friendships are. So each eigenvector of \(L\) is a normal mode, a pattern in which every member swings at one frequency, and its eigenvalue is \(\lambda = \omega^2\).

Show code

A = nx.to_numpy_array(G, weight=None)     # 1 for a friendship; Zachary's counts are ignored
d = A.sum(axis=1)
L = np.diag(d) - A
w, V = eigh(L)
V = V * np.sign(V[0])                     # eigh picks signs at random; orient every mode by node 0
w = np.clip(w, 0, None)                   # the zero comes out as round-off, near -1e-16
fiedler = V[:, 1]
print(f"trace of L = {np.trace(L):.0f} = 2 x {G.number_of_edges()} friendships")
print("lowest eigenvalues:", "  ".join(f"{x:.4f}" for x in w[:4]), f"   largest: {w[-1]:.2f}")

trace of L = 156 = 2 x 78 friendships
lowest eigenvalues: 0.0000  0.4685  0.9092  1.1250    largest: 18.14

Nobody is tied to a post. Lift the whole club by the same amount and no spring stretches, as the free CO₂ molecule of that tutorial slides along its axis at no cost: a mode with \(\lambda = 0\) and every entry equal. A network in two pieces could lift each piece on its own and would have two such modes. The club is one piece, and eigh returns one zero, then 0.4685 and 0.9092, and 18.14 at the top. The first mode that moves is the second one. Here are the first, the second, and the last, each member lifted by its entry in the eigenvector above a gray rest mark:

Show code

xy = np.array([pos[k] for k in G])
edges = np.array(G.edges)

def draw_mode(ax, k, lift=0.35):
    """The club with every member displaced by its entry of mode k, largest entry = lift."""
    shape = V[:, k] / np.abs(V[:, k]).max() * lift
    moved = xy + np.c_[np.zeros(n), shape]
    ax.plot(np.c_[xy[:, 0] - 0.05, xy[:, 0] + 0.05].T, np.c_[xy[:, 1], xy[:, 1]].T, color=MUTED, lw=1)
    for i, j in edges:
        ax.plot(moved[[i, j], 0], moved[[i, j], 1], color=MUTED, lw=0.6, alpha=0.6)
    ax.scatter(moved[:, 0], moved[:, 1], s=30, zorder=3, edgecolors="white", linewidths=0.6,
               c=[SECOND if s >= 0 else ACCENT for s in V[:, k]])
    lam = w[k]
    ax.text(0.0, -0.02, f"mode {k + 1}\nλ = {lam:.4g}, ω = {np.sqrt(lam):.3g}", transform=ax.transAxes,
            va="top", color=INK)
    ax.set_axis_off()
    return moved

fig, axes = plt.subplots(1, 3, figsize=(8, 3.3), sharex=True, sharey=True)
for ax, k in zip(axes, [0, 1, n - 1]):
    moved = draw_mode(ax, k)
hub = np.abs(V[:, -1]).argmax()           # the member that dominates the last mode
axes[2].annotate(f"{hub}", moved[hub], xytext=(5, 4), textcoords="offset points", color=INK)
plt.show()

stretch = np.abs(fiedler[edges[:, 0]] - fiedler[edges[:, 1]])
across = (fiedler[edges[:, 0]] > 0) != (fiedler[edges[:, 1]] > 0)   # the springs between the two signs
print(f"mode 2: {across.sum()} of {len(edges)} springs across, "
      f"{100 * (stretch[across]**2).sum() / (stretch**2).sum():.0f} % of the spring energy, "
      f"stretched {stretch[across].mean() / stretch[~across].mean():.1f} x as far as one inside")

Three normal modes of the karate club, members displaced up or down from gray rest marks. Mode 1, eigenvalue 0: everyone lifted together. Mode 2, eigenvalue 0.4685: one half up, the other down. Mode 34, eigenvalue 18.14: node 33 against its friends.

mode 2: 10 of 78 springs across, 33 % of the spring energy, stretched 2.3 x as far as one inside

In the second mode, at \(\lambda = 0.4685\) and \(\omega = 0.684\), the two halves of the club swing against each other, one up and one down. The ten springs between the two signs stretch on average 2.3 times as far as one inside a half, yet they hold only a third of the energy, because the 68 inside outnumber them almost seven to one. The last mode is the other extreme: the hub, node 33, moves against all 17 of its friends at once, the stiffest motion the network has. Watch the three swing in turn, each slowed to the same 2.4 s:

The karate club as masses on springs, swinging in three normal modes in turn. Mode 1: everyone moves up and down together. Mode 2: the two halves move against each other, the springs across the middle stretched. Mode 34: node 33 alone against its friends.

Show code

"""Club modes: three normal modes of Zachary's karate club as masses on springs.

Every member is a unit mass, every friendship a unit spring, and each mass moves only
up and down. The modes are the eigenvectors of the graph Laplacian L = D - A. Each of
three modes swings for one cycle, slowed to the same 2.4 s. Renders ../../assets/club-modes.gif:

    python scene.py
"""
from pathlib import Path

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.collections import LineCollection
from PIL import Image
from scipy.linalg import eigh

OUT = Path(__file__).resolve().parents[2] / "assets" / "club-modes.gif"
INK, ACCENT, SECOND, MUTED = "#1f2a44", "#c8553d", "#2a7f9e", "#8a8f98"
plt.rcParams.update({"font.size": 13})   # the title stays legible when the GIF shrinks to a phone

# ---- the club, its Laplacian, and three of its modes
G = nx.karate_club_graph()
pos = nx.spring_layout(G, seed=3, weight=None)
xy = np.array([pos[n] for n in G])
A = nx.to_numpy_array(G, weight=None)
w, V = eigh(np.diag(A.sum(axis=1)) - A)
V = V * np.sign(V[0])                    # eigh picks signs at random; node 0 (Mr. Hi) moves up
edges = np.array(G.edges)

MODES = [0, 1, len(w) - 1]               # the zero mode, the Fiedler mode, the stiffest mode
NAMES = ["mode 1: everyone together", "mode 2: half against half", f"mode {len(w)}: the hub against its friends"]
FRAMES_PER_MODE = 30                     # 2.4 s: a GIF stores 1/12 s as 80 ms
LIFT = 0.35                              # largest displacement, in layout units

# ---- figure, drawn once
fig, ax = plt.subplots(figsize=(7, 4.6), dpi=80, layout="constrained")
ax.set_axis_off()
ax.set(xlim=(xy[:, 0].min() - 0.1, xy[:, 0].max() + 0.1),
       ylim=(xy[:, 1].min() - LIFT - 0.05, xy[:, 1].max() + LIFT + 0.05))
ax.plot(np.c_[xy[:, 0] - 0.04, xy[:, 0] + 0.04].T, np.c_[xy[:, 1], xy[:, 1]].T,
        color=MUTED, lw=1)                # a rest mark under every member
springs = LineCollection([], colors="#b9bcc2", linewidths=1)   # MUTED, lightened: no alpha keeps the GIF small
ax.add_collection(springs)
masses = ax.scatter(xy[:, 0], xy[:, 1], s=60, zorder=3, edgecolors="white", linewidths=1)
title = ax.set_title("", loc="left", color=INK)

# ---- one frame: mode k, displaced by its shape times sin of the phase
def update(frame):
    i, step = divmod(frame, FRAMES_PER_MODE)
    k = MODES[i]
    shape = V[:, k] / np.abs(V[:, k]).max() * LIFT
    moved = xy + np.c_[np.zeros(len(xy)), shape * np.sin(2 * np.pi * step / FRAMES_PER_MODE)]
    springs.set_segments(moved[edges])
    masses.set_offsets(moved)
    masses.set_facecolor([SECOND if s >= 0 else ACCENT for s in V[:, k]])
    lam = w[k] if w[k] > 1e-9 else 0.0    # the zero mode comes out as round-off, about 1e-16
    title.set_text(f"{NAMES[i]}   λ = {lam:.4g}, ω = {np.sqrt(lam):.3g}")

# ---- render, and read back what was written
OUT.parent.mkdir(exist_ok=True)
FuncAnimation(fig, update, frames=len(MODES) * FRAMES_PER_MODE).save(OUT, writer=PillowWriter(fps=12))
plt.close(fig)
with Image.open(OUT) as im:              # Pillow merges repeated frames, so add up the stored delays
    delays = []
    for i in range(im.n_frames):
        im.seek(i)
        delays.append(im.info["duration"])
    print(f"{OUT.name}: {im.width} x {im.height} px, {im.n_frames} frames, "
          f"{sum(delays) / 1000:.1f} s, {OUT.stat().st_size / 1024:,.0f} kB")

Counting a cut with plus and minus one

The energy in the springs is half the sum of the squared stretches, which in matrix form reads

\[\tfrac12\,\mathbf x^{\mathsf T} L\,\mathbf x = \tfrac12 \sum_{\text{friendships } ij} (x_i - x_j)^2 .\]

For a normal mode \(L\mathbf x = \lambda\mathbf x\), so \(\lambda = \mathbf x^{\mathsf T}L\mathbf x / \mathbf x^{\mathsf T}\mathbf x\), the spring energy per squared amplitude. A low eigenvalue is a shape that stretches few springs for its size.

Now give the shape two values only, \(x_i = +1\) on one side of a split and \(-1\) on the other. A friendship inside a side stretches by 0. One across the split stretches by 2 and adds 4. So \(\mathbf x^{\mathsf T}L\mathbf x\) is four times the cut, the number of friendships the split crosses, and \(\mathbf x^{\mathsf T}\mathbf x = n\) for every such vector. For the split by the signs of the second eigenvector:

Show code

x = np.where(fiedler > 0, 1.0, -1.0)
hi = [k for k in G if fiedler[k] > 0]
print(f"x L x = {x @ L @ x:.0f}   4 x cut = 4 x {nx.cut_size(G, hi)}   x x = {x @ x:.0f}")
print(f"balanced +-1 vectors: {math.comb(n, n // 2):,}")

x L x = 40   4 x cut = 4 x 10   x x = 34
balanced +-1 vectors: 2,333,606,220

Forty, four times the ten friendships across the sign split. A good split is a \(\pm 1\) vector with a small \(\mathbf x^{\mathsf T}L\mathbf x\).

Minimized over all vectors, \(\mathbf x^{\mathsf T}L\mathbf x\) finds nothing: the zero vector wins, and so does the constant vector, which puts everyone on one side and cuts nothing. Two constraints rule them out. Fix the length, \(\mathbf x^{\mathsf T}\mathbf x = n\), as every \(\pm 1\) vector has it. And ask that the entries sum to zero, \(\mathbf x \perp \mathbf 1\), with \(\mathbf 1\) the vector of all ones and \(\perp\) meaning perpendicular, as the columns eigh returns are to each other. For a \(\pm 1\) vector a zero sum means equal head counts. For real entries it only balances the plus against the minus.

With \(\pm 1\) entries the problem is a search through 2.3 billion balanced vectors, every split into 17 and 17 counted once for each sign. Let the entries be any real numbers under the same two constraints, and it asks only for the smallest \(\mathbf x^{\mathsf T}L\mathbf x / \mathbf x^{\mathsf T}\mathbf x\) over \(\mathbf x \perp \mathbf 1\). Loosening a hard problem this way is a relaxation.

The answer is \(\lambda_2\), and the modes say why. Write \(\mathbf x = \sum_k c_k\mathbf v_k\) in the orthonormal eigenvectors. Then \(\mathbf x^{\mathsf T}L\mathbf x = \sum_k \lambda_k c_k^2\) and \(\mathbf x^{\mathsf T}\mathbf x = \sum_k c_k^2\), so the quotient is an average of the eigenvalues, each weighted by its share \(c_k^2\). In a connected network the only zero mode is \(\mathbf v_1\), every entry \(1/\sqrt n\), and \(\mathbf x \perp \mathbf 1\) means \(c_1 = 0\). The average then runs over \(\lambda_2\) and up, an average never falls below its smallest term, and it reaches \(\lambda_2\) when all the weight sits on mode 2. Here is the members' own split, as a \(\pm 1\) vector:

Show code

y = np.array([1.0 if side[k] == "Mr. Hi" else -1.0 for k in G])
c = V.T @ y                               # coefficients in the orthonormal modes
share = c**2 / (c**2).sum()
print(f"sum of entries = {y.sum():.0f}   y L y = {y @ L @ y:.0f} = 4 x {nx.cut_size(G, [k for k in G if y[k] > 0])}")
print(f"quotient = {y @ L @ y / (y @ y):.3f}   weighted average of eigenvalues = {(share * w).sum():.3f}")
print("share on modes 1 to 4:", "  ".join(f"{s:.2f}" for s in share[:4]))

sum of entries = 0   y L y = 44 = 4 x 11
quotient = 1.294   weighted average of eigenvalues = 1.294
share on modes 1 to 4: 0.00  0.66  0.07  0.11

The quotient is 1.294, the weighted average of the eigenvalues and above 0.4685, as it has to be. Of its weight, 66 % sits on mode 2 and none on mode 1: the members' own line is mostly the second mode. That eigenvector, the answer of the relaxation, is the Fiedler vector. Turning its real entries back into two sides is called rounding.

Formalization

The matrix of the springs is the graph Laplacian,

\[L = D - A ,\]

with \(A\) the adjacency matrix and \(D\) the diagonal matrix of the degrees. With weights on the edges, such as confidence scores, \(A\) becomes the weight matrix \(W\), \(D\) holds each node's total weight, its strength, and everything below holds unchanged. The energy identity, \(\mathbf x^{\mathsf T}L\mathbf x = \sum_{ij} w_{ij}(x_i - x_j)^2\) over the edges, is a sum of squares and never negative, which is what positive semidefinite means. It is zero only for an \(\mathbf x\) that is constant on every connected piece, so the number of zero eigenvalues is the number of pieces. Delete the ten friendships of the sign split:

Show code

W = nx.to_numpy_array(G)                  # Zachary's counts as weights, for this one check
print("laplacian_matrix(G) equals diag(W 1) - W:", np.allclose(nx.laplacian_matrix(G).toarray(), np.diag(W.sum(1)) - W))

H = G.copy()
H.remove_edges_from([(i, j) for i, j in G.edges if x[i] != x[j]])
wH = np.linalg.eigvalsh(nx.laplacian_matrix(H, weight=None).toarray())
print(f"zero eigenvalues: {(np.abs(wH) < 1e-9).sum()}   next: {wH[2]:.3f}   "
      f"components: {nx.number_connected_components(H)}   "
      f"algebraic_connectivity: {nx.algebraic_connectivity(H, weight=None):.1f}")
try:
    nx.fiedler_vector(H, weight=None, seed=0)
except nx.NetworkXError as err:
    print("fiedler_vector:", err)

laplacian_matrix(G) equals diag(W 1) - W: True
zero eigenvalues: 2   next: 0.775   components: 2   algebraic_connectivity: 0.0
fiedler_vector: graph is not connected.

Two zeros for two pieces, and NetworkX refuses a Fiedler vector. Protein interaction networks nearly always come in pieces: split each one you care about, usually the largest of nx.connected_components, on its own.

The algebraic connectivity bounds every balanced cut. Fiedler named \(\lambda_2\), the smallest Rayleigh quotient \(\mathbf x^{\mathsf T}L\mathbf x / \mathbf x^{\mathsf T}\mathbf x\) over \(\mathbf x \perp \mathbf 1\), the algebraic connectivity in 1973. A split into equal halves is a \(\pm 1\) vector with zero sum, so its quotient \(4\,\mathrm{cut}/n\) cannot fall below \(\lambda_2\):

\[\mathrm{cut} \ge \frac{n\,\lambda_2}{4} .\]

For the club that is 34 × 0.4685 / 4 = 3.98, so no split into 17 and 17 crosses fewer than four friendships; the members' own crosses 11. The bound is loose, but it says what \(\lambda_2\) measures: near zero, the network is nearly in two pieces. The normalized Laplacian below says how near is near.

Rounding, and what it does and does not find. By sign, the Fiedler vector puts 15 members on one side and 19 on the other, cuts ten friendships, and has nodes 2 and 8 on the wrong side. The 15/19 is the zero sum of real entries at work, not equal head counts. A split at the median forces 17 and 17 and recovers the members' own line exactly, but only because the 17 came from the answer. Your own network will not supply it.

Show code

def wrong_side(top):
    return [k for k in G if (k in top) != (side[k] == "Mr. Hi")]

top = {k for k in G if fiedler[k] > np.median(fiedler)}
print(f"bound for a 17/17 split: n lambda_2 / 4 = {n * w[1] / 4:.2f}")
print(f"sign:   {len(hi)}/{n - len(hi)}   wrong {wrong_side(set(hi))}   cut {nx.cut_size(G, hi)}")
print(f"median: {len(top)}/{n - len(top)}   wrong {wrong_side(top)}   cut {nx.cut_size(G, top)}")

bound for a 17/17 split: n lambda_2 / 4 = 3.98
sign:   15/19   wrong [2, 8]   cut 10
median: 17/17   wrong []   cut 11

Without the answer, a split needs a score fair to sides of any size, and the raw cut is not: cutting off node 11, who has one friend, crosses a single edge. The ratio cut divides the cut by the size of each side and adds the two, so a 15/19 split pays about an eighth of its edges and a lone member its one edge in full. The Fiedler vector relaxes this score too: \(\sqrt{|\bar S|/|S|}\) on a side \(S\) and \(-\sqrt{|S|/|\bar S|}\) on the rest sum to zero for any sizes, and their quotient is exactly the ratio cut. A sweep tries every cut point along the sorted Fiedler vector and keeps the lowest score.

Show code

def ratio_cut(S):
    cut = nx.cut_size(G, S, weight=None)
    return cut / len(S) + cut / (n - len(S))

a, b = len(hi), n - len(hi)
z = np.array([np.sqrt(b / a) if fiedler[k] > 0 else -np.sqrt(a / b) for k in G])
print(f"sign split: ratio cut {ratio_cut(hi):.2f}   quotient of its vector {z @ L @ z / (z @ z):.2f}   sum {z.sum():.0f}")
print(f"node 11 alone: cut {nx.cut_size(G, [11])}   ratio cut {ratio_cut([11]):.2f}")
order = np.argsort(-fiedler)
best = min((ratio_cut(order[:i].tolist()), i) for i in range(1, n))
S = sorted(order[:best[1]].tolist())
print(f"sweep: best ratio cut {best[0]:.2f}, cutting off {S} with {nx.cut_size(G, S)} friendships")

sign split: ratio cut 1.19   quotient of its vector 1.19   sum 0
node 11 alone: cut 1   ratio cut 1.03
sweep: best ratio cut 0.94, cutting off [4, 5, 6, 10, 16] with 4 friendships

Counted in edges, the loner beats the sign split ten to one. The ratio cut narrows that to 1.03 against 1.19 without reversing it, and the sweep beats both with 0.94 by cutting off five of Mr. Hi's members, nodes 4, 5, 6, 10, and 16, with four friendships.

What the normalized Laplacian changes. Give each member a mass equal to its degree, so that hubs are heavy. Newton's law becomes \(D\ddot{\mathbf x} = -L\mathbf x\), the modes solve \(L\mathbf x = \lambda D\mathbf x\), and the substitution \(\mathbf v = D^{1/2}\mathbf x\), the \(q_i = \sqrt{m_i}\,x_i\) that made CO₂'s matrix symmetric, turns it into the eigenproblem of the normalized Laplacian

\[L_{\mathrm{sym}} = I - D^{-1/2} A\, D^{-1/2} .\]

The two problems have the same eigenvalues, and their eigenvectors differ entry by entry by the positive factor \(\sqrt{d_i}\), so they have the same signs. Its relaxation approximates the normalized cut, the ratio cut with each side's head count replaced by the sum of its members' degrees: a side counts as large when it holds many friendships, not many members. The eigenvalues lie between 0 and 2 for every network, which gives near zero a scale.

Show code

Lsym = nx.normalized_laplacian_matrix(G, weight=None).toarray()
ws, Vs = eigh(Lsym)
vs = Vs[:, 1] * np.sign(Vs[0, 1])
wg, Vg = eigh(L, np.diag(d))              # the generalized problem L x = lambda D x
print(f"L_sym: lambda_2 = {ws[1]:.4f}   largest = {ws[-1]:.3f}   generalized lambda_2 = {wg[1]:.4f}")
print(f"|correlation| of sqrt(d) x with v: {abs(np.corrcoef(Vg[:, 1] * np.sqrt(d), vs)[0, 1]):.6f}")
print("wrong side by sign:", wrong_side({k for k in G if vs[k] > 0}))

ncut = lambda S: nx.normalized_cut_size(G, S, weight=None)
print(f"normalized cut: sign split {ncut(hi):.3f}   node 11 alone {ncut([11]):.3f}")
order_s = np.argsort(-vs)
best_s = min((ncut(order_s[:i].tolist()), i) for i in range(1, n))
S_s = set(order_s[:best_s[1]].tolist())
print(f"sweep: best normalized cut {best_s[0]:.3f} at {len(S_s)}/{n - len(S_s)}, cut {nx.cut_size(G, S_s)}, wrong {wrong_side(S_s)}")

shuffled = []
for seed in range(20):                    # every member keeps its number of friends, the friends change
    R = nx.Graph(G.edges)
    nx.double_edge_swap(R, nswap=500, max_tries=10**5, seed=seed)
    if nx.is_connected(R):
        shuffled.append(eigh(nx.normalized_laplacian_matrix(R).toarray())[0][1])
print(f"lambda_2 of L_sym in {len(shuffled)} shuffled clubs: {min(shuffled):.2f} to {max(shuffled):.2f}")

L_sym: lambda_2 = 0.1323   largest = 1.715   generalized lambda_2 = 0.1323
|correlation| of sqrt(d) x with v: 1.000000
wrong side by sign: [2, 8]
normalized cut: sign split 0.263   node 11 alone 1.006
sweep: best normalized cut 0.257 at 16/18, cut 10, wrong [8]
lambda_2 of L_sym in 20 shuffled clubs: 0.20 to 0.37

On the club \(\lambda_2 = 0.1323\) from both problems, and the sign split is the same, nodes 2 and 8. The largest eigenvalue is 1.715, against 18.14 for \(L\). Weighed by friendships, node 11 alone scores 1.006 against the sign split's 0.263, so the normalized cut punishes the loner where the ratio cut did not, and its sweep finds 0.257 at 16/18, with only node 8 on the wrong side. Twenty shuffled clubs, each member keeping its number of friends, give \(\lambda_2\) between 0.20 and 0.37, so the club's 0.13 is a real fault line. Shuffle your own network the same way to scale its \(\lambda_2\).

The sign split is a first look. Report the split from a sweep along the \(L_{\mathrm{sym}}\) eigenvector, with its normalized cut, whenever many nodes have only one or two links, as in protein interaction networks, because there the ratio cut rewards shaving off loners.

See it in code

Show code

from matplotlib.colors import LinearSegmentedColormap, TwoSlopeNorm

L_nx = nx.laplacian_matrix(G, weight=None)                # sparse; weight is read by default
w_nx, V_nx = eigh(L_nx.toarray())
f = V_nx[:, 1] * np.sign(V_nx[0, 1])                      # node 0, Mr. Hi, positive
f_nx = nx.fiedler_vector(G, weight=None, seed=0)
f_norm = nx.fiedler_vector(G, weight=None, normalized=True, seed=0)
wrong = wrong_side({k for k in G if f[k] > 0})
print(f"largest difference from the hand-built L: {np.abs(L_nx.toarray() - L).max():.1f}")
print(f"lambda_2: eigh {w_nx[1]:.4f}   algebraic_connectivity {nx.algebraic_connectivity(G, weight=None, seed=0):.4f}")
print(f"correlation of fiedler_vector with eigh's: {np.corrcoef(f_nx, f)[0, 1]:+.3f}")
print(f"normalized=True against the L_sym eigenvector: |correlation| {abs(np.corrcoef(f_norm, vs)[0, 1]):.3f}")
print("wrong side by sign:", wrong)

wW, VW = eigh(nx.laplacian_matrix(G).toarray())           # with Zachary's weights: L = D - W
fW = VW[:, 1] * np.sign(VW[0, 1])
print(f"weighted: lambda_2 = {wW[1]:.3f}   wrong side by sign: {wrong_side({k for k in G if fW[k] > 0})}")

cmap = LinearSegmentedColormap.from_list("fiedler", [ACCENT, "#d6d3cc", SECOND])
norm = TwoSlopeNorm(0, vmin=-np.abs(f).max(), vmax=np.abs(f).max())
fig, (ax, ax2) = plt.subplots(1, 2, figsize=(8.5, 4), width_ratios=[1.3, 1])
nx.draw_networkx_edges(G, pos, ax=ax, edge_color=MUTED, alpha=0.5, width=1)
for club, shape in SHAPES:
    nodes = [k for k in G if side[k] == club]
    nx.draw_networkx_nodes(G, pos, nodelist=nodes, node_shape=shape, ax=ax, node_size=80,
                           node_color=[cmap(norm(f[k])) for k in nodes],
                           edgecolors=[INK if k in wrong else "white" for k in nodes], linewidths=1.5)
for k, offset in zip(wrong, [(-9, 0), (0, -11)]):        # node 2 to the left, node 8 below: clear of their neighbors
    ax.annotate(f"{k}", pos[k], xytext=offset, textcoords="offset points", ha="center", va="center",
                color=INK, bbox=dict(fc="white", ec="none", pad=0.5))
shape_key = [plt.Line2D([], [], ls="", marker=m, ms=8, color=MUTED, label=h.get_label())   # shape only, no color
             for h, m in zip(handles, ["o", "s"])]
ax.legend(handles=shape_key, frameon=False, loc="lower left")
ax.set_axis_off()
fig.colorbar(plt.cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax, shrink=0.7, pad=0.0,
             label="entry of the Fiedler vector")

ranked = np.argsort(-f)
ax2.scatter(np.arange(1, n + 1), f[ranked], c=[cmap(norm(v)) for v in f[ranked]], s=36, zorder=3,
            edgecolors=[INK if k in wrong else "white" for k in ranked], linewidths=1.2)
ax2.axhline(0, color=MUTED, ls="--", lw=1)
ax2.axhline(np.median(f), color=MUTED, lw=1)
n_pos = int((f > 0).sum())
ax2.axvline(n_pos + 0.5, color=MUTED, ls="--", lw=1)
for r, k in enumerate(ranked, start=1):
    if k in wrong:
        ax2.annotate(f"{k}", (r, f[k]), xytext=(6, 4), textcoords="offset points", color=INK)
ax2.text(n_pos + 0.8, np.abs(f).max() * 0.9, "sign change", color=MUTED)
ax2.text(n, np.median(f) + 0.01, "median", color=MUTED, ha="right", va="bottom")
ax2.set(xlabel="rank", ylabel="entry of the Fiedler vector")
fig.tight_layout()
plt.show()

largest difference from the hand-built L: 0.0
lambda_2: eigh 0.4685   algebraic_connectivity 0.4685
correlation of fiedler_vector with eigh's: -1.000
normalized=True against the L_sym eigenvector: |correlation| 1.000
wrong side by sign: [2, 8]
weighted: lambda_2 = 1.187   wrong side by sign: [8]

Left: the karate club, each member colored by its entry of the Fiedler vector, blue positive, red negative, shape by the side joined; nodes 2 and 8 ringed. Right: the 34 entries sorted, 15 above zero and 19 below; nodes 2 and 8, ringed, are the first two below zero but above the median line.

All routes give \(\lambda_2 = 0.4685\), and NetworkX's sparse Laplacian agrees with the hand-built one entry for entry. fiedler_vector returns the same vector with the opposite sign, a correlation of −1.000, which is why the sign is fixed by a node and never trusted as returned. With normalized=True it returns the \(L_{\mathrm{sym}}\) eigenvector, whose signs are those of \(L\)'s. Drop weight=None and NetworkX reads Zachary's counts of shared contexts, \(L = D - W\): \(\lambda_2\) rises to 1.187 and only node 8 stays on the wrong side. That is the one change a reader with confidence-scored interactions makes.

In the figure, the diverging colors run from red through gray to blue, and the two ringed members sit near the gray middle. In the sorted entries they are ranks 16 and 17, just below zero and just above the median: the two places to round disagree on exactly these two. A protein network with thousands of nodes needs one change more: keep the sparse matrix, skip toarray(), and ask scipy.sparse.linalg.eigsh for the few smallest eigenvalues with sigma a little below zero, because the zero mode makes \(L\) singular, as the free edge in Sparse eigenvalues with scipy.sparse.linalg.eigsh does.

Where it shows up

Biology. Modules in protein interaction networks, sets of proteins that act together in a complex or a pathway, are found by spectral clustering with each interaction's confidence score as its weight. Compartments in a food web, groups of species that eat mostly among themselves, show up the same way, with interaction strengths as weights.
Machine learning and data. Spectral clustering (Ng, Jordan, and Weiss 2002) takes the first \(k\) eigenvectors as coordinates for each node or data point and runs k-means on them, which is how two halves become \(k\) modules; scikit-learn's SpectralClustering accepts a weighted adjacency matrix with affinity="precomputed". Laplacian eigenmaps use the same coordinates to draw high-dimensional data in two dimensions.
Image segmentation. Shi and Malik's normalized cuts (2000) make every pixel a node and the similarity of neighboring pixels a weight. The Fiedler vector of the normalized problem separates an object from its background.
Parallel solvers. Recursive spectral bisection (Pothen, Simon, and Liou 1990) splits a finite element mesh by its Fiedler vector, then splits each half again, so that processors share few edges. Multilevel partitioners such as METIS have since taken over, but the goal, a balanced cut with few edges, is the same.
Diffusion and consensus. \(\dot{\mathbf x} = -L\mathbf x\), first order instead of second, is heat flowing along the edges. Every mode decays as \(e^{-\lambda t}\), so \(\lambda_2\) sets how fast a heat pulse, a rumor, or a set of coupled oscillators evens out: in the club the last difference between the halves fades as \(e^{-0.4685\,t}\).
Physics and PDEs. The Laplacian of a grid graph is the five-point stencil of the drum in Sparse eigenvalues with scipy.sparse.linalg.eigsh, with the sign flipped and with free edges where the drum's are clamped. Its Fiedler vector is the lowest mode that is not constant: one half of the plate against the other.

In every case the work is the one done on the club: put the weights into \(L\), count the zero eigenvalues for pieces, and round the Fiedler vector by a sweep with a named score.

The question

The idea: the network as masses on springs

Counting a cut with plus and minus one

Formalization

See it in code

Where it shows up

Further reading