$k$-Way Spectral Clustering as a Graph Comes Apart

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

November 3, 2026

Normalized spectral clustering embeds vertex $i$ as the row $\bfy_i=(\phi_1(i),\ldots,\phi_k(i))$ of the vectors $\bfphi_\ell=\bfD^{-1/2}\bfv_\ell$, where $\bfv_1,\ldots,\bfv_k$ are eigenvectors of the normalized Laplacian $\bfL_{\mathrm{sym}}=\bfI-\bfD^{-1/2}\bfW\bfD^{-1/2}$ for its $k$ smallest eigenvalues, and then runs $k$-means on the rows. The lecture justified this in the ideal case: if the graph has exactly $k$ connected components, the eigenvalue $0$ of $\bfL_{\mathrm{sym}}$ has multiplicity $k$, and the rows collapse to $k$ points, one per component. For a graph that is only nearly disconnected, the claim was that the rows land in $k$ small blobs instead.

Let’s see it! We take a graph with three planted communities and many edges between them, then delete the edges between communities a few at a time until the graph falls apart into three components. Along the way we track the embedding, the spectrum of $\bfL_{\mathrm{sym}}$, and how well $k$-means on the embedding recovers the communities.

import itertools

import networkx as nx
import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from sklearn.cluster import KMeans

Figure style

# A muted ink colour and a translucent grid read correctly against both the
# light and the dark version of this site, so the figures need no theme switch.
INK = "#7c7c78"
GRID = "rgba(128, 128, 128, 0.28)"
COLOR = ["#2a78d6", "#eb6834", "#8e5cc4"]    # one colour per community
CROSS = "#2ca25f"                            # edges between communities


def style(fig, height=420, title=None, showlegend=None):
    """Apply the shared layout to a plot with axes."""
    fig.update_layout(
        template="simple_white",
        paper_bgcolor="rgba(0,0,0,0)",
        plot_bgcolor="rgba(0,0,0,0)",
        font=dict(color=INK, size=13),
        title=dict(text=title, x=0, xanchor="left", font=dict(size=15)),
        margin=dict(l=70, r=30, t=60 if title else 30, b=55),
        height=height,
        showlegend=showlegend,
        hovermode="closest",
    )
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig


def equal_aspect(fig, panels):
    """Make one unit the same length on both axes of the given subplots."""
    for i in panels:
        suffix = "" if i == 1 else str(i)
        fig.layout[f"yaxis{suffix}"].scaleanchor = f"x{suffix}"
        fig.layout[f"xaxis{suffix}"].constrain = "domain"


def segments(edges, pos):
    """Coordinates for drawing the edges as one plotly line trace."""
    xs, ys = [], []
    for i, j in edges:
        xs += [pos[i, 0], pos[j, 0], None]
        ys += [pos[i, 1], pos[j, 1], None]
    return xs, ys

The graph

Let’s consider three communities of $40$, $30$, and $20$ vertices. Each pair of vertices in the same community is joined with probability $0.25$, and each pair in different communities with probability $0.14$. All edges have weight $1$. The communities have unequal sizes on purpose, so that their volumes differ.

sizes = [40, 30, 20]
k = len(sizes)
community = np.repeat(np.arange(k), sizes)    # the planted community of each vertex
n = len(community)
same = community[:, None] == community[None, :]

rng = np.random.default_rng(0)
coin = np.triu(rng.random((n, n)), 1)
coin += coin.T                                # one symmetric coin flip per pair
A = np.where(same, coin < 0.25, coin < 0.14).astype(float)
np.fill_diagonal(A, 0)

# The edges between communities, in a random order in which we will delete them.
cross = [(i, j) for i in range(n) for j in range(i + 1, n) if A[i, j] and not same[i, j]]
cross = [cross[e] for e in rng.permutation(len(cross))]
print(f"{n} vertices, {int(A[same].sum() / 2)} edges inside communities, {len(cross)} between them")

90 vertices, 390 edges inside communities, 362 between them

Draw the graph

# Draw each community with a spring layout of its own edges, around a corner of a triangle.
angles = np.pi / 2 + 2 * np.pi * np.arange(k) / k
corners = 2.4 * np.column_stack([np.cos(angles), np.sin(angles)])
pos = np.zeros((n, 2))
for l in range(k):
    idx = np.flatnonzero(community == l)
    layout = nx.spring_layout(nx.from_numpy_array(A[np.ix_(idx, idx)]), seed=l)
    pos[idx] = corners[l] + np.sqrt(sizes[l] / sizes[0]) * np.array([layout[v] for v in range(len(idx))])
xspan, yspan = np.abs(pos).max(axis=0) + 0.2
vertex_color = np.array(COLOR)[community]

fig = make_subplots(rows=1, cols=2, column_widths=[0.55, 0.45], horizontal_spacing=0.08,
                    subplot_titles=["graph", "weight matrix W"])

# The graph: edges first, so the vertices are drawn on top.
inside = [(i, j) for i in range(n) for j in range(i + 1, n) if A[i, j] and same[i, j]]
for edges, color, width, opacity in [(inside, INK, 0.6, 0.35), (cross, CROSS, 0.8, 0.45)]:
    xs, ys = segments(edges, pos)
    fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines", line=dict(color=color, width=width),
                             opacity=opacity, hoverinfo="skip", showlegend=False), row=1, col=1)
fig.add_trace(go.Scatter(
    x=pos[:, 0], y=pos[:, 1], mode="markers", showlegend=False,
    marker=dict(size=8, color=vertex_color, line=dict(color="white", width=0.8)),
    customdata=np.column_stack([np.arange(n), np.array(sizes)[community], A.sum(axis=1)]),
    hovertemplate="vertex %{customdata[0]}<br>community of %{customdata[1]}"
                  "<br>degree %{customdata[2]}<extra></extra>",
), row=1, col=1)

