The Sphere

Layer normalization via root mean square (RMS) is used by every model. It divides each token's residual-stream vector by its L2 norm before passing it to the next layer. After this operation, every token vector satisfies . The token cloud lives on the unit sphere:

On the sphere, the natural measure of similarity between two points is their inner product:

where is the angle between the two vectors. The inner product is the cosine similarity when both vectors are unit length.

The range is :

• : the two tokens point in the same direction, so their representations are identical up to the normalization.
• : the two tokens are orthogonal, so their representations share no common direction.
• : the two tokens point in opposite directions, so they have maximally anti-correlated representations.

The Histogram

For a layer with token particles, there are distinct pairs. Compute every pairwise inner product and put them in a histogram. The x-axis runs from to . The y-axis is density and integrates to 1.

Reading the histogram across layers:

Early layers build Gaussian about 0. In high dimensions, random vectors on the sphere concentrate near the equator. By concentration of measure, most pairwise inner products between randomly-oriented high-dimensional unit vectors land close to 0.

Clustering in progress: values growing at 1. As layers push similar tokens together, a subset of pairs achieves high inner product. The histogram develops a second peak migrating toward 1 while the main mass stays near 0. Each such peak is a cluster: the pairs at high inner product are tokens that have been pulled together.

Full collapse: all values at 1. When all tokens have merged into a single point on the sphere, every pairwise inner product equals 1. The histogram is a single spike at 1. This is what the theory's long-run prediction (consensus / a single Dirac mass) looks like empirically.

Mass-Near-1

Tracking the full histogram across all layers, models, and prompts produces many plots. We We simplify the histogram into a scalar: Mass-near-1 is the fraction of pairs with inner product above 0.9.

This is a direct, threshold-based read on "how much of the token cloud has clustered so far." A layer where mass-near-1 goes from 0 to 0.4 in one step is a layer where 40% of pairs snapped into high agreement at once, signaling a merge event.

The pair extraction is one helper. It takes the upper triangle of G (each pair once, no self-pairs):

def pairwise_inner_products_from_gram(G: np.ndarray) -> np.ndarray:
    """Upper-triangle pairwise cosine similarities from a pre-computed Gram matrix."""
    n   = G.shape[0]
    idx = np.triu_indices(n, k=1)   # k=1 drops the diagonal (self-pairs = 1.0)
    return G[idx]

The histogram, the summary scalars, and mass-near-1 are then four lines in the per-layer loop:

ips                  = pairwise_inner_products_from_gram(G)
lr["ip_mean"]        = float(ips.mean())
lr["ip_std"]         = float(ips.std())
lr["ip_histogram"]   = np.histogram(ips, bins=50, range=(-1, 1))[0].tolist()
lr["ip_mass_near_1"] = float((ips > 0.9).mean())

The last line is the formula above, verbatim: (ips > 0.9) is the indicator , and .mean() over the upper triangle is the division by . The histogram uses 50 fixed bins on as a fixed range, so the bin edges are comparable across every layer, model, and prompt without renormalizing. This is the array the Figure-1 replication plots: results["layers"][li]["ip_histogram"] is fed directly to the bar chart, one panel per layer.

There is no separate code path for the histogram and for mass-near-1. They read the same ips array. The scalar is just the histogram with everything above the 0.9 bin edge summed.

Take the activation matrix (n tokens, d dimensions per token). Compute its singular values . Convert them to a probability distribution by normalizing, then take the exponential of the Shannon entropy:

This is the effective rank of Roy and Vetterli (2007).

Understanding Effective Rank

Entropy measures how spread-out a distribution is. A distribution concentrated on one value has entropy 0; a uniform distribution over values has entropy . Exponentiating maps entropy back to a "number of effective components."

Applied to the singular value spectrum:

All weight on one direction (): , entropy , effective rank . The cloud is essentially 1D. This is a line through the origin, meaning all tokens are nearly identical.

Spectrum flat (): entropy , effective rank . The cloud uses all directions equally.

In between: smooth interpolation. A cloud where a few directions dominate gets an effective rank near that count, not some artifact of whatever threshold you happened to pick.

The key advantage is that it's smooth and weights directions by their share of the spectrum. A direction carrying a negligible fraction of the total barely contributes, even though it would count as a "nonzero singular value" in a threshold-based count.

Reading Across Layers

Effective rank tends to start high at early layers (many directions contribute roughly equally) and fall as clustering progresses (token representations align, the spectrum concentrates in fewer directions). The endpoint depends on architecture:

• Full collapse: effective rank -> 1. All tokens have essentially the same representation.
• Partial collapse: effective rank stabilizes above 1. Multiple clusters means multiple persistent directions.
• No collapse: effective rank stays near . Tokens remain spread out.

Effective rank tracks the full distribution of the spectrum, not just the high-agreement pairs like the inner products from before. A model with two perfectly separated clusters of equal size will show mass-near-1 = 0 (no cross-cluster pairs near 1) but effective rank ≈ 2 (two principal directions). They're measuring different things and can complement each other from their local and global views.

The Degeneracy Gate

When effective rank drops low enough, the token cloud is nearly a point-mass. At this point:

• CKA between consecutive layers goes to 1. Any two representations that are both near-point-masses have nearly identical Gram matrices, so CKA can't distinguish them. A CKA value of 0.99 at a degenerate layer means "both layers are collapsed," not "the geometry is stable in a nontrivial way."
• Nearest-neighbor assignment becomes noise. When all tokens are nearly identical, which token is nearest to which is determined by floating-point rounding at the scale of . NN-stability values computed here are not interpretable.
• Spectral cluster counts are meaningless. The Laplacian eigengap method finds clusters in a graph; a near-point-mass graph has no structure to find.

Rather than report these unhelpful quantities in the analysis, we gate them on effective rank and suppress them below the threshold.

The metric itself:

from scipy.linalg import svdvals

def effective_rank_from_raw(activations: torch.Tensor) -> float:
    """
    Spectral entropy of singular values.

    SVD runs on raw activations, not L2-normed ones. L2 normalization sets
    every token's norm to 1, collapsing the inter-token scale variation that
    the singular values measure - svdvals(normed) gives a different answer.
    Named _from_raw to make the contract explicit at call sites.
    """
    sv      = svdvals(activations.numpy())
    sv      = sv[sv > 1e-10]              # drop numerical zeros
    sv_norm = sv / sv.sum()              # normalize to a probability distribution
    entropy = -np.sum(sv_norm * np.log(sv_norm + 1e-12))
    return float(np.exp(entropy))

