Winning the GPU MODE eigh competition

World's fastest batched symmetric eigendecomposition on a B200

17 October 2026

GPU MODE ran a competition on eigh: computing the eigenvalues and eigenvectors of a batch of symmetric matrices on an NVIDIA B200. The baseline is torch.linalg.eigh, which runs the 13 benchmark cases in 55.0 ms (geometric mean). My final submission runs them in 6.60 ms on the official grader, and it won.

I had a bit of kernel experience from earlier GPU MODE competitions but zero background in the math behind this problem. I entered for fun and because I wanted to try doing research with coding agents. The agents ran mostly on their own, and my job ended up being high-level direction.

The code is on GitHub: submission.py is the cleaned-up version (3,254 lines) and messy.py is the competition submission as it was scored.

Over 2 weeks 6.60 ms ended up being the fastest time I could achieve.

Eigendecomposition

What we are computing

A square matrix A sends a vector v to another vector Av. For most v the result points in a new direction. For a few special directions it stays on the same line, and Av = λv. Those directions are the eigenvectors, and the scale factors λ are the eigenvalues. When A is symmetric, all the eigenvalues are real and the eigenvectors can always be chosen orthogonal to each other. Stacking the unit eigenvectors as the columns of a matrix Q gives

A = Q diag(L) Qᵀ,   QᵀQ = I

You can read this right to left: Qᵀ rotates into the eigenvector basis, diag(L) scales each axis, and Q rotates back. The 2×2 case already shows the whole idea:

The input is a batch × n × n FP32 tensor. The output is Q and L in the torch.linalg.eigh convention, with eigenvectors in the columns of Q and eigenvalues in ascending order. Batched small eigendecompositions show up in optimizers that factor each layer’s second-moment matrix on every step, and the competition’s shapes come from that use case.

Eigenvectors are not unique. You can flip their sign, and when eigenvalues repeat, the eigenvectors can rotate freely inside the shared subspace. So the checker cannot compare against a reference answer. Instead it checks three identities in FP64, scaled by n·eps (FP32 epsilon), the first two relative to the matrix norm:

This was an important detail of the competition. The checker only asks that the answer is good enough, and “good enough” is loose enough that FP16 is allowed inside the solver. A lot of the speedups below come from doing work in FP16 and repairing the error once at the end.

These are the 13 ranked cases:

The score is a geometric mean, so making the tiny n32 case 30 percent faster helps exactly as much as making n2048 30 percent faster. That means no case can be ignored.

How torch does it

If we profile torch.linalg.eigh at n512, we see cuSOLVER running the classic three-stage algorithm:

1. Tridiagonalization. Turn A into a tridiagonal matrix T (nonzeros only on the diagonal and right next to it) with the same eigenvalues.
2. Tridiagonal eigensolver. Solve T, which is much easier. cuSOLVER uses divide and conquer here.
3. Back-transform. Map the eigenvectors of T back to eigenvectors of A.

The first stage is where most of the time is spent, about 60 percent of it, in kernels with names like cuds_symv_lo_direct and gemvx, which are matrix-vector products. It uses Householder reflectors. A reflector is a mirror that zeros out everything below the subdiagonal of one column. Applying it to both sides of the matrix keeps the matrix symmetric and keeps its eigenvalues. After n − 2 reflectors, one per column, the matrix is tridiagonal. Roughly, in NumPy:

for i in range(n - 2):
    v, tau = householder(A[i+1:, i])        # mirror that zeros A[i+2:, i]
    p = tau * (A[i+1:, i+1:] @ v)
    w = p - 0.5 * tau * (v @ p) * v
    A[i+1:, i+1:] -= np.outer(v, w) + np.outer(w, v)   # rank-2 update

The catch is what each reflector costs. Applying it to the rest of the matrix needs a matrix-vector product (symv) between the remaining part of the matrix and the reflector’s vector, and the next column can’t start until that’s done. So a matrix of size n needs n − 2 matrix-vector products, one after the other. Tensor cores don’t help here. A matrix-vector product does 2 flops per element it loads, so in FP16 that is 1 flop per byte, and a B200 needs on the order of 250 flops per byte to keep its FP16 tensor cores busy. Within this algorithm, the levers are making each step faster and overlapping other work with the chain. Most of the first half of this post is about those two levers.

