Unshredding an Image with One Eigenvector

Math 123: Mathematical Aspects of Data Analysis — Fall 2026

Published

November 4, 2026

Cut a photograph into single columns of pixels, shuffle the columns, and hand them back: can we put them in order again? Neighbouring columns of a photograph look almost the same, so the columns are points lying along a curve in a high-dimensional space, and the order we want is position along that curve. A one-dimensional Laplacian eigenmap gives each column a single number, and sorting the columns by that number puts them back in order. Using an eigenvector to put things in order, instead of to cluster or to plot them, is called seriation.

import base64
import io

import numpy as np
import plotly.graph_objects as go
from PIL import Image
from scipy.ndimage import binary_dilation
from scipy.sparse.csgraph import connected_components
from skimage import data
from sklearn.neighbors import kneighbors_graph

Figure style and image helpers

INK = "#7c7c78"                               # muted colors that read in light and dark mode
GRID = "rgba(128, 128, 128, 0.28)"
BLUE, ORANGE = "#2a78d6", "#eb6834"


def style(fig, height):
    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),
                      height=height, margin=dict(l=60, r=20, t=30, b=50), hovermode="closest")
    fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID)
    return fig


def to_rgba(image):
    """An H x W (grey) or H x W x 3 (colour) image in [0, 1] as opaque RGBA bytes."""
    rgb = np.repeat(image[..., None], 3, axis=2) if image.ndim == 2 else image
    alpha = np.ones(rgb.shape[:2] + (1,))
    return (np.clip(np.concatenate([rgb, alpha], axis=2), 0, 1) * 255).astype(np.uint8)


def image_uri(rgba):
    """Encode RGBA bytes as a data URI: JPEG for opaque photographs, PNG otherwise.

    An image is far smaller than the same picture as a heatmap of numbers.
    """
    buffer = io.BytesIO()
    if rgba[..., 3].min() == 255:
        Image.fromarray(rgba[..., :3], "RGB").save(buffer, "JPEG", quality=90)
        kind = "jpeg"
    else:
        Image.fromarray(rgba, "RGBA").save(buffer, "PNG")
        kind = "png"
    return f"data:image/{kind};base64," + base64.b64encode(buffer.getvalue()).decode()


def image_figure(images, labels, height, gap=16):
    """Show RGBA images of equal height side by side, with a label under each."""
    fig, left, centres = go.Figure(), 0, []
    for image in images:
        fig.add_trace(go.Image(source=image_uri(image), x0=left, hoverinfo="skip"))
        centres.append(left + (image.shape[1] - 1) / 2)
        left += image.shape[1] + gap
    fig.update_xaxes(tickvals=centres, ticktext=labels, showgrid=False, zeroline=False,
                     showline=False, ticks="")
    fig.update_yaxes(showticklabels=False, showgrid=False, zeroline=False, showline=False)
    return style(fig, height)


def adjacency_image(W):
    """Draw each edge of W as a small ink-coloured square on a faint grey background."""
    rgba = np.empty(W.shape + (4,), dtype=np.uint8)
    rgba[:] = (128, 128, 128, 30)
    rgba[binary_dilation(W > 0, np.ones((3, 3)))] = (124, 124, 120, 255)
    return rgba

The shredded photograph

We use the cat from skimage.data, a $300\times451$ colour photograph. Each column is $300$ pixels with three colour channels, so it is a vector in $\Rbb^{900}$, and the image is $451$ points in $\Rbb^{900}$. We shuffle the columns with a random permutation. From here on the method only sees shredded; we keep perm only to check the answer at the end.

photo = data.chelsea() / 255                  # (300, 451, 3), brightness in [0, 1]
n_rows, n_cols, _ = photo.shape

rng = np.random.default_rng(0)
perm = rng.permutation(n_cols)                # column j of shredded is column perm[j] of photo
shredded = photo[:, perm]