squared vs. unsquared. The concept box above wrote , i.e. normalizing the variances/eigenvalues of . The code normalizes the singular values directly, sv / sv.sum(). These are not the same function; the squared form puts more weight on the leading directions and reports a lower effective rank. The implemented version is the original Roy–Vetterli (2007) definition. If you want the variance-weighted version, square sv before the normalize line. The findings below are computed with the unsquared form as shown.

Two rank variants are now computed per layer. The raw rank drives every degeneracy gate; the normed rank (directional spread on the sphere, independent of residual norm growth) is kept for reporting:

# Raw: captures scale + directional collapse; used for all gates.
# Normed: directional spread on the sphere only; for reporting.
lr["effective_rank"]        = effective_rank_from_raw(activations)
lr["effective_rank_normed"] = effective_rank_from_normed(normed)

The gate is then applied at the call site. Effective rank is computed first, then used to decide whether CKA is even meaningful at this layer:

from core.config import DEGENERATE_RANK_THRESHOLD   # = 2

# CKA vs previous layer. suppressed when the cloud is a near-point-mass.
# Centering then produces noise-dominated vectors and the Frobenius norms
# collapse to near-zero, so the ratio is numerically meaningless.
if prev_normed is not None and lr["effective_rank"] >= DEGENERATE_RANK_THRESHOLD:
    lr["cka_prev"] = linear_cka(normed, prev_normed)
else:
    lr["cka_prev"] = float("nan")

The gate is a single shared constant, DEGENERATE_RANK_THRESHOLD = 2, used identically by CKA, NN-stability, and energy-drop suppression. It was previously split (CKA gated at >= 3.0, NN at >= 2.0); unifying it at 2 matches the geometric claim in the concept box above: at rank 2 the cloud is still genuinely 2-D on the sphere (two clusters separated by a great circle), so CKA and NN remain meaningful, while rank ≈ 1 is the point-mass regime where they become float-noise. Degenerate layers drop out of the plateau analysis rather than contributing false near-1 values.

The Clustering Problem

You have token vectors on the sphere at a given layer. You want to know: how many distinct groups are there, and which tokens belong to which group?

The answer must handle:

• Unknown number of groups which may vary per layer
• Groups that might not be spherical clusters of equal size
• Tokens that are in between or transitioning which don't belong to a group

DBSCAN (Density-Based Spatial Clustering of Applications with Noise) defines clusters as dense regions separated by sparse regions, and is capable of working well despite the challenges faced.

If tokens are genuinely clustering, there should be regions of high local density (the cluster cores) separated by low-density gaps. A point is a "core point" if it has at least minPts neighbors within distance . Clusters grow by linking neighboring core points. Points reachable from a core point but not core themselves are border points. Points in no core point's neighborhood are noise (label = ).

This handles variable-shape clusters naturally. Density is local, so elongated clusters work without issue. The noise label doesn't force boundary tokens into clusters.

The problem with vanilla DBSCAN is that it uses a flat density threshold everywhere. This fails when clusters have different densities. A threshold that captures a sparse cluster will merge two nearby dense clusters; a threshold that separates those dense clusters will dissolve the sparse one.

HDBSCAN: The Hierarchical Extension

HDBSCAN (McInnes and Healy, 2017) fixes this by building a full density hierarchy.

Step 1: Mutual reachability distance. For each point, compute its core distance , the distance to its -th nearest neighbor. The mutual reachability distance between points and is:

In sparse regions, distances are replaced by the larger core distances, so sparse points are pushed apart. In dense regions, points are already closer than their core distances, so nothing changes. The result is a distance metric robust to density variation.

Step 2: Minimum spanning tree (MST). Build the MST on the full mutual reachability distance matrix. This is the skeleton of the density structure.

Step 3: Build the condensed cluster tree. Simulate removing MST edges in order of increasing mutual reachability distance (equivalently: decreasing density threshold), tracking which clusters split off. This produces a dendrogram of cluster births and deaths.

Step 4: Extract stable clusters. Each branch has a "persistence," or how long it survives as a distinct cluster as the threshold rises. The algorithm selects the most persistent set of non-overlapping clusters. Branches that split off and die quickly are absorbed back as noise.

The result is cluster assignments for tokens in stable dense regions, a noise label () for tokens in low-density regions or transitions, and a cluster count that was discovered, not specified.

The Noise Label

During metastability, some tokens have committed to a cluster and others haven't yet. The noise label captures this. A layer where 30% of tokens are labeled noise may be a signal that the clustering process is mid-transition. Tracking the noise fraction across layers is itself a useful diagnostic.

On the Sphere

We apply HDBSCAN to the sphere-projected token vectors using cosine distance (). This is because the tokens live on and the paper's dynamics are stated in terms of angular relationships.

Cosine distance is precomputed from the normed activations, then handed to HDBSCAN as a precomputed metric:

from sklearn.metrics import pairwise_distances
import hdbscan

# 1 − ⟨x_i, x_j⟩ on the sphere; clip removes tiny negative round-off
cos_dist = np.clip(pairwise_distances(normed, metric="cosine"), 0, None)

hdb        = hdbscan.HDBSCAN(min_cluster_size=2, metric="precomputed")
hdb_labels = hdb.fit_predict(cos_dist.astype(np.float64))   # −1 = noise

# Cluster count excludes the noise label
n_clusters = len(set(hdb_labels)) - (1 if -1 in hdb_labels else 0)

Three things to read off this:

• metric="precomputed" means HDBSCAN never sees the vectors, only the cosine-distance matrix. The angular geometry is the only thing it clusters on, consistent with the paper's inner-product dynamics.
• min_cluster_size=2 is deliberately permissive: a "cluster" can be a single merged pair, which is what an early merge event looks like.
• n_clusters subtracts 1 if and only if the noise label is present. The noise tokens (-1) are counted separately, never folded into a cluster. The noise fraction, (hdb_labels == -1).mean() , is tracked per layer as the mid-transition signal described above.

The labels array is what feeds every downstream consumer: the trajectory tracker that records merges across layers (Group B), the per-pair semantic/artifact tagging, and the multi-scale nesting analysis all take hdb_labels as their starting point.

The cluster-count sweep runs both HDBSCAN and a k-means silhouette search, which makes the k-means failure modes directly observable in the tests:

def test_antipodal_kmeans_best_k_is_2(self, antipodal_normed):
    # Two tight clusters at opposite poles. Silhouette peaks at k=2.
    result = cluster_count_sweep(antipodal_normed)
    assert result["kmeans"]["best_k"] == 2