Picking the algorithm

The textbook advice for large matrices is a two-stage reduction: first reduce the dense matrix to a band (a few diagonals wide) using big matrix multiplies, then reduce the band to tridiagonal with a “bulge chase”. The first stage can use tensor cores, so it moves most of the work off the slow matrix-vector path.

Early on, three agents built three different approaches in parallel: two-stage, the classic one-stage reduction, and a spectral divide and conquer that I’ll come back to in the dead ends. At n512, one-stage did the reduction and back-transform in 85.5 ms, and two-stage took 190 ms. Of those 190 ms, 126 were the bulge chase alone. The chase itself could be made fast: with a narrow band it later got from 48 to 13.7 ms. But there is a trade-off. A narrow band makes the chase cheap but leaves many band reflectors to apply in the back-transform, and a wide band does the opposite, so no band width is cheap on both ends. The best joint budget was 34.8 ms against 29.9 ms for one-stage.

The reason the textbook advice doesn’t apply is the batch. Here we have 640 matrices at n512, so every SM already has its own matrix to work on, and the one-stage version’s slow matrix-vector product is not leaving the GPU idle. For the small shapes the first kernels were Jacobi solvers, which rotate pairs of rows and columns until the matrix is diagonal. Those were quick early wins at n32, n176 and n352, but for the big shapes the one-stage route won and Jacobi was removed.

What the reduction is up against

Every design decision in the reduction kernel comes back to the same chain: 510 dependent matrix-vector products per matrix at n512. For each one, the questions are how many bytes it reads, where they come from, and how many loads are in flight while it waits. Some arithmetic sets the stage.

Bytes. Each column’s matrix-vector product reads the whole remaining matrix, so over a full reduction the kernel reads about n³/3 values. At n512 and batch 640 in FP32 that is about 115 GB, which at B200’s 8 TB/s nominal, about 7 achievable, is a 16 ms floor. The first kernel took 55 ms, and the profiler put DRAM throughput at 15 percent (part of the traffic was served from L2). So it wasn’t moving too many bytes, it was waiting on them. That says the first lever is more loads in flight, and the second is fewer bytes, since a smaller working set also stays in cache.

Cache. The remaining matrix at n512 starts at 512 KB in FP16. With one CTA per SM on 148 SMs, the machine is working on about 76 MB at once, which is already in the range where the B200’s L2 (126 MB on paper, but split across two dies) starts missing. This came up in a surprising way: a version with two CTAs per SM, each holding a packed 262 KB triangle of the matrix, had twice as many warps hiding latency and was slower, because the L2 hit rate fell from 77 to 20 percent. It is also why FP8 kept getting proposed. An FP8 copy would have been half the size and fit, but FP8 doesn’t have the precision, and that attempt died four times.

SMs. The chain is serial within one matrix, so the only way to add parallelism to one matrix is to give it more SMs, and the only way to split it is across a thread-block cluster that has to join at least once per column. The batch decides. At n512 there are 640 matrices for 148 SMs, so one CTA per matrix fills the machine 4.3 times over. At n1024 there are 60, so one CTA per matrix would leave more than half the GPU idle, and each matrix is spread over a 2-CTA cluster. At n2048 there are 8, so the cluster is 8 CTAs. Each step up pays a cluster-wide join per column, which is what the barrier section below is about.

The panel width, 32 columns, was a guess that held. An agent swept it near the end and every big shape was already at its optimum.

Making the reduction fast

The one-stage reduction kernel is LAPACK’s latrd in CUDA: it reduces a panel of 32 columns at a time, then updates the rest of the matrix with one matrix multiply per panel. The first working version spent 54.5 ms per n512 batch in this kernel. One agent then took it to 24.2 ms in five steps, and I like this sequence because every step is a classic GPU lesson.

A profile confirmed the arithmetic above: the kernel used about 15 percent of memory bandwidth, and most of its time went into waiting on shared memory. It also flagged two inefficiencies, bank conflicts and uncoalesced global loads.

Padding shared memory (54.5 → 48.8 ms). Shared memory is split into 32 banks, and if several threads of a warp hit different words in the same bank, their reads happen one after another. The FP16 panel was stored with a row stride of 32 elements, 64 bytes, so reading down a column hit only two banks, a 16-way conflict. Padding the stride to 33 spreads consecutive rows across all 32 banks.