F = shredded.transpose(1, 0, 2).reshape(n_cols, -1)   # one row per column of pixels
print(F.shape)

(451, 900)

Plot the photograph and the shredded version

image_figure([to_rgba(photo), to_rgba(shredded)], ["photograph", "shredded"], height=380).show()

Figure 1: Left: the photograph. Right: the same $451$ columns of pixels in a random order.

The graph and its eigenmap

We join each column to its $4$ nearest neighbours in $\Rbb^{900}$, in Euclidean distance, and keep an edge if either end chose the other. The weights are all $1$. The one-dimensional eigenmap is the single coordinate $\bfphi_2=\bfD^{-1/2}\bfv_2$, computed exactly as in Lecture 16 but without the constant first column.

def laplacian_eigenmap(W, m):
    """Coordinates phi_2, ..., phi_{m+1} and the eigenvalues lambda_1, ..., lambda_{m+1}."""
    d = W.sum(axis=1)
    L_sym = np.eye(len(d)) - W / np.sqrt(np.outer(d, d))
    lam, V = np.linalg.eigh(L_sym)                       # increasing order
    return lam[:m + 1], V[:, 1:m + 1] / np.sqrt(d)[:, None]


def knn_graph(F, k):
    A = kneighbors_graph(F, k, mode="connectivity").toarray()
    return np.maximum(A, A.T)                            # symmetrize: i ~ j if either chose the other


W = knn_graph(F, k=4)
print("connected components:", connected_components(W, directed=False)[0])

lam, Phi = laplacian_eigenmap(W, m=1)
phi = Phi[:, 0]
order = np.argsort(phi)                       # sort the shredded columns by phi_2
restored = shredded[:, order]
print("lambda_1, lambda_2 =", np.round(lam, 6))

connected components: 1
lambda_1, lambda_2 = [0.0e+00 6.6e-05]

Plot the restored photograph

image_figure([to_rgba(photo), to_rgba(restored)], ["photograph", "sorted by φ₂"], height=380).show()

Figure 2: Left: the photograph. Right: the shredded columns, sorted by the eigenmap coordinate $\phi_2$.

To score the result, we check how many consecutive columns in the restored image were also next to each other in the photograph. The few that were not are swaps of two or three columns, which are invisible.

original_position = perm[order]               # where each restored column came from
steps = np.abs(np.diff(original_position))
print(f"consecutive columns that were neighbours in the photograph: {np.mean(steps == 1):.1%}")
print("sizes of the other jumps:", np.unique(steps[steps > 1]))

consecutive columns that were neighbours in the photograph: 96.7%
sizes of the other jumps: [2 3]

Why it works

By the smoothest-embedding theorem, $\bfphi_2$ is the non-constant function on the vertices with the smallest energy $\sum_{i<j}w_{ij}(\phi_i-\phi_j)^2$, subject to the same constraints as before. A function has small energy when it changes little across every edge. Our graph is close to a path: each column is joined only to the few columns on either side of it in the photograph. The smoothest non-constant function on a path rises steadily from one end to the other. For the unnormalized Laplacian of a path on $n$ vertices it is exactly the half cosine $\cos\big(\pi(j+\tfrac12)/n\big)$, and the normalized one is very close to it. Sorting needs much less than that: any function that increases along the path puts the vertices in path order.

Figure 3 plots $\phi_2$ against the position each column had in the photograph, which the method never saw.

Plot phi_2 against the original position

j = np.arange(n_cols)
half_cosine = -np.cos(np.pi * (j + 0.5) / n_cols)
sign = np.sign(np.corrcoef(phi, half_cosine[perm])[0, 1])
scale = np.abs(phi).max() / np.abs(half_cosine).max()