def test_collapsed_agglomerative_count_is_1(self, collapsed_normed):
    # Fully collapsed cloud: every cosine distance ~1e-6, below all thresholds.
    # Must return a single cluster.
    result = cluster_count_sweep(collapsed_normed)
    assert result["agglomerative"][_MID_THRESH] == 1

K-means recovers the right answer only when you already know to look for and the clusters happen to be the clean, equal-sized, antipodal case it assumes. HDBSCAN gets the count without that assumption, which is why it, and not k-means, is the authoritative membership source throughout.

The Object Being Compared

At layer you have an activation matrix (sphere-projected). At layer you have , same shape. You want one number saying how similar these two representations are geometrically.

The natural object is the representational similarity matrix (RSM): the Gram matrix , whose entry is . Two representations are considered the same geometry if their RSMs agree, regardless of how the underlying coordinates are oriented. This is exactly the inner-product matrix from Group A. CKA is built on top of the same object.

HSIC

The Hilbert-Schmidt Independence Criterion measures whether two sets of features are statistically related. For linear kernels and , the unnormalized linear HSIC reduces to

after centering. F is the Frobenius norm. This gives us the alignment between the two RSMs. it is large when pairs of tokens that are close in are also close in .

Normalization → CKA

HSIC on its own scales with the magnitudes of and , so it is not comparable across layers. CKA normalizes it by the self-similarities:

The result is in : 1 means the two RSMs are identical up to an orthogonal transform and isotropic scaling. 0 means they are orthogonal.

Invariances

This is the property that makes CKA the right plateau instrument:

• Invariant to orthogonal transformations. If layer is layer rotated, CKA = 1. The residual stream rotating between layers does not register as change.
• Invariant to isotropic scaling. Uniformly scaling all activations does not change CKA.
• Non-invariant to arbitrary invertible linear maps. A measure invariant to any invertible linear transform would call almost any two full-rank representations identical, which is not helpful for detecting persistence. CKA only calls representations the same if their relational geometry is the same.

So a CKA plateau near 1 means: across these layers, the model is rotating and rescaling the token cloud but not reorganizing which tokens are close to which. That is precisely the metastable picture we are seeking. We want to detect when the clusters are fixed and the dynamics are idling.

Reading It Across Layers

• Flat near 1: consecutive layers are representationally identical: a metastable plateau.
• Sharp single-step drop: the representation reorganized between two layers: the end of a plateau, typically a merge event.
• Low and noisy: no stable structure.

The metric is the normalized HSIC formula:

def linear_cka(X: np.ndarray, Y: np.ndarray) -> float:
    """
    CKA(X, Y) = ||Y^T X||_F^2 / (||X^T X||_F * ||Y^T Y||_F)
    X, Y : (n_tokens, d). already L2-normed; both mean-centered internally.
    """
    X = X - X.mean(axis=0, keepdims=True)
    Y = Y - Y.mean(axis=0, keepdims=True)
    # ||Y^T X||_F^2 = tr(X^T Y Y^T X) = ||X.T @ Y||_F^2
    YtX       = Y.T @ X                          # (d, d)
    numerator = float(np.sum(YtX ** 2))
    XtX_norm  = float(np.linalg.norm(X.T @ X, "fro"))
    YtY_norm  = float(np.linalg.norm(Y.T @ Y, "fro"))
    denom     = XtX_norm * YtY_norm
    if denom < 1e-12:
        return float("nan")
    return float(np.clip(numerator / denom, 0.0, 1.0))

In the per-layer loop it is called only when the layer is non-degenerate (the Group A gate), so collapsed layers return nan rather than a spurious 1.0:

from core.config import DEGENERATE_RANK_THRESHOLD   # = 2

if prev_normed is not None and lr["effective_rank"] >= DEGENERATE_RANK_THRESHOLD:
    lr["cka_prev"] = linear_cka(normed, prev_normed)
else:
    lr["cka_prev"] = float("nan")

A CKA plateau is then just a flat run of this series. The plateau detector (see B.3) is run on the CKA values with a tight tolerance, because CKA near 1 should be very flat inside a real metastable window:

cka_pairs    = [(r["layer"], r["cka_prev"]) for r in layers
                if not np.isnan(r["cka_prev"])]   # drop layer 0 + degenerate layers