Reading a row instead of a column (48.8 → 45.9 ms). One step read column i of a row-major matrix, which is a strided access. But the matrix is symmetric, so column i is the same as row i, and reading the row is a coalesced access. It also used fewer registers.

Two rows per warp (45.9 → 36.0 ms). In the matrix-vector product, each warp handled one row at a time: load, wait for the load, reduce. Handling two rows at once means twice as many loads are in flight before the warp has to wait. The final form of the loop, with NRW rows per warp:

float acc[NRW] = {};
#pragma unroll
for (int t = 0; t < TCW; ++t)          // column pairs of the tile
  #pragma unroll
  for (int q = 0; q < NRW; ++q)        // NRW independent rows -> NRW loads in flight
    acc[q] = fhfma2_acc(tile[(warp + q * WARPS) * TC + 32 * t + lane], vh[t], acc[q]);

Four rows per warp (36.0 → 30.3 ms). Same idea. The agent swept 2, 3, 4, 6 and 8. Four was best; eight ran out of registers and spilled.

An FP16 copy of the matrix (30.3 → 24.2 ms). The matrix-vector product reads the whole remaining matrix for every column. Keeping a second copy in FP16 for just this product halves the bytes read. This is safe because the first thing the pipeline does is divide each matrix by its largest entry, so every value is in a range where FP16’s 11 bits of precision are enough. bf16 would have been about 8× less accurate, and FP16’s extra precision is what made this pass the checker. The diagonal and off-diagonal of T still come out in FP32.

The kernel got a lot more work later, but this sequence was one of the biggest jumps of the whole campaign.

Against the floor: 24.2 ms is 1.5× the 16 ms FP32 bandwidth floor. But the FP16 copy halves the bytes, so the floor is now about 8 ms, and the kernel was 3× off it. By the end of the competition the n512 reduction was around 14 to 15 ms, under 2× that floor. The remaining gap is latency, not bandwidth, and we never closed it.

One instruction instead of three

Later, the matrix-vector loop was reading the FP16 copy of the matrix and an FP16 copy of the vector, and adding up the products in FP32, since summing hundreds of terms in FP16 loses too much accuracy. That meant converting each value to FP32 before the multiply. For each pair of values the loop issued six instructions: two loads, two conversions and two FMAs.

The agent working on this had actually been sent to try packed FP16 math that also accumulates in FP16. Its accuracy simulation showed that would blow the error budget at n512. While looking for another way, it found that B200 has a mixed-precision FMA, fma.rn.f32.f16: two FP16 inputs, product added into an FP32 accumulator, one instruction. There is no intrinsic for it, so it’s inline PTX:

__device__ __forceinline__ float fhfma(unsigned short a, unsigned short b, float c) {
  float r;
  asm("{.reg .b16 ha, hb;\n\t"
      "mov.b16 ha, %1;\n\t"
      "mov.b16 hb, %2;\n\t"
      "fma.rn.f32.f16 %0, ha, hb, %3;}\n"
      : "=f"(r) : "h"(a), "h"(b), "f"(c));
  return r;
}
__device__ __forceinline__ float fhfma2_acc(__half2 a, __half2 b, float c) {
  unsigned au = *(unsigned*)&a, bu = *(unsigned*)&b;
  c = fhfma(au & 0xffff, bu & 0xffff, c);
  c = fhfma(au >> 16,    bu >> 16,    c);
  return c;
}

Six instructions became three. In the SASS this is a single FHFMA per element with no unpack instructions at all, and the shared-memory instruction-queue (MIO throttle) stall cycles dropped 53 percent.

It doesn’t cost any accuracy either. Two FP16 numbers have 11 bits of precision each, so their product fits exactly in FP32’s 24 bits. The only rounding is the FP32 add, which is the same as before. n512 dense went from 38.7 to 35.7 ms, n1024 from 29.3 to 26.6 ms, and n2048 from 55.5 to 47.9 ms, with the checker’s margins unchanged.

There is a nice follow-up. The agent was later asked to build the FP16-accumulate version anyway, and it measured no faster, slightly slower at n2048. After this change the loop was waiting on memory, not on instruction issue, so deleting more instructions bought nothing.