fig = go.Figure([
    go.Scatter(x=j, y=sign * scale * half_cosine, mode="lines", name="half cosine",
               line=dict(color=INK, dash="dash")),
    go.Scatter(x=perm, y=phi, mode="markers", name="φ₂", marker=dict(size=5, color=BLUE),
               hovertemplate="original column %{x}<br>φ₂ = %{y:.4f}<extra></extra>"),
])
fig.update_xaxes(title="position in the original photograph")
fig.update_yaxes(title="φ₂")
style(fig, height=380).update_layout(legend=dict(x=0.02, y=0.98)).show()

Figure 3: The eigenmap coordinate $\phi_2$ of each column, against the position of that column in the original photograph. The dashed curve is a half cosine, rescaled to match.

Seriation has a second picture. In the shredded order the edges of the graph are scattered all over the adjacency matrix. In the order given by $\phi_2$, every edge joins two columns that are close in the new order, so the matrix becomes a narrow band along the diagonal.

Plot the adjacency matrix in both orders

image_figure([adjacency_image(W), adjacency_image(W[np.ix_(order, order)])],
             ["shredded order", "sorted by φ₂"], height=480).show()

Figure 4: The adjacency matrix of the $4$-nearest-neighbour graph, with each edge drawn as a small square. Left: rows and columns in the shredded order. Right: the same matrix with rows and columns sorted by $\phi_2$.

Left or right?

If $\bfphi_2$ is an eigenmap, so is $-\bfphi_2$, since a Laplacian eigenmap is only determined up to rotation and reflection, and in one dimension that just means the sign. Sorting by $-\bfphi_2$ gives the mirror image of the photograph. This is not a flaw of the method: the shredded columns carry no information about which way is left, so no method can tell a photograph from its mirror image. Our cat came out the right way round by luck of the sign; with default_rng(3) it comes out mirrored.

When it fails: a brick wall

The method relies on the graph being close to a path, which in turn relies on far-apart columns looking different. A brick wall breaks that assumption. Its pattern repeats, so a column of one brick looks very much like the matching column of a brick far away, and the nearest-neighbour graph joins them.

wall = data.brick() / 255                     # (512, 512), greyscale
wall_perm = rng.permutation(wall.shape[1])
wall_shredded = wall[:, wall_perm]

W_wall = knn_graph(wall_shredded.T, k=4)
_, Phi_wall = laplacian_eigenmap(W_wall, m=1)
wall_order = np.argsort(Phi_wall[:, 0])
wall_restored = wall_shredded[:, wall_order]

steps = np.abs(np.diff(wall_perm[wall_order]))
print(f"consecutive columns that were neighbours in the photograph: {np.mean(steps == 1):.1%}")

consecutive columns that were neighbours in the photograph: 16.2%

Plot the brick wall

image_figure([to_rgba(wall), to_rgba(wall_shredded), to_rgba(wall_restored)],
             ["photograph", "shredded", "sorted by φ₂"], height=340).show()

Figure 5: The same procedure on a photograph of a brick wall. Left: the photograph. Middle: shredded. Right: sorted by $\phi_2$.

The result is not noise: columns that were neighbours in the photograph are still joined in the graph, so pieces of brick texture reappear. But the graph also has many edges between distant columns, so it is no longer close to a path, and the smoothest function on it has no reason to increase along the original order.