cka_series   = [v for _, v in cka_pairs]
plateaus_cka = detect_plateaus(cka_series, window=2, tol=0.02

The end of a plateau is flagged separately as the sharpest single-step drop, where consecutive-layer CKA falls the most is the clearest signal that a metastable window has ended:

diffs    = np.diff(valid_vs)
drop_pos = int(np.argmin(diffs))
if diffs[drop_pos] < -0.05:
    severity = "SHARP" if diffs[drop_pos] < -0.15 else "MILD"
    # marks L → L+1 as the plateau-ending reorganization

The Measure

For each token at a layer, its nearest neighbor is . This is the token it points most directly towards on the sphere. NN-stability between layer and is the fraction of tokens whose nearest neighbor index is identical:

1.0 means every token kept the same nearest neighbor, signifying a perfectly locked plateau. 0.0 means every token's nearest neighbor changed, meaning the cloud is still reorganizing.

Why It Is Not Redundant With CKA

CKA is invariant to rotation and measures aggregate relational geometry. Two failure modes it cannot see:

• Partner swapping within a stable geometry. Imagine two tightly bound pairs of tokens, A↔B and C↔D, that across a layer become A↔C and B↔D. The overall pairwise-similarity structure can stay nearly identical (same distances, different identities), so CKA barely moves, but every token changed its nearest neighbor: NN-stability drops to 0. This is genuine reorganization that a metastability claim must rule out.
• The discrete vs. continuous gap. CKA can drift smoothly while NN identities snap. NN-stability is a step function on identities; it marks the exact layers where the discrete cluster graph rewires.

The NN Functional Graph

Because each token has exactly one nearest neighbor, the NN map defines a functional graph (every node has out-degree 1). The cycles of this graph are the stable atoms: a mutual pair is a 2-cycle, and longer cycles are tightly bound groups. At plateau layers, these cycles are the operational definition of a "stable cluster" used in the token-membership reporting. They are tagged as semantic (members are distinct, structurally parallel tokens) or duplicate (members are the same token string). The semantic-vs-duplicate split keeps positional repeats from being counted as meaningful clusters.

The Degeneracy Problem

When the cloud collapses to a near-point-mass, every token is nearly equidistant from every other, and which one is "nearest" is decided by floating-point noise at the level. NN-stability computed there is meaningless. So, like CKA, it is suppressed below the rank threshold.

The nearest neighbor of each token is one masked argmax on the Gram matrix:

def nearest_neighbor_indices(G: np.ndarray) -> np.ndarray:
    """nn[i] = argmax_{j≠i} G[i, j], from a pre-computed Gram matrix."""
    G_masked = G.copy()
    np.fill_diagonal(G_masked, -np.inf)   # exclude self (G[i,i] = 1)
    return np.argmax(G_masked, axis=1).astype(np.int32)

Stability is the fraction unchanged from the previous layer, computed inline in the loop (layer 0 has no predecessor, so it is None):

nn = nearest_neighbor_indices(G)
lr["nn_indices"] = nn.tolist()
if prev_nn is not None:
    lr["nn_stability"] = float(np.mean(nn == prev_nn))   # the indicator-mean formula
else:
    lr["nn_stability"] = None
prev_nn = nn

When NN-stability plateaus are extracted, degenerate layers are dropped first:

nn_stab_defined = [
    (i, v) for i, v in enumerate(nn_stab_vals)
    if v is not None and layers[i]["effective_rank"] >= DEGENERATE_RANK_THRESHOLD
]
nn_series = [v for _, v in nn_stab_defined]
nn_plateaus = detect_plateaus(nn_series, window=2, tol=0.02)

The semantic/duplicate split is read off the NN functional graph at plateau layers: cycles whose members are distinct token strings are semantic; cycles of identical strings are positional duplicates. That tagging is what lets the report distinguish a real structural cluster from the same word appearing three times.

The Flatness Criterion

A window of layers is a plateau if the signal's relative span across it stays below a tolerance. For a segment of values , the criterion is

Using relative span (dividing by the mean) makes the criterion scale-appropriate per signal. The detector starts from a minimum width, then greedily extends the window as long as the flatness criterion still holds, and records (start, end, mean).

Per-Signal Tolerances

Because the signals live on different scales, each gets its own tolerance:

The Operational Plateau Set vs. the Joint Criterion

This is worth stating precisely, because the intent and the implementation differ.

The authoritative plateau set stored on each run (plateau_layers, step P1-7) is computed from mass-near-1 alone, with tol 0.10. That single signal defines the windows that downstream analyses (semantic coherence at plateau layers, the two-timescale numerator) actually use.

Separately, a multi-signal vote is computed: for each layer, count how many of {mass, rank, spectral-k, hdbscan-k, fiedler, CKA} have a plateau covering it, and flag layers covered by ≥2. This stricter "joint" set is used in the flagged-anomalies cross-check, not as the operational plateau set.

The detector is one greedy scan, identical for every signal:

def detect_plateaus(values: list, window: int = 2, tol: float = 0.05) -> list:
    """Contiguous windows where relative span < tol. Returns (start, end, mean)."""
    plateaus = []
    n = len(values)
    i = 0
    while i < n - window:
        segment = values[i:i + window + 1]
        span    = max(segment) - min(segment)
        ref     = abs(np.mean(segment)) + 1e-8
        if span / ref < tol:                       # the flatness criterion
            j = i + window
            while j < n - 1:                       # greedily extend
                extended = values[i:j + 2]
                if (max(extended) - min(extended)) / (abs(np.mean(extended)) + 1e-8) < tol:
                    j += 1
                else:
                    break
            plateaus.append((i, j, float(np.mean(values[i:j + 1]))))
            i = j + 1
        else:
            i += 1
    return plateaus

The operational plateau set is mass-near-1 only (this is the plateau_layers every downstream consumer reads):

# Post-loop (P1-7): the authoritative plateau set
mass1    = [r["ip_mass_near_1"] for r in results["layers"]]
plateaus = detect_plateaus(mass1, window=2, tol=0.10)
plateau_layer_set = set()
for s, e, _ in plateaus:
    for l in range(s, e + 1):
        plateau_layer_set.add(l)
results["plateau_layers"] = sorted(plateau_layer_set)

The multi-signal vote is computed separately, as a stricter cross-check (note it is not what plateau_layers uses):

layer_count = Counter()
for group in [plateaus_mass, plateaus_rank, plateaus_spk, plateaus_hdb, plateaus_fied]:
    for s, e, _ in group:
        for l in range(s, e + 1):
            layer_count[l] += 1
for s, e, _ in plateaus_cka:                       # CKA indexed by position → remap to layer
    for idx in range(s, e + 1):
        layer_count[cka_pairs[idx][0]] += 1
multi = [l for l, c in sorted(layer_count.items()) if c >= 2]   # ≥2 signals agree

Prompt sensitivity turns the onset layer into the weight-level/content-driven verdict:

onset_vals = [first_plateau_onset(run) for run in runs_of_this_model]
sd = float(np.std(onset_vals))
classification = "weight-level" if sd < 2.0 else "content-driven"

The Matching Problem

At layer HDBSCAN found some clusters; at layer it found some clusters. Which layer- cluster is which layer- cluster? They are the same cluster if they contain mostly the same tokens. The overlap measure is Jaccard:

computed over token membership, with HDBSCAN noise tokens (label ) excluded.

Hungarian Assignment

Given the full matrix of Jaccard overlaps between every layer- cluster and every layer- cluster, the optimal one-to-one matching is the assignment that maximizes total overlap. That is the Hungarian algorithm (linear_sum_assignment), run on the negated overlap matrix because the routine minimizes cost. Matches below a minimum Jaccard (0.1) are discarded as too weak to be continuations.

So the description is both "Jaccard overlap" and "Hungarian matching": Hungarian is the assignment procedure, Jaccard is the cost it optimizes.

Births, Deaths, Merges

After the optimal matching:

• Match: a layer- cluster paired to a layer- cluster -> the cluster continued.
• Birth: a layer- cluster with no match in -> new structure appeared.
• Death: a layer- cluster with no match in -> structure dissolved.
• Merge: an unmatched layer- cluster that nonetheless has strong overlap with a layer- cluster that is already matched to a different layer- cluster. Two layer- clusters mapping into one layer- cluster express the slow timescale.

Trajectory Chains

Matched clusters are chained across layers into trajectories. Each trajectory is a sequence of (layer, cluster_id) pairs with a birth layer, a death (or end) layer, and a lifespan. Trajectory lifespans are the per-cluster version of the plateau width: a long-lived trajectory is a cluster that survived a long metastable stretch before merging.

Why This, Not Spectral-k, for ALBERT-xlarge

Group A explained that ALBERT-xlarge's Laplacian has a dominant zero mode that pins the spectral cluster count at regardless of real structure. A merge-detection method based on spectral- drops would therefore report zero merges there, manifesting an artifact. Trajectory tracking works entirely from HDBSCAN token membership, which has no such failure mode, so its merge counts are the authoritative ones for ALBERT-xlarge. Any nMerges = 0 in a spectral-sourced column for that model is the spectral artifact, not an absence of merges.

The overlap matrix is Jaccard over membership sets, noise excluded:

def _jaccard_overlap_matrix(labels_a, labels_b):
    ids_a = sorted(set(labels_a) - {-1})
    ids_b = sorted(set(labels_b) - {-1})
    sets_a = {c: set(np.where(labels_a == c)[0]) for c in ids_a}
    sets_b = {c: set(np.where(labels_b == c)[0]) for c in ids_b}
    overlap = np.zeros((len(ids_a), len(ids_b)))
    for i, ca in enumerate(ids_a):
        for j, cb in enumerate(ids_b):
            inter = len(sets_a[ca] & sets_b[cb])
            union = len(sets_a[ca] | sets_b[cb])
            overlap[i, j] = inter / union if union else 0.0
    return overlap, ids_a, ids_b

The matching is Hungarian on the negated overlap, then merges are the unmatched-prev-into-matched-curr case:

from scipy.optimize import linear_sum_assignment

# maximum-weight assignment → minimize negated overlap (padded to square)
cost = np.zeros((size, size)); cost[:n_prev, :n_curr] = -overlap
row_ind, col_ind = linear_sum_assignment(cost)

matches = []
for r, c in zip(row_ind, col_ind):
    if r < n_prev and c < n_curr and overlap[r, c] >= min_jaccard:   # min_jaccard = 0.1
        matches.append((ids_prev[r], ids_curr[c], float(overlap[r, c])))

# merge: an unmatched prev cluster overlapping a curr cluster that is already matched
for up in unmatched_prev:
    best_j = int(np.argmax(overlap[ids_prev.index(up), :]))
    target = ids_curr[best_j]
    if overlap[..., best_j] >= min_jaccard and target in matched_curr:
        # record (prev clusters that fed `target`, target) as one merge event

Trajectories are chains of matched (layer, cluster_id) tips, advanced one transition at a time. The run-level summary is what the report prints:

"summary": {
    "total_births":  total_births,
    "total_deaths":  total_deaths,
    "total_merges":  total_merges,
    "max_alive":     max_alive,
    "n_trajectories": len(traj_info),
    "mean_lifespan":  float(np.mean([t["lifespan"] for t in traj_info])),
    "max_lifespan":   max(t["lifespan"] for t in traj_info),
}

The Measure

An attention matrix row is a probability distribution: and (softmax output). Its Shannon entropy is

Averaging over the rows gives a per-head scalar. The bounds are interpretable:

• Uniform attention ( for all ): . The token attends equally to everything. Maximum entropy. No routing structure.
• Hard attention (each token attends to exactly one token): . Minimum entropy. Fully concentrated routing.

So entropy runs from (no structure) down to (maximally concentrated), and watching it fall across layers tells you attention is sharpening.

The Connection to β

In the paper's dynamics, attention weights are . The inverse temperature controls sharpness: at attention is uniform (maximum entropy); as attention becomes a hard argmax (zero entropy), and low entropy corresponds to high effective . Therefor, attention entropy is a monotone proxy for the effective the model is operating at.

This is why entropy is the right opener for Group C. The metastability conjecture in the paper is a large- phenomenon. If measured attention entropy is high and flat (effective near zero), there is no reason to expect metastability and the geometric clustering would need a different explanation. If entropy falls with depth, the model is moving into the regime where the theory's prediction applies, and the geometric clustering from Group A has a mechanism behind it.

Why per-head, then averaged

Entropy is computed per head because certain heads may concentrate, some may stay diffuse, and an average over heads would conceal this information. Group C reports the per-head values and uses the mean only as a summary.

The entire instrument is one vectorized function over all heads at once:

def attention_entropy(attn_matrix: torch.Tensor) -> np.ndarray:
    """
    Shannon entropy of each attention row, averaged over tokens.
    attn_matrix : (n_heads, n_tokens, n_tokens)
    Returns (n_heads,), mean entropy per head.
    """
    attn              = attn_matrix.numpy()
    log_attn          = np.log(attn + 1e-12)                 # 1e-12 guards log(0)
    entropy_per_token = -(attn * log_attn).sum(axis=-1)      # (n_heads, n_tokens)
    return entropy_per_token.mean(axis=-1)                    # (n_heads,)

-(attn * log_attn).sum(axis=-1) is for every row of every head; .mean(axis=-1) averages over the rows. The two reference points are pinned by tests: uniform attention returns and identity (each token attends only to itself) returns , confirming the scale runs from no-structure to fully-concentrated as described.

Why Make It Doubly Stochastic

Raw attention is row-stochastic (each row sums to 1), but not column-stochastic, as some tokens receive far more attention than others. The paper's Section 3.3 connects the idealized dynamics to a doubly stochastic object (both rows and columns sum to 1), which is the symmetric, balanced form a gradient flow would produce. The gap between raw attention and its doubly stochastic version is itself a diagnostic ("row/col balance" below): zero means attention is already balanced, large means it is far from the idealized form.

Sinkhorn–Knopp produces the doubly stochastic matrix by alternately normalizing rows and columns until both converge. It is the standard, provably convergent way to get there.

The Fiedler Value

Treat the doubly stochastic matrix as a weighted graph (after symmetrizing, ). The graph Laplacian encodes its connectivity. The eigenvalues of the normalized Laplacian start at (always, it is the trivial constant mode) and increase. The second one, , is the Fiedler value (algebraic connectivity):

• : the graph is nearly disconnected. It has near-separate components. In attention terms, tokens are routed into clusters that barely attend across to each other. This is the metastable signature.
• large: the graph is well-connected. Every token can reach every other through attention. Tokens are mixing.

The number of near-zero eigenvalues equals the number of near-disconnected components, which is why a related eigengap reading gives a cluster count (below).

Why This Is the Group C Payoff

The Group C question is whether attention knows about the geometric clusters. The Fiedler value makes that testable directly: if low-Fiedler layers (cluster-separated attention) coincide with the layers where Group A finds geometric clusters and Group B finds plateaus, then attention and geometry are telling the same story. The Spearman cross-check below quantifies exactly this coincidence.

Cluster Count From the Eigengap

The same spectrum gives a cluster count: for a -cluster structure, the largest eigenvalues of sit near 1 with a drop below them. The largest gap in the descending eigenvalue sequence locates without a hard threshold.

Sinkhorn–Knopp is alternating row/column normalization to convergence, batched across all heads:

def sinkhorn_normalize_batched(A, max_iter=SINKHORN_MAX_ITER, tol=SINKHORN_TOL):
    """A: (n_heads, n, n) raw attention → doubly stochastic per head."""
    P = np.clip(A.astype(np.float64), 1e-12, None)
    for _ in range(max_iter):
        P_prev = P.copy()
        P /= P.sum(axis=2, keepdims=True)   # row-normalise all heads
        P /= P.sum(axis=1, keepdims=True)   # col-normalise all heads
        if np.abs(P - P_prev).max() < tol:
            break
    return P

The Fiedler value is of the normalized Laplacian of the symmetrized matrix:

def fiedler_value(P: np.ndarray) -> float:
    """λ₂ of the normalised Laplacian. λ₂≈0 → cluster-separated; large → mixing."""
    P_sym       = (P + P.T) / 2
    L           = laplacian(P_sym, normed=True)
    k           = min(3, P.shape[0] - 1)
    eigenvalues = eigh(L, eigvals_only=True, subset_by_index=[0, k - 1])
    return float(eigenvalues[1]) if len(eigenvalues) > 1 else 0.0

The cluster count reads the eigengap, falling back to a hard threshold only when there is no real gap:

def sinkhorn_cluster_count(P, min_gap_ratio: float = 0.1) -> int:
    # k largest eigenvalues near 1; largest gap in the descending sequence = k.
    # Fallback to hard >0.5 count when largest gap < min_gap_ratio*(λmax−λmin),
    # i.e. on near-uniform (post-collapse) matrices with no genuine structure.
    ...

The per-layer summary records the mean Fiedler, the per-head Fiedler list, the eigengap cluster count, and the row/col balance (distance of raw attention from doubly stochastic):

result = {
    "fiedler_mean":                float(np.mean(fiedler_vals)),
    "fiedler_per_head":            fiedler_vals,
    "sinkhorn_cluster_count_mean": float(np.mean(cluster_counts)),
    "row_col_balance_mean":        float(np.mean(row_col_balance)),
}

Caveat carried from the plateau detector. Fiedler plateaus are found with detect_plateaus(fiedler, tol=0.05), and that routine measures relative flatness (span / (|mean| + 1e-8)). The cluster-separated state is exactly Fiedler ≈ 0, where the denominator collapses and the relative criterion can fail to register a genuinely flat-near-zero window.

Does low Fiedler co-occur with high cluster count?

The single most direct answer to "does attention know about the clusters" is the correlation between the attention-graph Fiedler value and the geometric (HDBSCAN) cluster count, layer by layer. If attention is implementing the clustering, low Fiedler (separated attention) should co-occur with high cluster count (multi-cluster geometry), which is a negative correlation.

fied_hdb_pairs = [
    (r["sinkhorn"]["fiedler_mean"],
     r["clustering"]["hdbscan"]["n_clusters"])
    for r in layers
    if "sinkhorn" in r and not np.isnan(r["clustering"]["hdbscan"]["n_clusters"])
]
rho, pval = spearmanr([p[0] for p in fied_hdb_pairs],
                      [p[1] for p in fied_hdb_pairs])
# rho < -0.4 → attention Fiedler tracks geometric cluster count (signal)
# rho ≈ 0    → attention and geometry are telling different stories

A Spearman is the threshold the report uses to call this an interpretable signal; would mean the attention graph and the token geometry are decoupled, which would undercut the mechanistic claim even though both individually show clustering.

The Classification

For each head, collect its Fiedler value at every active-phase layer, take the mean, and bin it:

• CLUSTER: mean Fiedler < 0.3 -> consistently routes tokens into separated clusters.
• MIXED: 0.3–0.7 -> variable behavior across layers.
• MIXING: > 0.7 -> consistently allows free mixing.

The Active-Phase Restriction

Classification is restricted to layers where effective rank ≥ 10. The reason is a saturation artifact: once tokens collapse to a near-point-mass (rank below ~10), there is only one cluster, the doubly stochastic attention matrix is nearly uniform, and the Laplacian has no gap, so every head's Fiedler value tends to 1. Including those layers would pull every head's mean toward 1.0 and label everything MIXING regardless of its real role. Excluding them keeps the classification meaningful. (Note this rank-10 threshold is specific to this analysis and separate from the rank-2 degeneracy gate used for CKA/NN: it is a larger cutoff because Fiedler saturation sets in well before full point-mass collapse. It is currently hardcoded in the profiler and would be cleaner in config.)

STABLE vs VARIABLE Across Prompts

A head that lands in the same class for every prompt is STABLE. Its routing role is a property of the weights, input-independent. A head whose class changes with the prompt is VARIABLE, meaning it is content driven. The cross-prompt consistency table reports this per head per model

The Causal-Mask Confound

This is the central caution of Group C. GPT-2 attention is causally masked (lower-triangular). Sinkhorn-normalizing and symmetrizing a lower-triangular matrix forces a low-connectivity graph regardless of content, which would manufacture low Fiedler, and therefore a "100% STABLE-CLUSTER" reading across every prompt, as an artifact of the mask rather than a fact about the weights.

The main analysis computes raw per-head Fiedler, and for causal models also the mask-only baseline and each head's deviation from it:

fiedler_vals = [fiedler_value(P_all[h]) for h in range(n_heads)]
result["fiedler_per_head"] = fiedler_vals

# Fix 3a: causal-mask baseline subtraction
if is_causal:
    baseline = causal_fiedler_baseline(n)               # Fiedler of uniform-causal attn
    result["fiedler_baseline"]           = baseline
    result["fiedler_per_head_deviation"] = [round(f - baseline, 6) for f in fiedler_vals]

causal_fiedler_baseline is the content-free reference the mask produces on its own:

def causal_fiedler_baseline(n: int) -> float:
    A_base = _uniform_causal_attention(n)   # uniform within the lower triangle
    P_base = sinkhorn_normalize(A_base)
    return fiedler_value(P_base)            # what the mask alone forces

The BERT control applies an artificial causal mask to a bidirectional model to test whether masking alone collapses heads to CLUSTER:

if apply_causal_control:
    causal_mask = np.tril(np.ones((n, n)))
    attn_masked = attn * causal_mask[None, :, :]            # zero upper triangle
    attn_masked /= attn_masked.sum(axis=2, keepdims=True)   # renormalise rows
    P_ctrl = sinkhorn_normalize_batched(attn_masked)
    result["fiedler_causal_control_per_head"] = [fiedler_value(P_ctrl[h]) for h in range(n_heads)]

The classifier then chooses the signal based on model type, deviation for causal models, raw Fiedler otherwise, and restricts to the active phase:

# in _per_head_fiedler_profile, active phase = effective_rank >= 10
# causal models classify on fiedler_per_head_deviation; others on raw fiedler_per_head
classification = ("CLUSTER" if mean < 0.3 else
                  "MIXED"   if mean < 0.7 else
                  "MIXING")

The Energy

For tokens on the sphere with Gram matrix , the interaction energy at inverse temperature is

Read the structure: every term is large when tokens and point the same way (inner product near ) and small when they point apart. So is large when tokens are aligned and small when they are spread out. Maximal clustering, every token identical, maximizes it. The analytical reference values make this concrete: for a fully collapsed cloud , for two antipodal clusters , for a uniform spread , and , so collapsed > antipodal > uniform.

The Prediction: Monotone Increasing

The idealized dynamics ascend this energy. As depth increases and tokens cluster, should rise at every layer. Plotted against layer index, it should be a non-decreasing curve. This is the gradient-ascent picture: the model is climbing toward consensus, and the energy is the height it has climbed.

Violations of Energy Monotonicity

A violation is a layer where decreases: . Mechanistically this is sharp. For the energy to fall, some pairs must have their terms shrink, which means those tokens moved apart. That is a locally repulsive violation of energy.

The pure-attraction dynamics in the paper cannot do this. Attention only pulls tokens together; it has no repulsive term. Under the idealized flow, the energy is a monotonic Lyapunov function. So an energy drop is not a small quantitative deviation; it is qualitative evidence of a force the theory does not contain. Something in the trained layer pushed two tokens apart. Group D's job is to detect that, and D.2's job is to find the tokens it happened to.

The energy is computed for all in one vectorized pass over the cached Gram matrix:

def interaction_energies_batched(G: np.ndarray, beta_values: list) -> dict:
    """E_beta = (1 / 2β n²) Σ_ij exp(β ⟨x_i, x_j⟩), for every β at once."""
    n        = G.shape[0]
    betas    = np.asarray(beta_values, dtype=np.float64)   # (B,)
    exp_G    = np.exp(betas[:, None, None] * G[None])       # (B, n, n)
    sums     = exp_G.sum(axis=(1, 2))                       # (B,)
    energies = sums / (2.0 * betas * n * n)                 # (B,)
    return {float(beta): float(e) for beta, e in zip(beta_values, energies)}

In the loop it is one line per layer:

lr["energies"] = interaction_energies_batched(G, beta_values)   # {β: E_β}

The violation test is applied to one -series across layers:

ENERGY_VIOLATION_REL_TOL = 1e-3   # a drop counts only if > 0.1% of |E_prev|

def energy_violation_severity(energies, rel_tol=ENERGY_VIOLATION_REL_TOL):
    arr   = np.array(energies, dtype=np.float64)
    diffs = np.diff(arr)
    ref   = np.maximum(np.abs(arr[:-1]), 1e-12)   # |E| at the preceding layer
    rel_drop = -diffs / ref                       # positive = energy fell
    viol_mask   = rel_drop > rel_tol
    return dict(
        violation_layers = [i + 1 for i, v in enumerate(viol_mask) if v],
        n_violations     = int(viol_mask.sum()),
        max_severity     = float(rel_drop[viol_mask].max()) if viol_mask.any() else 0.0,
        sum_severity     = float(rel_drop[viol_mask].sum()),
        total_rel_change = float((arr[-1] - arr[0]) / max(abs(arr[0]), 1e-12)),
    )

total_rel_change is the net climb from first to last layer (the energy does rise overall but isn't monotone); violation_layers are the steps where it dropped against the trend. The reporting layer runs this for every also checks whether the violation layers coincide with the merge events to see whether the repulsive moves happen at merges or independently of them.

Per-Pair Contribution

The energy is a sum over pairs, so the change in energy between layer and decomposes exactly into per-pair contributions:

A pair with contributed to the energy drop. Its tokens moved apart (lower inner product at than at ). The most-negative are the pairs most responsible for the violation. Sorting all pairs by ascending and taking the top few gives the localization.

The Structural-Token Question

Each token in a top pair is checked against a set of structural tokens (i.e., [CLS], [SEP], <s>, </s> , padding, single-character punctuation, so on). The interpretation could be:

• Repulsion driven by structural tokens: the energy drop is likely a positional/structural effect due to special tokens being pushed out of clusters as bookkeeping, not semantic content.
• Repulsion driven by semantic tokens: content-bearing tokens being pushed apart is a stronger signal, meaning the repulsive force acts on meaning-geometry.

Degeneracy Suppression

Once the cloud collapses to a near-point-mass, the inner products are all , the terms are nearly equal, and the per-pair deltas are floating-point noise. Energy "violations" detected there are not real, so violation layers in the degenerate regime are suppressed and the pair localization is not reported for them.

The public wrappers are thin (one normalizes raw tensors, one takes the pre-normed arrays from the loop) and both call a shared core:

def energy_drop_pairs_from_normed(normed_before, normed_after, beta, top_k=10):
    """Token pairs (i,j) sorted by Δ ascending (most negative first), where
       Δ = [exp(β⟨x_i,x_j⟩_after) − exp(β⟨x_i,x_j⟩_before)] / (2β n²)."""
    return _energy_drop_pairs_core(normed_before, normed_after, beta, top_k)

The core is the pairwise delta, upper triangle only, sorted to surface the most repulsive pairs:

def _energy_drop_pairs_core(normed_before, normed_after, beta, top_k):
    n = normed_before.shape[0]
    if n < 2:
        return []
    G_before = normed_before @ normed_before.T
    G_after  = normed_after  @ normed_after.T
    delta = (np.exp(beta * G_after) - np.exp(beta * G_before)) / (2 * beta * n * n)
    iu = np.triu_indices(n, k=1)                       # i < j only
    pairs = [(int(i), int(j), float(delta[i, j])) for i, j in zip(*iu)]
    pairs.sort(key=lambda t: t[2])                     # ascending → most negative first
    return pairs[:top_k]

In the loop, this is computed only at violation layers and only when the layer is non-degenerate. Using prev_normed:

# prev_normed still points to layer L here (Fix 8 deferred its update)
lr["energy_drop_pairs"] = {
    beta: energy_drop_pairs_from_normed(prev_normed, normed, beta, top_k=10)
    for beta in beta_values
} if (is_violation and lr["effective_rank"] >= DEGENERATE_RANK_THRESHOLD) else {}

The reporting layer maps the top pairs back to token strings and annotates structural ones:

SPECIAL_TOKENS = {"[CLS]","[SEP]","<s>","</s>","<pad>","[PAD]","<|endoftext|>","Ġ","▁"}
PUNCT_CHARS    = set(".,!?;:'\"-–—()[]{}…/\\")
# each top pair printed as:  'tok_i'[FLAG] ↔ 'tok_j'[FLAG]  with δ
# [FLAG] ∈ {[CLS], [SEP], [PUNCT]} marks structural/special-token repulsion

The Definition

We compute everything from one real-prompt run:

• formation_time = the layer index of the first stable multi-cluster plateau (the first plateau where the mean cluster count is ≥ 1.5, i.e. genuinely more than one cluster). This is when the fast process finishes, signaling that clusters exist and are holding.
• merge_time = the number of layers from that first plateau until collapse (or, if no collapse, until the structure decays). This is how long the slow process runs.
• separation = merge_time / formation_time, both measured on the same run.

Thresholds: ≥ 3.0 = CONFIRMED, 1.5–3.0 = WEAK, < 1.5 = NO SEPARATION.

The cluster count is read from HDBSCA, so the multi-cluster test does not depend on the spectral method that degenerates for ALBERT-xlarge. The repeated-token control is kept, but only as a standalone "collapse speed diagnostic."

The ratio is defined as:

def _single_trajectory_separation(layers_data) -> dict:
    """Two-timescale separation from ONE real-prompt run.
       formation_time = first stable multi-cluster plateau onset
       merge_time     = layers from that plateau to collapse/decay
       separation     = merge_time / formation_time"""
    # cluster count per layer: HDBSCAN preferred, spectral k fallback
    k_series = [layer_cluster_count(l) for l in layers_data]

    # multi-cluster plateaus: detect_plateaus on k, keep windows with mean k >= 1.5
    plats = detect_plateaus(k_series, window=2, tol=0.5)
    multi = [(s, e) for (s, e, mean_k) in plats if mean_k >= 1.5]
    if not multi:
        return {"separation": 0.0, "verdict": "NO SEPARATION"}

    formation_time = multi[0][0]                       # onset of first multi-cluster plateau
    collapse_layer = first_collapse_or_decay(k_series) # k → 1 (or structure decays)
    merge_time     = collapse_layer - formation_time

    separation = merge_time / formation_time if formation_time > 0 else 0.0
    verdict = ("CONFIRMED"      if separation >= 3.0 else
               "WEAK"           if separation >= 1.5 else
               "NO SEPARATION")
    return {"formation_time": formation_time, "merge_time": merge_time,
            "separation": separation, "verdict": verdict}

The repeated-tokens runs remain in the report under a COLLAPSE SPEED DIAGNOSTICS heading, explicitly marked as not a metastability test.

Singular Values vs. Eigenvalues

For the value matrix (square, , for all three architectures' full projection), two spectra carry different information:

• Singular values give the magnitudes of how much stretching is seen in each orthogonal direction on the linear map. They track total spectral energy but discard sign.
• Eigenvalues carry sign and phase. A real eigenvalue with is an attractive direction (the map pulls along it); is repulsive; a complex eigenvalue is rotational. This sign structure is exactly the attractive/repulsive tension that Group D's energy violations implied must exist somewhere, and 's eigenvalues are where Phase 2 (next post) looks for it.

Spectral Radius as Collapse Rate

The spectral radius is the dominant per-layer amplification. For an iterated linear map, structure along the dominant eigendirection grows like over layers. So controls how many layers it takes to collapse: a back-of-envelope estimate for ALBERT-xlarge, , is iterations to reach high mass, which is why it has not collapsed at 48 and matches its observed slow dynamics.

The Theorem 6.1 Test

Theorem 6.1 predicts higher → faster convergence. The decisive comparison is ALBERT-base () vs. ALBERT-xlarge (). If the theorem governs, the higher-dimensional xlarge should collapse faster. It collapses far slower: ALBERT-base collapses by ~24 iterations, ALBERT-xlarge has not collapsed at 48. The dimension prediction is backwards for this pair. The variable that does order them correctly is , lower spectral radius, slower collapse, so the spectral radius, not , is the governing quantity. That is a direct falsification of the Theorem 6.1 prediction in trained models.

Why This Resolves E.1's Contradiction

Order the models by rather than by depth and the contradiction dissolves. BERT-base has (the fast-clean regime) despite being only 12 layers. GPT-2-medium has (24 layers, no separation); GPT-2-large has (36 layers, separation). The "depth threshold between 24 and 36 layers" is really a spectral-radius threshold that happens to fall between medium and large within the GPT-2 family. Depth is a within-family proxy that breaks the moment BERT (shallow but low-) is included. The ratio table looked contradictory because it was sorted by the wrong axis.

The value matrix is pulled per architecture (per-layer for BERT/GPT-2, one shared matrix for ALBERT), then the sign-aware eigenspectrum is computed:

# extraction differs by architecture; GPT-2 slices V out of the fused c_attn
if "gpt2" in model_name:
    for i, block in enumerate(model.h):
        d = block.attn.c_attn.weight.shape[1]
        v = block.attn.c_attn.weight[:, 2*d//3:].detach().cpu().float().numpy()  # last third = V
        v_matrices.append((f"layer_{i}", v))
# (bert: layer.attention.self.value.weight per layer;
#  albert: one shared attn.value.weight)

for name, V in v_matrices:
    sv = svdvals(V)                                # magnitudes (sign-blind)
    if V.shape[0] == V.shape[1]:
        eigs       = np.linalg.eigvals(V)          # complex — carries sign + phase
        real_arr   = np.real(eigs)
        is_complex = np.abs(np.imag(eigs)) > 0.01 * (np.abs(real_arr) + 1e-8)
        spec_radius = float(np.abs(eigs).max())    # ← collapse-rate predictor
        frac_pos    = float((real_arr > 0).mean()) # attractive directions
        frac_neg    = float((real_arr < 0).mean()) # repulsive directions
        frac_complex = float(is_complex.mean())    # rotational directions

spec_radius, eig_frac_pos_real, eig_frac_neg_real, and eig_frac_complex are stored per layer in v_eigenspectrum.json for cross-referencing against the Group A/B plateau and collapse locations.

Methodological seam (Phase 2 refinement). This Phase 1 instrument eigendecomposes alone. Phase 2 argues the physically meaningful operator is the composed OV circuit , the actual residual-stream-to-residual-stream map, and that eigendecomposing in isolation omits the output projection. Phase 2's extract_ov_circuit replaces this. The Phase 1 spectral-radius numbers are directionally robust (the ordering and the Theorem 6.1 falsification hold), but the precise operator is refined there; the Phase 1 figures should be read as the -only proxy, not the final OV-circuit values.