Removing barriers

For n1024 and n2048, one matrix is split across a cluster of 2 or 8 CTAs, which can read each other’s shared memory. The kernel synchronized the cluster three times per column with cluster.sync(). At n512 the same kernel ran with a cluster of one CTA.

An agent sampled where the kernel was stalling, line by line in the compiled assembly (SASS). At n512, cluster.sync() on a cluster of one CTA compiled into a full GPU-scope memory fence plus the cluster barrier hardware, and it was 12.2 percent of the kernel’s samples, while synchronizing with nobody. The fix is a compile-time branch:

if constexpr (C == 1) {
  __syncthreads();   // every buffer is CTA-local, there is nobody to sync with
} else {
  cluster.sync();    // release fence (lowered to MEMBAR.GPU) + cluster rendezvous, three times per column
}

Same output bit for bit, and n512 dense went from 32.4 to 29.2 ms. The other n512 cases moved by the same amount.

On the real clusters we do need to synchronize, but not three times per column. The agent replaced the three cluster.sync() calls with a single handoff. Each CTA writes its partial results straight into its peers’ shared memory with st.async, which counts the bytes into the receiver’s mbarrier as they land:

// peer address of the same smem offset in CTA `rank`
__device__ unsigned mapa_u32(unsigned saddr, int rank) {
  unsigned r; asm("mapa.shared::cluster.u32 %0, %1, %2;" : "=r"(r) : "r"(saddr), "r"(rank)); return r;
}
// store into peer smem; completion is reported to the peer's mbarrier, no fence needed
__device__ void st_async_f32(unsigned dst, float v, unsigned mbar) {
  asm volatile("st.async.weak.shared::cluster.mbarrier::complete_tx::bytes.f32 [%0], %1, [%2];"
               :: "r"(dst), "f"(v), "r"(mbar) : "memory");
}

Two things make this cheaper. A cluster barrier carries a cluster-scope release fence, and on sm_100 ptxas can only lower that to a full MEMBAR.GPU; the st.async stores report completion through the mbarrier’s byte count and need no fence. And the three rendezvous per column become one: each CTA still waits once for its peers’ bytes, but no longer three times. n1024 dense went from 22.1 to 20.5 ms, and n2048 from 40.5 to 36.0 ms.

FP16 eigenvectors, fixed at the end

The back-transform multiplies the tridiagonal eigenvectors by all the reflectors from the first stage. Many reflectors can be packed into one matrix of the form I − V T Vᵀ (the “compact-WY” form), which turns the back-transform into ordinary matrix multiplies that tensor cores can run. B200 has fast FP16 tensor cores but no tensor-core path for full FP32. It does have TF32, and that was tried first: with a single TF32 pass the orthogonality residual came out at 370 to 684 against a limit of 100, and the three-pass TF32 split that reaches FP32 accuracy won only n1024, by 0.3 ms, and lost n512 by 1.5 ms.

Doing the back-transform in FP16 is fast, but the eigenvectors come out less orthogonal than the checker allows. The fix is one Newton-Schulz step:

def ns_purify_fp16(Z):
    Zh = Z.half()
    G = fp16_baddbmm(Zh.transpose(1, 2), Zh, alpha=-0.5)   # G = -0.5 Zᵀ Z, fp32 out
    G.diagonal(dim1=1, dim2=2).add_(1.5)                  # M = 1.5 I - 0.5 ZᵀZ
    return fp16_baddbmm(Zh, G.half())                     # Q = Z M

If Z is almost orthogonal, this pushes it much closer: the error after the step is roughly the square of the error before, which is why one step is enough. It is just two more FP16 matrix multiplies. Here is what it did to the orthogonality residual on the hardest input family at each size, where the checker’s limit is 100:

The result was more orthogonal than the old FP32 path (the n512 orthogonality residual went from 50.8 to 17.2), and faster too: n512 went from 51.1 to 47.5 ms, n1024 from 41.3 to 39.2 ms, and n2048 from 63.0 to 61.3 ms. The eigenvalues are never stored in FP16.

Much of the second half of the campaign repeated this pattern. When a kernel’s output was going to be rounded to FP16 by the next kernel anyway, the producer wrote FP16 directly. Each of those was worth 0.3 to 2.7 percent on a shape.