# The weight matrix: one square dot per nonzero entry, colored like the edge it stands for.
row, col = np.nonzero(A)
fig.add_trace(go.Scatter(
    x=col, y=row, mode="markers", showlegend=False,
    marker=dict(symbol="square", size=3, color=np.where(same[row, col], vertex_color[row], CROSS)),
    hovertemplate="W[%{y}, %{x}] = 1<extra></extra>",
), row=1, col=2)
for b in np.cumsum(sizes)[:-1] - 0.5:
    fig.add_vline(x=b, line=dict(color=INK, width=1, dash="dash"), row=1, col=2)
    fig.add_hline(y=b, line=dict(color=INK, width=1, dash="dash"), row=1, col=2)

# Legend entries.
for l in range(k):
    fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name=f"community of {sizes[l]}",
                             marker=dict(size=10, color=COLOR[l])), row=1, col=1)
fig.add_trace(go.Scatter(x=[None], y=[None], mode="lines", name="edge between communities",
                         line=dict(color=CROSS, width=2)), row=1, col=1)

fig.update_xaxes(visible=False, range=[-xspan, xspan], row=1, col=1)
fig.update_yaxes(visible=False, range=[-yspan, yspan], row=1, col=1)
fig.update_xaxes(range=[-0.5, n - 0.5], title="j", showgrid=False, row=1, col=2)
fig.update_yaxes(range=[n - 0.5, -0.5], title="i", showgrid=False, row=1, col=2)
equal_aspect(fig, [1, 2])
style(fig, height=460, showlegend=True)
fig.update_annotations(font=dict(size=13, color=INK))
fig.update_layout(legend=dict(orientation="h", x=0, y=-0.12, xanchor="left", yanchor="top"))
fig.show()

Figure 1: The graph with every edge. Left: each community drawn around its own corner, with edges inside communities in grey and edges between communities in green. Right: the weight matrix $\bfW$ with vertices ordered by community, a dot at $(i,j)$ for each edge; the dashed lines separate the communities.

Between communities there are almost as many edges as inside them. In Figure 1 the green edges all but hide the three communities, and the off-diagonal blocks of $\bfW$ are only a little sparser than the diagonal ones. A vertex has about $0.25$ of its own community as neighbours and $0.14$ of each other one, so for the smallest community most neighbours of a vertex are outsiders.

The function below keeps only the first $r$ edges of cross and deletes the rest. At $r=0$ the graph is three separate communities.

def graph_with(r):
    """The weight matrix keeping only the first r edges between communities."""
    W = A.copy()
    for i, j in cross[r:]:
        W[i, j] = W[j, i] = 0
    return W


print("smallest degree with no edges between communities:", int(graph_with(0).sum(axis=1).min()))

smallest degree with no edges between communities: 1

No vertex is ever isolated, so $\bfD$ stays invertible as in the lecture.

The embedding

The embedding follows the algorithm, with one simplification. The first eigenvector $\bfphi_1=\mathbf{1}/\sqrt{\operatorname{vol}(V)}$ is constant, so every row $\bfy_i$ has the same first coordinate. Dropping that coordinate changes no distance between rows, and so changes nothing about $k$-means. For $k=3$ we can therefore draw the rows exactly in the plane, as the points $(\phi_2(i),\phi_3(i))$.

When the graph is disconnected, $0$ is a repeated eigenvalue of $\bfL_{\mathrm{sym}}$, and an eigensolver may return a basis of its eigenspace that does not contain $\bfv_1=\bfD^{1/2}\mathbf{1}/\sqrt{\operatorname{vol}(V)}$. The code avoids this with the same trick as in the algebraic connectivity example: it moves the eigenvalue of $\bfv_1$ out of the way, so the eigenvectors it returns first are orthogonal to $\bfv_1$.

def spectral_embedding(W, k):
    """The eigenvalues lam_1 <= lam_2 <= ... of L_sym, and the rows of [phi_2, ..., phi_k]."""
    d = W.sum(axis=1)
    L_sym = np.eye(len(d)) - W / np.sqrt(np.outer(d, d))
    lam = np.linalg.eigvalsh(L_sym)               # increasing
    v1 = np.sqrt(d / d.sum())
    # The spectrum of L_sym lies in [0, 2]. Moving the eigenvalue of v_1 from 0 to 3 puts it
    # last, so the first k - 1 eigenvectors returned are v_2, ..., v_k, all orthogonal to v_1.
    V = np.linalg.eigh(L_sym + 3 * np.outer(v1, v1))[1][:, :k - 1]
    return lam, V / np.sqrt(d)[:, None]           # phi_l = D^{-1/2} v_l


def errors(labels):
    """Which vertices k-means puts in the wrong community, under the best matching of labels."""
    wrong = [np.array(p)[labels] != community for p in itertools.permutations(range(k))]
    return min(wrong, key=np.sum)

Three components

Start at the end, with every edge between communities deleted.

W0 = graph_with(0)
lam0, Y0 = spectral_embedding(W0, k)
print("smallest eigenvalues of L_sym:", np.round(lam0[:5], 4) + 0.0)  # + 0.0 turns -0.0 into 0.0
print("distinct rows of [phi_2, phi_3]:", len(np.unique(Y0.round(10), axis=0)))

smallest eigenvalues of L_sym: [0.     0.     0.     0.4373 0.5128]
distinct rows of [phi_2, phi_3]: 3

As promised, $\lambda_1=\lambda_2=\lambda_3=0$, and the $90$ rows take only $3$ distinct values.

Plot the eigenvectors and the rows

fig = make_subplots(rows=1, cols=2, column_widths=[0.6, 0.4], horizontal_spacing=0.1,
                    subplot_titles=["φ₂ and φ₃ by vertex", "rows (φ₂(i), φ₃(i))"])
for c, name in enumerate(["φ₂", "φ₃"]):
    fig.add_trace(go.Scatter(x=np.arange(n), y=Y0[:, c], mode="markers+lines", name=name,
                             marker=dict(size=5, color=[INK, "#2ca25f"][c]),
                             line=dict(color=[INK, "#2ca25f"][c])), row=1, col=1)
for l in range(k):
    fig.add_trace(go.Scatter(x=Y0[community == l, 0], y=Y0[community == l, 1], mode="markers",
                             name=f"community of {sizes[l]}",
                             marker=dict(size=12, color=COLOR[l], line=dict(color="white", width=1))),
                  row=1, col=2)