--- title: "Unshredding an Image with One Eigenvector" subtitle: "Math 123: Mathematical Aspects of Data Analysis — Fall 2026" date: 2026-11-04 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 --- Cut a photograph into single columns of pixels, shuffle the columns, and hand them back: can we put them in order again? Neighbouring columns of a photograph look almost the same, so the columns are points lying along a curve in a high-dimensional space, and the order we want is position along that curve. A one-dimensional [Laplacian eigenmap](../../lectures/lecture-17.qmd#def-laplacian-eigenmap) gives each column a single number, and sorting the columns by that number puts them back in order. Using an eigenvector to put things in order, instead of to cluster or to plot them, is called *seriation*. ```{python} import base64 import io import numpy as np import plotly.graph_objects as go from PIL import Image from scipy.ndimage import binary_dilation from scipy.sparse.csgraph import connected_components from skimage import data from sklearn.neighbors import kneighbors_graph ``` ```{python} #| code-fold: true #| code-summary: "Figure style and image helpers" INK = "#7c7c78" # muted colors that read in light and dark mode GRID = "rgba(128, 128, 128, 0.28)" BLUE, ORANGE = "#2a78d6", "#eb6834" def style(fig, height): 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), height=height, margin=dict(l=60, r=20, t=30, b=50), hovermode="closest") fig.update_xaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) fig.update_yaxes(gridcolor=GRID, linecolor=GRID, zerolinecolor=GRID) return fig def to_rgba(image): """An H x W (grey) or H x W x 3 (colour) image in [0, 1] as opaque RGBA bytes.""" rgb = np.repeat(image[..., None], 3, axis=2) if image.ndim == 2 else image alpha = np.ones(rgb.shape[:2] + (1,)) return (np.clip(np.concatenate([rgb, alpha], axis=2), 0, 1) * 255).astype(np.uint8) def image_uri(rgba): """Encode RGBA bytes as a data URI: JPEG for opaque photographs, PNG otherwise. An image is far smaller than the same picture as a heatmap of numbers. """ buffer = io.BytesIO() if rgba[..., 3].min() == 255: Image.fromarray(rgba[..., :3], "RGB").save(buffer, "JPEG", quality=90) kind = "jpeg" else: Image.fromarray(rgba, "RGBA").save(buffer, "PNG") kind = "png" return f"data:image/{kind};base64," + base64.b64encode(buffer.getvalue()).decode() def image_figure(images, labels, height, gap=16): """Show RGBA images of equal height side by side, with a label under each.""" fig, left, centres = go.Figure(), 0, [] for image in images: fig.add_trace(go.Image(source=image_uri(image), x0=left, hoverinfo="skip")) centres.append(left + (image.shape[1] - 1) / 2) left += image.shape[1] + gap fig.update_xaxes(tickvals=centres, ticktext=labels, showgrid=False, zeroline=False, showline=False, ticks="") fig.update_yaxes(showticklabels=False, showgrid=False, zeroline=False, showline=False) return style(fig, height) def adjacency_image(W): """Draw each edge of W as a small ink-coloured square on a faint grey background.""" rgba = np.empty(W.shape + (4,), dtype=np.uint8) rgba[:] = (128, 128, 128, 30) rgba[binary_dilation(W > 0, np.ones((3, 3)))] = (124, 124, 120, 255) return rgba ``` ## The shredded photograph {#sec-data} We use the cat from `skimage.data`, a $300\times451$ colour photograph. Each column is $300$ pixels with three colour channels, so it is a vector in $\Rbb^{900}$, and the image is $451$ points in $\Rbb^{900}$. We shuffle the columns with a random permutation. From here on the method only sees `shredded`; we keep `perm` only to check the answer at the end. ```{python} photo = data.chelsea() / 255 # (300, 451, 3), brightness in [0, 1] n_rows, n_cols, _ = photo.shape rng = np.random.default_rng(0) perm = rng.permutation(n_cols) # column j of shredded is column perm[j] of photo shredded = photo[:, perm] F = shredded.transpose(1, 0, 2).reshape(n_cols, -1) # one row per column of pixels print(F.shape) ``` ```{python} #| label: fig-shredded #| fig-cap: "Left: the photograph. Right: the same $451$ columns of pixels in a random order." #| fig-alt: "Two images side by side. The left is a close-up photograph of a tabby cat's face. The right is the same image with its columns shuffled, which looks like vertical streaks of brown, black and cream with no recognizable shape." #| code-fold: true #| code-summary: "Plot the photograph and the shredded version" image_figure([to_rgba(photo), to_rgba(shredded)], ["photograph", "shredded"], height=380).show() ``` ## The graph and its eigenmap {#sec-eigenmap} We join each column to its $4$ nearest neighbours in $\Rbb^{900}$, in Euclidean distance, and keep an edge if either end chose the other. The weights are all $1$. The one-dimensional eigenmap is the single coordinate $\bfphi_2=\bfD^{-1/2}\bfv_2$, computed exactly as in [Lecture 16](../../lectures/lecture-16.qmd#alg-normalized-spectral-clustering) but without the constant first column. ```{python} def laplacian_eigenmap(W, m): """Coordinates phi_2, ..., phi_{m+1} and the eigenvalues lambda_1, ..., lambda_{m+1}.""" d = W.sum(axis=1) L_sym = np.eye(len(d)) - W / np.sqrt(np.outer(d, d)) lam, V = np.linalg.eigh(L_sym) # increasing order return lam[:m + 1], V[:, 1:m + 1] / np.sqrt(d)[:, None] def knn_graph(F, k): A = kneighbors_graph(F, k, mode="connectivity").toarray() return np.maximum(A, A.T) # symmetrize: i ~ j if either chose the other W = knn_graph(F, k=4) print("connected components:", connected_components(W, directed=False)[0]) lam, Phi = laplacian_eigenmap(W, m=1) phi = Phi[:, 0] order = np.argsort(phi) # sort the shredded columns by phi_2 restored = shredded[:, order] print("lambda_1, lambda_2 =", np.round(lam, 6)) ``` ```{python} #| label: fig-restored #| fig-cap: "Left: the photograph. Right: the shredded columns, sorted by the eigenmap coordinate $\\phi_2$." #| fig-alt: "Two images side by side, both showing the same close-up of a tabby cat's face. The right one, rebuilt from the shuffled columns, is indistinguishable from the original on the left." #| code-fold: true #| code-summary: "Plot the restored photograph" image_figure([to_rgba(photo), to_rgba(restored)], ["photograph", "sorted by φ₂"], height=380).show() ``` To score the result, we check how many consecutive columns in the restored image were also next to each other in the photograph. The few that were not are swaps of two or three columns, which are invisible. ```{python} original_position = perm[order] # where each restored column came from steps = np.abs(np.diff(original_position)) print(f"consecutive columns that were neighbours in the photograph: {np.mean(steps == 1):.1%}") print("sizes of the other jumps:", np.unique(steps[steps > 1])) ``` ## Why it works {#sec-why} By [the smoothest-embedding theorem](../../lectures/lecture-17.qmd#thm-laplacian-eigenmap), $\bfphi_2$ is the non-constant function on the vertices with the smallest energy $\sum_{i<j}w_{ij}(\phi_i-\phi_j)^2$, subject to the same constraints as before. A function has small energy when it changes little across every edge. Our graph is close to a *path*: each column is joined only to the few columns on either side of it in the photograph. The smoothest non-constant function on a path rises steadily from one end to the other. For the unnormalized Laplacian of a path on $n$ vertices it is exactly the half cosine $\cos\big(\pi(j+\tfrac12)/n\big)$, and the normalized one is very close to it. Sorting needs much less than that: any function that increases along the path puts the vertices in path order. @fig-phi-vs-position plots $\phi_2$ against the position each column had in the photograph, which the method never saw. ```{python} #| label: fig-phi-vs-position #| fig-cap: "The eigenmap coordinate $\\phi_2$ of each column, against the position of that column in the original photograph. The dashed curve is a half cosine, rescaled to match." #| fig-alt: "A scatter plot of 451 points with original column position from 0 to 450 on the horizontal axis and phi 2 on the vertical axis. The points form a single smooth S-shaped curve that increases from left to right, lying almost on a dashed half-cosine curve." #| code-fold: true #| code-summary: "Plot phi_2 against the original position" j = np.arange(n_cols) half_cosine = -np.cos(np.pi * (j + 0.5) / n_cols) sign = np.sign(np.corrcoef(phi, half_cosine[perm])[0, 1]) scale = np.abs(phi).max() / np.abs(half_cosine).max() fig = go.Figure([ go.Scatter(x=j, y=sign * scale * half_cosine, mode="lines", name="half cosine", line=dict(color=INK, dash="dash")), go.Scatter(x=perm, y=phi, mode="markers", name="φ₂", marker=dict(size=5, color=BLUE), hovertemplate="original column %{x}<br>φ₂ = %{y:.4f}<extra></extra>"), ]) fig.update_xaxes(title="position in the original photograph") fig.update_yaxes(title="φ₂") style(fig, height=380).update_layout(legend=dict(x=0.02, y=0.98)).show() ``` Seriation has a second picture. In the shredded order the edges of the graph are scattered all over the adjacency matrix. In the order given by $\phi_2$, every edge joins two columns that are close in the new order, so the matrix becomes a narrow band along the diagonal. ```{python} #| label: fig-adjacency #| fig-cap: "The adjacency matrix of the $4$-nearest-neighbour graph, with each edge drawn as a small square. Left: rows and columns in the shredded order. Right: the same matrix with rows and columns sorted by $\\phi_2$." #| fig-alt: "Two square images. The left one is covered in scattered dots with no pattern. The right one has all its dots in a thin band along the diagonal from top left to bottom right." #| code-fold: true #| code-summary: "Plot the adjacency matrix in both orders" image_figure([adjacency_image(W), adjacency_image(W[np.ix_(order, order)])], ["shredded order", "sorted by φ₂"], height=480).show() ``` ### Left or right? If $\bfphi_2$ is an eigenmap, so is $-\bfphi_2$, since a [Laplacian eigenmap](../../lectures/lecture-17.qmd#def-laplacian-eigenmap) is only determined up to rotation and reflection, and in one dimension that just means the sign. Sorting by $-\bfphi_2$ gives the mirror image of the photograph. This is not a flaw of the method: the shredded columns carry no information about which way is left, so no method can tell a photograph from its mirror image. Our cat came out the right way round by luck of the sign; with `default_rng(3)` it comes out mirrored. ## When it fails: a brick wall {#sec-bricks} The method relies on the graph being close to a path, which in turn relies on far-apart columns looking different. A brick wall breaks that assumption. Its pattern repeats, so a column of one brick looks very much like the matching column of a brick far away, and the nearest-neighbour graph joins them. ```{python} wall = data.brick() / 255 # (512, 512), greyscale wall_perm = rng.permutation(wall.shape[1]) wall_shredded = wall[:, wall_perm] W_wall = knn_graph(wall_shredded.T, k=4) _, Phi_wall = laplacian_eigenmap(W_wall, m=1) wall_order = np.argsort(Phi_wall[:, 0]) wall_restored = wall_shredded[:, wall_order] steps = np.abs(np.diff(wall_perm[wall_order])) print(f"consecutive columns that were neighbours in the photograph: {np.mean(steps == 1):.1%}") ``` ```{python} #| label: fig-bricks #| fig-cap: "The same procedure on a photograph of a brick wall. Left: the photograph. Middle: shredded. Right: sorted by $\\phi_2$." #| fig-alt: "Three greyscale images. The left shows a wall of bricks laid in a staggered pattern. The middle is the same image with shuffled columns, which looks like noise. The right, the attempted reconstruction, has pieces of brick texture, but its mortar lines run diagonally and break off abruptly, and it does not look like the original wall." #| code-fold: true #| code-summary: "Plot the brick wall" image_figure([to_rgba(wall), to_rgba(wall_shredded), to_rgba(wall_restored)], ["photograph", "shredded", "sorted by φ₂"], height=340).show() ``` The result is not noise: columns that were neighbours in the photograph are still joined in the graph, so pieces of brick texture reappear. But the graph also has many edges between distant columns, so it is no longer close to a path, and the smoothest function on it has no reason to increase along the original order.