Launch overhead and CUDA graphs

The pipeline launches dozens of kernels per call, and for the smaller shapes the launch gaps were a real fraction of the time. CUDA graphs fix this: you record the sequence of launches once and replay it with one call.

Graph replay was worth 16 percent at n512, 5 percent at n1024, 17 percent at n2048 and 26 percent at n176. Graphs also made one more trick easy. At n512 the batch is split into three chunks whose graphs form one DAG, so one chunk’s tridiagonal solve runs while another chunk is still being reduced.

The final pipeline

Putting it together, the big shapes run this:

1. Prescale. Divide each matrix by its largest entry and write the FP16 copy next to the FP32 one.
2. Tridiagonalization. The latrd kernel from above, 32 columns per panel. TMA, the hardware’s bulk-copy engine, streams the FP16 matrix into shared memory in tiles without going through registers, signalling each arrival on an mbarrier. A ring of four tiles keeps the next one arriving while the current one is consumed; measured, the ring delivers about 121 bytes per SM per cycle, more than the product consumes. After each panel, a tensor-core kernel updates the rest of the matrix. The last 128 to 288 columns, where the panel machinery no longer pays off, run as an unblocked lane-per-row reduction from shared memory (latrd_tail_kernel). At n1024 and n2048 a cluster of 2 or 8 CTAs shares one matrix.
3. Tridiagonal divide and conquer. Split T in half recursively until the leaves are 32 rows, solve each leaf by bisection and inverse iteration with one warp per leaf and one lane per eigenvalue, then merge back up. A merge (merge_build_kernel) solves one scalar equation per eigenvalue, the secular equation, in FP32 with a final FP64 polish, and deflation skips the eigenvalues the merge doesn’t move.
4. Back-transform in FP16 on tensor cores.
5. One Newton-Schulz step.
6. Graph replay of all of the above.

Here is how the batch maps onto the 148 SMs for each shape:

And this is where the time goes. Even at the end, the tridiagonalization, with its chain of matrix-vector products, is most of the runtime on every shape.

The small shapes have their own routes, and the arithmetic is different there because the batch is smaller than the machine. n32 has 20 matrices, so occupancy is irrelevant and the whole solve is one kernel, syevd_block: one 128-thread block per matrix, all three stages in one launch. Warp 0 does the serial parts (the reflector and the matrix-vector product) and all 128 threads do the rank-2 update. Launching it directly measured 4.6 percent faster than graph replay, because a graph launch of a single kernel costs more than a direct launch. n176 and n352 have 40 matrices each. A cluster of 3 CTAs per matrix gives 120 CTAs, which fits in one wave of 148 SMs; 4 per matrix would be 160 and need two waves. Both run a fast FP16 route and then check every matrix’s residual against a threshold tighter than the checker’s, re-solving the failures in FP32.

The n512 “clustered” case has a route of its own. Its eigenvalues sit in two tight clusters, 170 near one value and 342 near another. The submission checks each batch for this structure, and when it’s there, it uses a solver built for exactly that shape of spectrum. If the check fails, the batch goes through the general route. This took the case from 18.8 ms to 4.68 ms, the largest single gain late in the competition. Real matrices rarely look like this, though, so the route matters a lot less outside the benchmark.

These are the final per-case times against torch:

What didn’t work

The most time went into trying to replace the classic algorithm with something more novel. None of it beat the three-stage algorithm that torch already uses, just made fast.

Spectral divide and conquer splits the spectrum in half using the matrix sign function, which you can compute with only matrix multiplies, so it maps well onto tensor cores. It passed the checker, but at n1024 it ran in 176 ms against torch’s 101 ms. A later version built on the QDWH polar decomposition measured 106 ms at n512 against 75 ms for our route at the time, because the polar step alone took 62 ms instead of the 14 ms we had budgeted for it.

The two-stage reduction was built four times, including the ELPA-style schedule, and it lost every time on the joint budget above. The chase’s serial depth, 1850 to 2000 rounds, also barely depends on the band width.

Tensor cores for the matrix-vector product looked promising in a standalone test, but in the real kernel they were 3.7 percent slower at n512 and 15 percent slower at n1024. A matrix-vector product gives the tensor core almost nothing to do, and the MMA units ran at 1.78 percent utilization.