fig.update_xaxes(title="vertex", row=1, col=1)
fig.update_xaxes(title="φ₂(i)", row=1, col=2)
fig.update_yaxes(title="φ₃(i)", row=1, col=2)
equal_aspect(fig, [2])
style(fig, height=380, showlegend=True)
fig.update_annotations(font=dict(size=13, color=INK))
fig.show()

Figure 2: The graph with no edges between communities. Left: the eigenvectors $\bfphi_2$ and $\bfphi_3$, by vertex; vertices $0$–$39$, $40$–$69$, and $70$–$89$ are the three communities. Right: the rows $(\phi_2(i),\phi_3(i))$, colored by community; all $90$ of them sit on three points.

The left panel is the point of the remark in lecture that we look at rows, not individual eigenvectors. Each $\bfphi_\ell$ is constant on every community, but neither one is the indicator of a community. It cannot be: $\bfphi_2$ and $\bfphi_3$ are $\bfD$-orthogonal to the constant $\bfphi_1$, that is, $\sum_i\deg(i)\phi_\ell(i)=0$, so each takes both signs. Which combinations of indicators the eigensolver returns is an accident of the basis. The rows, on the right, do not care about the basis: a different orthonormal basis of the eigenspace would only rotate or reflect the triangle.

The three points are not just distinct, they are well spread out. Write $\bfPhi=[\bfphi_1,\ldots,\bfphi_k]$ and $\bfH=[\mathbf{1}_{C_1},\ldots,\mathbf{1}_{C_k}]$, so that $\bfPhi=\bfH\bfC$ as in lecture. The $\bfphi_\ell=\bfD^{-1/2}\bfv_\ell$ satisfy $\bfPhi^\top\bfD\bfPhi=\bfV^\top\bfV=\bfI$, and $\bfH^\top\bfD\bfH=\operatorname{diag}(\operatorname{vol}(C_1),\ldots,\operatorname{vol}(C_k))$, so \[ \bfI=\bfC^\top\bfH^\top\bfD\bfH\bfC=\bfC^\top\operatorname{diag}(\operatorname{vol}(C_\ell))\,\bfC. \] Thus $\operatorname{diag}(\operatorname{vol}(C_\ell))^{1/2}\bfC$ is an orthogonal matrix, and its rows $\sqrt{\operatorname{vol}(C_\ell)}\,\bfc_\ell$ are orthonormal. The point $\bfc_\ell$ that component $C_\ell$ lands on therefore has length $1/\sqrt{\operatorname{vol}(C_\ell)}$, and the $k$ points are mutually orthogonal in $\Rbb^k$, so \[ \|\bfc_\ell-\bfc_m\|_2^2=\frac1{\operatorname{vol}(C_\ell)}+\frac1{\operatorname{vol}(C_m)}. \] Dropping the constant first coordinate preserves these distances, so we can check them in the plane.

points = np.array([Y0[community == l][0] for l in range(k)])
vol = np.array([W0[community == l].sum() for l in range(k)])
for a, b in itertools.combinations(range(k), 2):
    print(f"communities of {sizes[a]} and {sizes[b]}:  "
          f"squared distance {np.sum((points[a] - points[b]) ** 2):.6f},  "
          f"1/vol + 1/vol = {1 / vol[a] + 1 / vol[b]:.6f}")

communities of 40 and 30:  squared distance 0.006358,  1/vol + 1/vol = 0.006358
communities of 40 and 20:  squared distance 0.011792,  1/vol + 1/vol = 0.011792
communities of 30 and 20:  squared distance 0.013434,  1/vol + 1/vol = 0.013434

The smallest community has the smallest volume and so lands farthest from the others.

Deleting the edges between communities

Now run the whole process, from all $362$ edges between communities down to none. The embedding is only defined up to a rotation or reflection of the plane, and near the end the eigensolver is free to pick any basis of a nearly repeated eigenspace. To keep the picture from spinning as the slider moves, each embedding is rotated to best match the next one, by orthogonal Procrustes. A rotation changes no distances, so it does not affect $k$-means.

def align(Y, Y_ref):
    """The rotation or reflection of the rows of Y closest to Y_ref (orthogonal Procrustes)."""
    U, _, Vt = np.linalg.svd(Y.T @ Y_ref)
    return Y @ U @ Vt


remaining = np.linspace(len(cross), 0, 21).round().astype(int)
runs = []
for r in remaining[::-1]:                     # from r = 0 upward, so each step aligns to the last
    W = graph_with(r)
    lam, Y = spectral_embedding(W, k)
    if runs:
        Y = align(Y, runs[-1]["Y"])
    labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(Y)
    runs.append(dict(r=r, lam=lam, Y=Y, wrong=errors(labels)))
runs = runs[::-1]                             # back to the slider order, many edges to none

for run in runs[::4]:
    print(f"{run['r']:3d} edges between communities:  lam_2..lam_5 = {np.round(run['lam'][1:5], 3)},"
          f"  k-means errors {run['wrong'].sum():2d}")

362 edges between communities:  lam_2..lam_5 = [0.584 0.596 0.605 0.626],  k-means errors 50
290 edges between communities:  lam_2..lam_5 = [0.53  0.569 0.585 0.594],  k-means errors 13
217 edges between communities:  lam_2..lam_5 = [0.466 0.51  0.58  0.585],  k-means errors  6
145 edges between communities:  lam_2..lam_5 = [0.347 0.4   0.544 0.574],  k-means errors  2
 72 edges between communities:  lam_2..lam_5 = [0.189 0.251 0.51  0.535],  k-means errors  1
  0 edges between communities:  lam_2..lam_5 = [-0.     0.     0.437  0.513],  k-means errors  0

In the figure below, the slider deletes edges between communities from left to right. The left panel is the graph, laid out as in Figure 1, with the edges between communities in green. The middle panel is the embedding, and the right panel the ten smallest eigenvalues $\lambda_\ell$ of $\bfL_{\mathrm{sym}}$, with the first $k=3$ in blue. Vertices are colored by their planted community, and a cross marks a vertex that $k$-means on the embedding puts in the wrong cluster.

Draw the slider figure

