BLAS Level 4

The BLAS hierarchy has been climbing an arithmetic-intensity ladder — flops performed per byte of memory traffic — since 1979. Level 1 (Lawson–Hanson–Kincaid–Krogh, 1979) standardized vector operations like axpy and dot: O(n) work on O(n) data, intensity O(1). Level 2 (1988) added matrix–vector operations like gemv: O(n²) on O(n²), still intensity O(1) — enough work per call to keep a Cray's vector pipelines full. Level 3 (1990) added matrix–matrix operations like gemm: O(n³) work on O(n²) data, the first level with intensity O(n) — the level where cache blocking lets each byte be reused ~n times. That single property is what made LAPACK and the cache-based machines of the late 1980s fit each other.

The pattern is explicit in the literature: each level raised flops-per-byte to track hardware whose compute outgrew its memory. Level 4 continues it. We define it as the operations that perform Θ(n⁴) arithmetic while touching only Θ(n²) data — intensity Θ(n²), quadratically above gemm.

Between 2016 and 2024, dense half-precision matrix throughput rose about 106× (P100's 21 TFLOPS to B200's 2,250 TFLOPS) while HBM bandwidth rose only about 10.5× (0.73 → 7.7 TB/s). The machine balance point — the flop/byte ratio above which a kernel is compute-bound — marched from ~29 to ~295, and to ~584 (FP8) and ~1,169 (FP4) at the precisions Blackwell is built to sell. Any kernel that cannot supply hundreds of flops per byte is bandwidth-bound silicon.

The canonical skeleton is

with an associative — and normally commutative — reduction ⊕, so GPU tiles can stream sources through shared memory while a bounded output-accumulator state stays register-resident. A stored 4-tensor T_ijkl fails on the spot: it is n⁴ data, and contracting it runs at memory speed. The whole value of Level 4 is that the n⁴ never exists in memory at all.

Scrub the problem size and watch the two regimes diverge. The materialized tensor grows as n⁴; the implicit generator's footprint stays at Θ(n²).

Two associativity properties do all the work. Commutativity + associativity of ⊕ license partitioning the source axis into arbitrary blocks, tree-combining partial states, and splitting the reduction across thousands of GPU cores — none of which is legal for an operation that must see all sources at once. This is exactly the algebra a semiring or a monoid provides, and it is why the accumulator can be a small, register-resident state rather than a full row of the interaction matrix.

The accumulator is a commutative monoid. For softmax it is the triple (m, ℓ, o) = (running max score, rescaled denominator, rescaled weighted value), merged between any two blocks by

m = max(m₁, m₂)
ℓ = e^(m₁−m)·ℓ₁ + e^(m₂−m)·ℓ₂
o = e^(m₁−m)·o₁ + e^(m₂−m)·o₂

Swap that monoid for sum and you get a kernel matrix–vector product; for (min, argmin) you get k-nearest-neighbors; for logsumexp you get a Sinkhorn optimal-transport step or a kernel-density estimate. Same dataflow, same register discipline, different three-line merge. That interchangeability is the thesis of the whole report — and the reason a single primitive can stand behind a dozen named operations.

// B4MR — the Level-4 map-reduce. Sources never leave the tile pipeline.
b4mr(
    targets  X[M, d],   sources Z[N, d],   values V[N, e],   theta,
    F:     (x_p, z_q, v_q; theta) -> partial,   // implicit interaction, c_F flops
    OPLUS: (state, partial) -> state,           // associative, bounded O(1) state
    PHI:   (state) -> out,                       // per-target finalize
    rank:  { VECTOR | SCALAR | OP1 }              // output shape / epilogue
) -> Y

// rank = VECTOR : Y[M,m]  — PAIRREDUCE / implicit-operator apply (attention, kernel sums, N-body)
// rank = SCALAR : Y[b]    — QFORM: fuse a final DOT / LSE / histogram, zero extra traffic
// rank = OP1    : Y[n,n]  — OUTERACC: each pair contributes a rank-1 update u_q u_qᵀ

// Monoids ship as a tested enum:
//   SUM · LOGSUMEXP (online m,ℓ state) · MIN/MAX (+ argmin/argmax) · PROD
// F ships from a fixed generator catalog (Gaussian, Coulomb, inner-product-then-monoid);
// a JIT tier exposes user F for the open family — labeled library-level, not the ABI.

Verification forced one refinement to the original five-candidate proposal. The all-pairs core is genuinely a family — PAIRREDUCE is the general primitive, implicit-kernel apply is its coordinate/law-generated form, and the scalar QFORM is its terminal epilogue, not a separate operation. What robustly earns its own entry point beyond the core is the non-separable trilinear lift (TRIATTN) and the already-shipping scientific instance (the SYM2E Fock build). The shift- and group-indexed members and the time-indexed recurrence are structured-index or serial-composition variants whose high-value instances fall to FFT, doubling, or low-rank shortcuts.

Over an ordinary ring, the FFT demolishes this in O(n² log n) — direct n⁴ is indefensible past moderate n. But the (min,+) and (max,+) tropical semirings have no additive inverses, so no transform applies; under the standard MinConv conjecture the arbitrary-cost tropical case admits nothing substantially below Θ(n⁴). It is the only candidate in the roster carrying a hardness certificate — and the reason a semiring flag, not just a faster convolution, is the right abstraction. The same move GraphBLAS made for sparse matrix multiply.

b4tri earns a separate entry point because its reduction index is an ordered pair of sources and, under softmax, the trilinear score provably does not factor — so the direct n³·d form is forced for the exact variant. b4conv is standardized only for semirings without additive inverses, where FFT and NTT do not apply; for the ordinary ring case the library must dispatch to cuFFT/cuDNN. Deliberately excluded: the fused recurrence (a serial fold with Θ(n) critical-path depth, dominated by repeated squaring) and span-DP (a recurrence over its own outputs) — both ship as library-level recipes, not primitives.

Try the ranking yourself. The four axes are the discriminators the skeptics scored; drag the weights to your own priorities and watch the order re-settle. The default weighting reproduces the panel's ranking exactly.

The gates read: (1) implicit generation — are the n⁴ interactions generated from Θ(n²) data, never stored? (2) one-shot fold — is it a single associative reduction over stored inputs, not a serial composition or a recurrence over its own outputs? (3) shortcut-proof at scale — does real workload survive FFT / FMM / low-rank / doubling, or do shortcuts dominate almost everywhere? (4) bounded state — does a fixed-size accumulator stream on the GPU without atomics or serial depth? The rejections each fail a different gate, which is precisely their instructive value: Level 4 is a one-shot associative fold over stored inputs — not a serial composition, not a shortcut-dominated structured index, not a recurrence over its own outputs.

The orthogonal effort is tensor-contraction standardization — cuTENSOR, TAPP, GETT, TBLIS, the Tensor Contraction Engine. But a general contraction of stored order-4 tensors materializes O(n⁴) data, which is precisely what the disqualification rule forbids; those operations remain BLAS-3-like in flops-per-byte-of-materialized-data. Level 4 is the genuinely higher rung: the n⁴ exists only in flight.