FP8 was tried in four forms. Halving the matrix copy again would have brought the L2 working set under the cache size, but FP8’s error is larger than the gaps between eigenvalues at n512, and repairing that cost more than it saved. MRRR, another tridiagonal eigensolver, was numerically fine, but just computing its eigenvalues cost as much as our whole divide-and-conquer tail.

Some dead ends came back to life. The direct n32 solver was rejected at experiment 130 because an agent believed its orthogonalization step ran every time. It doesn’t, and the same solver shipped at experiment 241, 33 percent faster than the Jacobi kernel it replaced.

Running the agents

How the system ran

An orchestrator session ran the loop: it read the results log, decided what to work on next, wrote a mission for it, and launched a subagent in its own git worktree off main. When the subagent came back, the orchestrator checked the result and merged it if it passed. Several subagents ran at once, all sharing one B200.

The results log shows how little of this was me. Of its 836 entries, 749 were written by subagent sessions, 57 by the orchestrator session, and 30 carry no author. None were written by me directly.

The subagents came in a few flavors. Kernel agents owned a set of shapes and could rewrite anything on their path. Architecture explorers were asked for designs outside the current one, which is where the spectral and two-stage attempts came from. A precision optimizer widened the checker margins without costing speed, and a simplifier made the code shorter without costing speed. Research scouts read papers and code and came back with a ranked list of ideas.

Running without a human only works if the tooling catches the mistakes a human would have caught. Most of the scripts below exist because something went wrong without them.

Sharing one GPU

Several agents benchmarking on one GPU at the same time will corrupt each other’s timings and run the GPU out of memory. So every GPU command from every worktree went through one script, gpu_run.sh, which holds a single file lock and runs one job at a time. It also loads the grader’s exact CUDA toolchain and re-pins the GPU clocks after every job, because Nsight drops the clock pin and the next benchmark would otherwise run at half speed.

Even this had incidents. A nested lock froze the queue for six minutes until it got a re-entrancy guard. An agent wrapped the script in its own timeout, which killed jobs that were still waiting in the queue.

Knowing whether something is faster

Early on, a good change was worth 10 or 20 percent and measurement noise didn’t matter. Late in the competition most wins were under 1 percent, and the sub-millisecond cases jitter 2 to 6 percent from run to run, so a single measurement can’t resolve them. A search that keeps whichever candidate measured fastest then starts keeping lucky measurements instead of better code, and it still looks like progress from the outside.

A win had to beat the baseline in every round of a three-round alternating A/B test at pinned clocks. The small shapes later got a stricter rule and were measured with a cold L2 cache, the way the grader measures. That rule came from one n176 change that measured 4.6 percent faster in one setup and 4.8 percent slower in another, on a case too short to resolve either way.

Before anything merged into main, one script ran the full check on the merged code: all 39 correctness tests, the benchmark with the grader’s own timing function run until its standard error is 0.1 percent, a mode that reuses inputs across iterations and rechecks every one (added after an uninitialized buffer produced NaNs only in the grader’s benchmark mode).

The last piece was a per-case ratchet. Once, a small per-call environment check costing 21 µs passed every test. It made n32 16 percent slower, but in the geometric mean that was smaller than the noise. After that, any single case more than 2 percent slower than its best-ever time (4 percent for the sub-millisecond n32) failed the merge.

Remembering what failed

With hundreds of experiments and agents that start every session fresh, the same bad idea gets tried again and again. So every measured result, win or loss, went into the log with the shapes it applied to, the change in time, the mechanism, and a pointer to the evidence. A small script compiled each agent’s brief from that log, so an agent working on n1024 started with everything that had already been tried on n1024.

The subtle part was that a rejected idea is only rejected under certain conditions. An FP16 copy of the matrix for the cluster kernels was first rejected as slower. A later agent went back to it and wrote:

The prior “fp16 shadow for cluster symv” KILL … conflated TWO changes: fp16 bytes AND ld.global.cs evict-first. The evict-first was the culprit — it DEFEATS the L2 reuse this latency-bound path depends on.