# The vertex positions pos, and xspan, yspan, vertex_color, are those of the graph figure above.
n_lam = 10
emb_span = 1.1 * max(np.abs(run["Y"]).max() for run in runs)

fig = make_subplots(rows=1, cols=3, column_widths=[0.37, 0.37, 0.26], horizontal_spacing=0.06,
                    subplot_titles=["graph", "embedding (φ₂(i), φ₃(i))", "eigenvalues λₗ of L_sym"])

# Traces drawn at every step: the edges inside communities, and legend entries.
xs, ys = segments([(i, j) for i in range(n) for j in range(i + 1, n) if A[i, j] and same[i, j]], pos)
fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines", line=dict(color=INK, width=0.6), opacity=0.35,
                         hoverinfo="skip", showlegend=False), row=1, col=1)
for l in range(k):
    fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name=f"community of {sizes[l]}",
                             marker=dict(size=10, color=COLOR[l])), row=1, col=1)
fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name="misassigned by k-means",
                         marker=dict(size=10, color=INK, symbol="x")), row=1, col=1)
n_static = len(fig.data)

# Four traces per step: edges between communities, graph vertices, embedding, eigenvalues.
for run in runs:
    symbol = np.where(run["wrong"], "x", "circle")
    xs, ys = segments(cross[:run["r"]], pos)
    fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines", line=dict(color=CROSS, width=0.8), opacity=0.45,
                             hoverinfo="skip", showlegend=False), row=1, col=1)
    fig.add_trace(go.Scatter(
        x=pos[:, 0], y=pos[:, 1], mode="markers", showlegend=False,
        marker=dict(size=8, color=vertex_color, symbol=symbol, line=dict(color="white", width=0.8)),
        customdata=np.column_stack([np.arange(n), np.array(sizes)[community]]),
        hovertemplate="vertex %{customdata[0]}<br>community of %{customdata[1]}<extra></extra>",
    ), row=1, col=1)
    fig.add_trace(go.Scatter(
        x=run["Y"][:, 0], y=run["Y"][:, 1], mode="markers", showlegend=False,
        marker=dict(size=8, color=vertex_color, symbol=symbol, line=dict(color="white", width=0.8)),
        customdata=np.arange(n),
        hovertemplate="vertex %{customdata}<br>φ₂ = %{x:.4f}<br>φ₃ = %{y:.4f}<extra></extra>",
    ), row=1, col=2)
    fig.add_trace(go.Bar(
        x=np.arange(1, n_lam + 1), y=run["lam"][:n_lam], showlegend=False,
        marker=dict(color=[COLOR[0] if l < k else INK for l in range(n_lam)]),
        hovertemplate="λ%{x} = %{y:.3f}<extra></extra>",
    ), row=1, col=3)


def visible(step):
    return [True] * n_static + [s == step for s in range(len(runs)) for _ in range(4)]


def title(run):
    return (f"{run['r']} edges between communities:  "
            f"λ₄ − λ₃ = {run['lam'][3] - run['lam'][2]:.2f},  k-means errors: {run['wrong'].sum()}")


for i, flag in enumerate(visible(0)):
    fig.data[i].visible = flag
steps = [dict(method="update", label=str(run["r"]),
              args=[{"visible": visible(s)}, {"title.text": title(run)}])
         for s, run in enumerate(runs)]
fig.update_xaxes(visible=False, range=[-xspan, xspan], row=1, col=1)
fig.update_yaxes(visible=False, range=[-yspan, yspan], row=1, col=1)
fig.update_xaxes(range=[-emb_span, emb_span], title="φ₂(i)", row=1, col=2)
fig.update_yaxes(range=[-emb_span, emb_span], title="φ₃(i)", row=1, col=2)
fig.update_xaxes(range=[0.4, n_lam + 0.6], dtick=1, title="ℓ", row=1, col=3)
fig.update_yaxes(range=[0, 1.05 * max(run["lam"][n_lam - 1] for run in runs)], row=1, col=3)
equal_aspect(fig, [1, 2])
style(fig, height=560, title=title(runs[0]), showlegend=True)
fig.update_annotations(font=dict(size=13, color=INK))
fig.update_layout(
    bargap=0.3,
    margin=dict(b=40),
    legend=dict(orientation="h", x=0, y=-0.2, xanchor="left", yanchor="top"),
    sliders=[dict(active=0, steps=steps, pad=dict(t=95), font=dict(color=INK),
                  currentvalue=dict(prefix="edges between communities: ", font=dict(color=INK)))],
)
fig.show()

Figure 3: Normalized spectral clustering with $k=3$ as the edges between three planted communities are deleted; use the slider to change how many are left. Left: the graph, with edges between communities in green. Middle: the rows $(\phi_2(i),\phi_3(i))$ of the embedding. Right: the ten smallest eigenvalues of $\bfL_{\mathrm{sym}}$. Colors are the planted communities, and crosses are vertices that $k$-means on the embedding misassigns.

Move the slider from left to right.

With all $362$ edges there is no cluster structure for the spectrum to find. The eigenvalues $\lambda_2,\lambda_3,\lambda_4,\ldots$ are all about the same, the embedding is a single cloud, and $k$-means gets more than half the vertices wrong.
As edges are deleted, $\lambda_2$ and $\lambda_3$ fall toward $0$ while $\lambda_4$ falls much more slowly, and the cloud pulls apart into three blobs, one per community. By the time half the edges between communities are gone, $k$-means makes only a few mistakes, and those are vertices with many neighbours in another community.
At $0$ edges, $\lambda_2=\lambda_3=0$ exactly, and each blob has shrunk to one of the three points of Figure 2.

The blobs shrink gradually, not all at once. Each remaining edge between communities pulls its two endpoints toward each other in the embedding, so a vertex with several such edges sits between two blobs until the last of them is deleted.

--- title: "$k$-Way Spectral Clustering as a Graph Comes Apart" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-11-03 number-sections: false jupyter: python3 format: html: toc: true toc-depth: 2 code-line-numbers: true code-tools: true out-width: "100%" column: page execute: echo: true draft: false --- [Normalized spectral clustering](../../lectures/lecture-16.qmd#alg-normalized-spectral-clustering) embeds vertex $i$ as the row $\bfy_i=(\phi_1(i),\ldots,\phi_k(i))$ of the vectors $\bfphi_\ell=\bfD^{-1/2}\bfv_\ell$, where $\bfv_1,\ldots,\bfv_k$ are eigenvectors of the normalized Laplacian $\bfL_{\mathrm{sym}}=\bfI-\bfD^{-1/2}\bfW\bfD^{-1/2}$ for its $k$ smallest eigenvalues, and then runs $k$-means on the rows. The lecture justified this in the ideal case: if the graph has exactly $k$ connected components, [the eigenvalue $0$ of $\bfL_{\mathrm{sym}}$ has multiplicity $k$](../../lectures/lecture-16.qmd#prp-normalized-ideal-embedding), and the rows collapse to $k$ points, one per component. For a graph that is only *nearly* disconnected, the claim was that the rows land in $k$ small blobs instead. Let's see it! We take a graph with three planted communities and many edges between them, then delete the edges between communities a few at a time until the graph falls apart into three components. Along the way we track the embedding, the spectrum of $\bfL_{\mathrm{sym}}$, and how well $k$-means on the embedding recovers the communities. ```{python} import itertools import networkx as nx import numpy as np import plotly.graph_objects as go from plotly.subplots import make_subplots from sklearn.cluster import KMeans ``` ```{python} #| code-fold: true #| code-summary: "Figure style" # A muted ink colour and a translucent grid read correctly against both the # light and the dark version of this site, so the figures need no theme switch. INK = "#7c7c78" GRID = "rgba(128, 128, 128, 0.28)" COLOR = ["#2a78d6", "#eb6834", "#8e5cc4"] # one colour per community CROSS = "#2ca25f" # edges between communities def style(fig, height=420, title=None, showlegend=None): """Apply the shared layout to a plot with axes.""" fig.update_layout( template="simple_white", paper_bgcolor="rgba(0,0,0,0)", plot_bgcolor="rgba(0,0,0,0)", font=dict(color=INK, size=13), title=dict(text=title, x=0, xanchor="left", font=dict(size=15)), margin=dict(l=70, r=30, t=60 if title else 30, b=55), height=height, showlegend=showlegend, hovermode="closest", ) fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) return fig def equal_aspect(fig, panels): """Make one unit the same length on both axes of the given subplots.""" for i in panels: suffix = "" if i == 1 else str(i) fig.layout[f"yaxis{suffix}"].scaleanchor = f"x{suffix}" fig.layout[f"xaxis{suffix}"].constrain = "domain" def segments(edges, pos): """Coordinates for drawing the edges as one plotly line trace.""" xs, ys = [], [] for i, j in edges: xs += [pos[i, 0], pos[j, 0], None] ys += [pos[i, 1], pos[j, 1], None] return xs, ys ``` ## The graph {#sec-graph} Let's consider three communities of $40$, $30$, and $20$ vertices. Each pair of vertices in the same community is joined with probability $0.25$, and each pair in different communities with probability $0.14$. All edges have weight $1$. The communities have unequal sizes on purpose, so that their volumes differ. ```{python} sizes = [40, 30, 20] k = len(sizes) community = np.repeat(np.arange(k), sizes) # the planted community of each vertex n = len(community) same = community[:, None] == community[None, :] rng = np.random.default_rng(0) coin = np.triu(rng.random((n, n)), 1) coin += coin.T # one symmetric coin flip per pair A = np.where(same, coin < 0.25, coin < 0.14).astype(float) np.fill_diagonal(A, 0) # The edges between communities, in a random order in which we will delete them. cross = [(i, j) for i in range(n) for j in range(i + 1, n) if A[i, j] and not same[i, j]] cross = [cross[e] for e in rng.permutation(len(cross))] print(f"{n} vertices, {int(A[same].sum() / 2)} edges inside communities, {len(cross)} between them") ``` ```{python} #| label: fig-graph #| fig-cap: "The graph with every edge. Left: each community drawn around its own corner, with edges inside communities in grey and edges between communities in green. Right: the weight matrix $\\bfW$ with vertices ordered by community, a dot at $(i,j)$ for each edge; the dashed lines separate the communities." #| fig-alt: "Left panel: three clusters of vertices, 40 blue at the top, 30 orange at the lower left, and 20 purple at the lower right, each with grey edges inside it, and a dense web of green edges running between every pair of clusters. Right panel: a 90 by 90 grid of dots for the weight matrix, split by dashed lines into a 3 by 3 grid of blocks. The three diagonal blocks are dotted in blue, orange, and purple, and the off-diagonal blocks in green, only a little more sparsely than the diagonal ones." #| code-fold: true #| code-summary: "Draw the graph" # Draw each community with a spring layout of its own edges, around a corner of a triangle. angles = np.pi / 2 + 2 * np.pi * np.arange(k) / k corners = 2.4 * np.column_stack([np.cos(angles), np.sin(angles)]) pos = np.zeros((n, 2)) for l in range(k): idx = np.flatnonzero(community == l) layout = nx.spring_layout(nx.from_numpy_array(A[np.ix_(idx, idx)]), seed=l) pos[idx] = corners[l] + np.sqrt(sizes[l] / sizes[0]) * np.array([layout[v] for v in range(len(idx))]) xspan, yspan = np.abs(pos).max(axis=0) + 0.2 vertex_color = np.array(COLOR)[community] fig = make_subplots(rows=1, cols=2, column_widths=[0.55, 0.45], horizontal_spacing=0.08, subplot_titles=["graph", "weight matrix W"]) # The graph: edges first, so the vertices are drawn on top. inside = [(i, j) for i in range(n) for j in range(i + 1, n) if A[i, j] and same[i, j]] for edges, color, width, opacity in [(inside, INK, 0.6, 0.35), (cross, CROSS, 0.8, 0.45)]: xs, ys = segments(edges, pos) fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines", line=dict(color=color, width=width), opacity=opacity, hoverinfo="skip", showlegend=False), row=1, col=1) fig.add_trace(go.Scatter( x=pos[:, 0], y=pos[:, 1], mode="markers", showlegend=False, marker=dict(size=8, color=vertex_color, line=dict(color="white", width=0.8)), customdata=np.column_stack([np.arange(n), np.array(sizes)[community], A.sum(axis=1)]), hovertemplate="vertex %{customdata[0]} community of %{customdata[1]}" " degree %{customdata[2]}<extra></extra>", ), row=1, col=1) # The weight matrix: one square dot per nonzero entry, colored like the edge it stands for. row, col = np.nonzero(A) fig.add_trace(go.Scatter( x=col, y=row, mode="markers", showlegend=False, marker=dict(symbol="square", size=3, color=np.where(same[row, col], vertex_color[row], CROSS)), hovertemplate="W[%{y}, %{x}] = 1<extra></extra>", ), row=1, col=2) for b in np.cumsum(sizes)[:-1] - 0.5: fig.add_vline(x=b, line=dict(color=INK, width=1, dash="dash"), row=1, col=2) fig.add_hline(y=b, line=dict(color=INK, width=1, dash="dash"), row=1, col=2) # Legend entries. for l in range(k): fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name=f"community of {sizes[l]}", marker=dict(size=10, color=COLOR[l])), row=1, col=1) fig.add_trace(go.Scatter(x=[None], y=[None], mode="lines", name="edge between communities", line=dict(color=CROSS, width=2)), row=1, col=1) fig.update_xaxes(visible=False, range=[-xspan, xspan], row=1, col=1) fig.update_yaxes(visible=False, range=[-yspan, yspan], row=1, col=1) fig.update_xaxes(range=[-0.5, n - 0.5], title="j", showgrid=False, row=1, col=2) fig.update_yaxes(range=[n - 0.5, -0.5], title="i", showgrid=False, row=1, col=2) equal_aspect(fig, [1, 2]) style(fig, height=460, showlegend=True) fig.update_annotations(font=dict(size=13, color=INK)) fig.update_layout(legend=dict(orientation="h", x=0, y=-0.12, xanchor="left", yanchor="top")) fig.show() ``` Between communities there are almost as many edges as inside them. In @fig-graph the green edges all but hide the three communities, and the off-diagonal blocks of $\bfW$ are only a little sparser than the diagonal ones. A vertex has about $0.25$ of its own community as neighbours and $0.14$ of each other one, so for the smallest community most neighbours of a vertex are outsiders. The function below keeps only the first $r$ edges of `cross` and deletes the rest. At $r=0$ the graph is three separate communities. ```{python} def graph_with(r): """The weight matrix keeping only the first r edges between communities.""" W = A.copy() for i, j in cross[r:]: W[i, j] = W[j, i] = 0 return W print("smallest degree with no edges between communities:", int(graph_with(0).sum(axis=1).min())) ``` No vertex is ever isolated, so $\bfD$ stays invertible as in the lecture. ## The embedding {#sec-embedding} The embedding follows [the algorithm](../../lectures/lecture-16.qmd#alg-normalized-spectral-clustering), with one simplification. The first eigenvector $\bfphi_1=\mathbf{1}/\sqrt{\operatorname{vol}(V)}$ is constant, so every row $\bfy_i$ has the same first coordinate. Dropping that coordinate changes no distance between rows, and so changes nothing about $k$-means. For $k=3$ we can therefore draw the rows exactly in the plane, as the points $(\phi_2(i),\phi_3(i))$. When the graph is disconnected, $0$ is a repeated eigenvalue of $\bfL_{\mathrm{sym}}$, and an eigensolver may return a basis of its eigenspace that does not contain $\bfv_1=\bfD^{1/2}\mathbf{1}/\sqrt{\operatorname{vol}(V)}$. The code avoids this with the same trick as in the [algebraic connectivity example](../algebraic-connectivity/algebraic-connectivity.qmd#sec-bridges): it moves the eigenvalue of $\bfv_1$ out of the way, so the eigenvectors it returns first are orthogonal to $\bfv_1$. ```{python} def spectral_embedding(W, k): """The eigenvalues lam_1 <= lam_2 <= ... of L_sym, and the rows of [phi_2, ..., phi_k].""" d = W.sum(axis=1) L_sym = np.eye(len(d)) - W / np.sqrt(np.outer(d, d)) lam = np.linalg.eigvalsh(L_sym) # increasing v1 = np.sqrt(d / d.sum()) # The spectrum of L_sym lies in [0, 2]. Moving the eigenvalue of v_1 from 0 to 3 puts it # last, so the first k - 1 eigenvectors returned are v_2, ..., v_k, all orthogonal to v_1. V = np.linalg.eigh(L_sym + 3 * np.outer(v1, v1))[1][:, :k - 1] return lam, V / np.sqrt(d)[:, None] # phi_l = D^{-1/2} v_l def errors(labels): """Which vertices k-means puts in the wrong community, under the best matching of labels.""" wrong = [np.array(p)[labels] != community for p in itertools.permutations(range(k))] return min(wrong, key=np.sum) ``` ### Three components {#sec-three-components} Start at the end, with every edge between communities deleted. ```{python} W0 = graph_with(0) lam0, Y0 = spectral_embedding(W0, k) print("smallest eigenvalues of L_sym:", np.round(lam0[:5], 4) + 0.0) # + 0.0 turns -0.0 into 0.0 print("distinct rows of [phi_2, phi_3]:", len(np.unique(Y0.round(10), axis=0))) ``` As promised, $\lambda_1=\lambda_2=\lambda_3=0$, and the $90$ rows take only $3$ distinct values. ```{python} #| label: fig-three-components #| fig-cap: "The graph with no edges between communities. Left: the eigenvectors $\\bfphi_2$ and $\\bfphi_3$, by vertex; vertices $0$–$39$, $40$–$69$, and $70$–$89$ are the three communities. Right: the rows $(\\phi_2(i),\\phi_3(i))$, colored by community; all $90$ of them sit on three points." #| fig-alt: "Left panel: two step-shaped curves over vertex index 0 to 89, each constant on the three blocks of vertices 0 to 39, 40 to 69, and 70 to 89. Each curve takes both positive and negative values, so neither is the indicator of a single community. Right panel: a scatter plot with only three visible points, one blue, one orange, one purple, forming a triangle around the origin." #| code-fold: true #| code-summary: "Plot the eigenvectors and the rows" fig = make_subplots(rows=1, cols=2, column_widths=[0.6, 0.4], horizontal_spacing=0.1, subplot_titles=["φ₂ and φ₃ by vertex", "rows (φ₂(i), φ₃(i))"]) for c, name in enumerate(["φ₂", "φ₃"]): fig.add_trace(go.Scatter(x=np.arange(n), y=Y0[:, c], mode="markers+lines", name=name, marker=dict(size=5, color=[INK, "#2ca25f"][c]), line=dict(color=[INK, "#2ca25f"][c])), row=1, col=1) for l in range(k): fig.add_trace(go.Scatter(x=Y0[community == l, 0], y=Y0[community == l, 1], mode="markers", name=f"community of {sizes[l]}", marker=dict(size=12, color=COLOR[l], line=dict(color="white", width=1))), row=1, col=2) fig.update_xaxes(title="vertex", row=1, col=1) fig.update_xaxes(title="φ₂(i)", row=1, col=2) fig.update_yaxes(title="φ₃(i)", row=1, col=2) equal_aspect(fig, [2]) style(fig, height=380, showlegend=True) fig.update_annotations(font=dict(size=13, color=INK)) fig.show() ``` The left panel is the point of the [remark in lecture](../../lectures/lecture-16.qmd#more-than-two-clusters) that we look at rows, not individual eigenvectors. Each $\bfphi_\ell$ is constant on every community, but neither one is the indicator of a community. It cannot be: $\bfphi_2$ and $\bfphi_3$ are $\bfD$-orthogonal to the constant $\bfphi_1$, that is, $\sum_i\deg(i)\phi_\ell(i)=0$, so each takes both signs. Which combinations of indicators the eigensolver returns is an accident of the basis. The rows, on the right, do not care about the basis: a different orthonormal basis of the eigenspace would only rotate or reflect the triangle. The three points are not just distinct, they are well spread out. Write $\bfPhi=[\bfphi_1,\ldots,\bfphi_k]$ and $\bfH=[\mathbf{1}_{C_1},\ldots,\mathbf{1}_{C_k}]$, so that $\bfPhi=\bfH\bfC$ as in lecture. The $\bfphi_\ell=\bfD^{-1/2}\bfv_\ell$ satisfy $\bfPhi^\top\bfD\bfPhi=\bfV^\top\bfV=\bfI$, and $\bfH^\top\bfD\bfH=\operatorname{diag}(\operatorname{vol}(C_1),\ldots,\operatorname{vol}(C_k))$, so $$ \bfI=\bfC^\top\bfH^\top\bfD\bfH\bfC=\bfC^\top\operatorname{diag}(\operatorname{vol}(C_\ell))\,\bfC. $$ Thus $\operatorname{diag}(\operatorname{vol}(C_\ell))^{1/2}\bfC$ is an orthogonal matrix, and its rows $\sqrt{\operatorname{vol}(C_\ell)}\,\bfc_\ell$ are orthonormal. The point $\bfc_\ell$ that component $C_\ell$ lands on therefore has length $1/\sqrt{\operatorname{vol}(C_\ell)}$, and the $k$ points are mutually orthogonal in $\Rbb^k$, so $$ \|\bfc_\ell-\bfc_m\|_2^2=\frac1{\operatorname{vol}(C_\ell)}+\frac1{\operatorname{vol}(C_m)}. $$ Dropping the constant first coordinate preserves these distances, so we can check them in the plane. ```{python} points = np.array([Y0[community == l][0] for l in range(k)]) vol = np.array([W0[community == l].sum() for l in range(k)]) for a, b in itertools.combinations(range(k), 2): print(f"communities of {sizes[a]} and {sizes[b]}: " f"squared distance {np.sum((points[a] - points[b]) ** 2):.6f}, " f"1/vol + 1/vol = {1 / vol[a] + 1 / vol[b]:.6f}") ``` The smallest community has the smallest volume and so lands farthest from the others. ## Deleting the edges between communities {#sec-slider} Now run the whole process, from all $362$ edges between communities down to none. The embedding is only defined up to a rotation or reflection of the plane, and near the end the eigensolver is free to pick any basis of a nearly repeated eigenspace. To keep the picture from spinning as the slider moves, each embedding is rotated to best match the next one, by [orthogonal Procrustes](https://en.wikipedia.org/wiki/Orthogonal_Procrustes_problem). A rotation changes no distances, so it does not affect $k$-means. ```{python} def align(Y, Y_ref): """The rotation or reflection of the rows of Y closest to Y_ref (orthogonal Procrustes).""" U, _, Vt = np.linalg.svd(Y.T @ Y_ref) return Y @ U @ Vt remaining = np.linspace(len(cross), 0, 21).round().astype(int) runs = [] for r in remaining[::-1]: # from r = 0 upward, so each step aligns to the last W = graph_with(r) lam, Y = spectral_embedding(W, k) if runs: Y = align(Y, runs[-1]["Y"]) labels = KMeans(n_clusters=k, n_init=10, random_state=0).fit_predict(Y) runs.append(dict(r=r, lam=lam, Y=Y, wrong=errors(labels))) runs = runs[::-1] # back to the slider order, many edges to none for run in runs[::4]: print(f"{run['r']:3d} edges between communities: lam_2..lam_5 = {np.round(run['lam'][1:5], 3)}," f" k-means errors {run['wrong'].sum():2d}") ``` In the figure below, the slider deletes edges between communities from left to right. The left panel is the graph, laid out as in @fig-graph, with the edges between communities in green. The middle panel is the embedding, and the right panel the ten smallest eigenvalues $\lambda_\ell$ of $\bfL_{\mathrm{sym}}$, with the first $k=3$ in blue. Vertices are colored by their planted community, and a cross marks a vertex that $k$-means on the embedding puts in the wrong cluster. ```{python} #| label: fig-slider #| fig-cap: "Normalized spectral clustering with $k=3$ as the edges between three planted communities are deleted; use the slider to change how many are left. Left: the graph, with edges between communities in green. Middle: the rows $(\\phi_2(i),\\phi_3(i))$ of the embedding. Right: the ten smallest eigenvalues of $\\bfL_{\\mathrm{sym}}$. Colors are the planted communities, and crosses are vertices that $k$-means on the embedding misassigns." #| fig-alt: "Three panels with a slider below from 362 edges between communities down to 0. Left: three clusters of vertices in blue, orange, and purple, joined by green edges that thin out as the slider moves right. Middle: a scatter of the embedded vertices. With all 362 edges it is one mixed cloud with many crosses; as edges are deleted it separates into three blobs, one per color, that shrink to three single points at 0 edges. Right: ten bars for the eigenvalues of the normalized Laplacian. At first all bars after the first are near 0.6; as edges are deleted the second and third bars fall to 0 while the fourth falls only to about 0.44, opening a gap." #| code-fold: true #| code-summary: "Draw the slider figure" # The vertex positions pos, and xspan, yspan, vertex_color, are those of the graph figure above. n_lam = 10 emb_span = 1.1 * max(np.abs(run["Y"]).max() for run in runs) fig = make_subplots(rows=1, cols=3, column_widths=[0.37, 0.37, 0.26], horizontal_spacing=0.06, subplot_titles=["graph", "embedding (φ₂(i), φ₃(i))", "eigenvalues λₗ of L_sym"]) # Traces drawn at every step: the edges inside communities, and legend entries. xs, ys = segments([(i, j) for i in range(n) for j in range(i + 1, n) if A[i, j] and same[i, j]], pos) fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines", line=dict(color=INK, width=0.6), opacity=0.35, hoverinfo="skip", showlegend=False), row=1, col=1) for l in range(k): fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name=f"community of {sizes[l]}", marker=dict(size=10, color=COLOR[l])), row=1, col=1) fig.add_trace(go.Scatter(x=[None], y=[None], mode="markers", name="misassigned by k-means", marker=dict(size=10, color=INK, symbol="x")), row=1, col=1) n_static = len(fig.data) # Four traces per step: edges between communities, graph vertices, embedding, eigenvalues. for run in runs: symbol = np.where(run["wrong"], "x", "circle") xs, ys = segments(cross[:run["r"]], pos) fig.add_trace(go.Scatter(x=xs, y=ys, mode="lines", line=dict(color=CROSS, width=0.8), opacity=0.45, hoverinfo="skip", showlegend=False), row=1, col=1) fig.add_trace(go.Scatter( x=pos[:, 0], y=pos[:, 1], mode="markers", showlegend=False, marker=dict(size=8, color=vertex_color, symbol=symbol, line=dict(color="white", width=0.8)), customdata=np.column_stack([np.arange(n), np.array(sizes)[community]]), hovertemplate="vertex %{customdata[0]} community of %{customdata[1]}<extra></extra>", ), row=1, col=1) fig.add_trace(go.Scatter( x=run["Y"][:, 0], y=run["Y"][:, 1], mode="markers", showlegend=False, marker=dict(size=8, color=vertex_color, symbol=symbol, line=dict(color="white", width=0.8)), customdata=np.arange(n), hovertemplate="vertex %{customdata} φ₂ = %{x:.4f} φ₃ = %{y:.4f}<extra></extra>", ), row=1, col=2) fig.add_trace(go.Bar( x=np.arange(1, n_lam + 1), y=run["lam"][:n_lam], showlegend=False, marker=dict(color=[COLOR[0] if l < k else INK for l in range(n_lam)]), hovertemplate="λ%{x} = %{y:.3f}<extra></extra>", ), row=1, col=3) def visible(step): return [True] * n_static + [s == step for s in range(len(runs)) for _ in range(4)] def title(run): return (f"{run['r']} edges between communities: " f"λ₄ − λ₃ = {run['lam'][3] - run['lam'][2]:.2f}, k-means errors: {run['wrong'].sum()}") for i, flag in enumerate(visible(0)): fig.data[i].visible = flag steps = [dict(method="update", label=str(run["r"]), args=[{"visible": visible(s)}, {"title.text": title(run)}]) for s, run in enumerate(runs)] fig.update_xaxes(visible=False, range=[-xspan, xspan], row=1, col=1) fig.update_yaxes(visible=False, range=[-yspan, yspan], row=1, col=1) fig.update_xaxes(range=[-emb_span, emb_span], title="φ₂(i)", row=1, col=2) fig.update_yaxes(range=[-emb_span, emb_span], title="φ₃(i)", row=1, col=2) fig.update_xaxes(range=[0.4, n_lam + 0.6], dtick=1, title="ℓ", row=1, col=3) fig.update_yaxes(range=[0, 1.05 * max(run["lam"][n_lam - 1] for run in runs)], row=1, col=3) equal_aspect(fig, [1, 2]) style(fig, height=560, title=title(runs[0]), showlegend=True) fig.update_annotations(font=dict(size=13, color=INK)) fig.update_layout( bargap=0.3, margin=dict(b=40), legend=dict(orientation="h", x=0, y=-0.2, xanchor="left", yanchor="top"), sliders=[dict(active=0, steps=steps, pad=dict(t=95), font=dict(color=INK), currentvalue=dict(prefix="edges between communities: ", font=dict(color=INK)))], ) fig.show() ``` Move the slider from left to right. - *With all $362$ edges* there is no cluster structure for the spectrum to find. The eigenvalues $\lambda_2,\lambda_3,\lambda_4,\ldots$ are all about the same, the embedding is a single cloud, and $k$-means gets more than half the vertices wrong. - *As edges are deleted*, $\lambda_2$ and $\lambda_3$ fall toward $0$ while $\lambda_4$ falls much more slowly, and the cloud pulls apart into three blobs, one per community. By the time half the edges between communities are gone, $k$-means makes only a few mistakes, and those are vertices with many neighbours in another community. - *At $0$ edges*, $\lambda_2=\lambda_3=0$ exactly, and each blob has shrunk to one of the three points of @fig-three-components. The blobs shrink gradually, not all at once. Each remaining edge between communities pulls its two endpoints toward each other in the embedding, so a vertex with several such edges sits between two blobs until the last of them is deleted.