With the eviction hint, n2048 took 212 ms. FP16 on its own took n2048 from 118 to 107 ms, and n1024 from 68.3 to 63.7 ms. The compiler version, the kernel generation and the accuracy budget all changed during the competition too, so every rejection in the log recorded the conditions it was measured under. When one of those changed, old rejections were up for retesting, and quite a few of the wins between experiments 105 and 160 came from retesting.

What agents are good and bad at with a profiler

The thing that impressed me most was how well the agents profiled. Without being told how, they ran Nsight Compute, sampled stalls line by line in the SASS, found the hot spots and fixed them, much like a human engineer would. The cluster barrier fix above came straight out of an agent’s own line-by-line stall sampling. They also found things I would never have looked for. This is how the agent from the mixed-precision FMA section announced its find in the notes:

THE INSTRUMENT DISCOVERY: sm_100 has FHFMA — PTX fma.rn.f32.f16 (accepted by ptxas 12.9 AND cu13 for plain sm_100): a SCALAR mixed FMA f32acc += f16*f16 in ONE SASS instruction.

Where they struggled was interpreting what they saw. Profilers report symptoms, and turning a symptom into “these bytes are waiting on that thing” takes first-principles reasoning about the hardware. A few examples from the notes:

• Nsight flagged 82 percent “excessive” uncoalesced global accesses in the reduction kernel. They were stores that nothing waited for, so they were off the critical path and fixing them wouldn’t help.
• Warps that wait by sleeping get sampled, but not on the line that is waiting. One wait took its slow path on 97 to 99 percent of columns and showed zero samples at that line in the profile.
• An agent spent an experiment removing an integer division from the address math. The division ran once per tile on a single thread, off the hot path, and the change made n512 slower, 27.6 to 27.9 ms.
• A very early entry in the log calls the n512 reduction “a measured memory-latency floor; every tuning axis dead” when the whole n512 case took 42 ms. That case finished at 17.8 ms on the grader.

That last kind of mistake was the most expensive. Agents are quick to declare a floor, and a floor that sounds rigorous stops everyone from trying. Around experiment 380 I went back over all the “floors” and found that several had been calibrated against a guessed target rather than a measurement. Reopening them gave several more wins.

Did it overfit?

An autonomous search optimizes exactly what you measure, so the real question is what happens on inputs the benchmark doesn’t contain. Structure-specific routes, like the clustered-case route above, were allowed only if they detect the structure on the actual input and fall back to the general route when the check fails.

Halfway through, an agent built a wider check: it ran the shipped solver on 22 input families, many of them not in the benchmark, over about ten seeds each. It found real overfitting. At n2048, a deflation threshold in the divide and conquer passed both official test seeds, but nearby seeds of the same families failed the checker by up to 2.7×. The agent’s verdict:

The two OFFICIAL n2048 test seeds … BOTH pass at 1e9 — pure luck of benign draws.

Moving the threshold from 1e9 to 1e8 fixed it and cost nothing measurable (40.5 to 40.6 ms). The same check found that the n176 Jacobi kernel’s stopping tolerance had been tuned so that the easy benchmark input passed, and that the n32 kernel overflowed on inputs with very large entries because it skipped the prescale step. Both of those fixes cost speed, 1 to 3 percent, and both Jacobi kernels were replaced by other routes before the end.

The benchmark still missed things. The production port below fails one adversarial “arrowhead” input from the LAPACK test set at n352.

What I would do again

Narrow missions worked much better than broad ones. An agent asked to improve the geomean will happily trade a 5 percent loss on one shape for an 8 percent win on another. An agent asked to make n1024 mixed faster with zero regression elsewhere separates the code paths so its change can’t touch other shapes, and the ratchet enforced that.

Every claimed win had to pass the full check on the merged code, not just in the agent’s worktree. The one regression that slipped through did so because one check lived in the merge script and not in the shared check script.

The most useful view for me was a small dashboard with a bar per case showing the first and newest time. The few times I stepped in, it usually started from a bar that hadn’t moved.

The code

After the competition, I rewrote the shipped kernel to be more readable and stable in a full training run. It’s 3,254 lines instead of the 13,119-line submission file. This was achieved by deleting a bunch of shape specific branches. It measures 7.8 ms on the same 13 cases.

The code is on GitHub.

Thanks

Thanks to Verda for generously providing a B200, and to GPU MODE for organizing the competition.