Is A·B·C cheaper as one three-way product?

Computing A·B·C normally means two matrix multiplications, (A·B)·C. With Strassen, that costs 7 + 7 = 14 scalar multiplications for 2×2 matrices. Can an algorithm that sees all three matrices at once do better? The answer depends on what you count as the cost.

The common thread: every gain comes from deciding what to do with the intermediate product A·B or B·C. See Takeaways. Follow-up: what if matrices were triangles, pentagons or hexagons?

1 · Fused 3-way products: rank exactly N⁴

Pair the output with a fourth matrix D. Then A·B·C is the 4-linear form TN(A,B,C,D) = trace(ABCD) = Σ AijBjkCklDli. A "3-way algorithm" writes T_N as a sum of R terms, each a product of three linear forms (one each in A, B, C) scattered into the outputs. The (A,C | B,D) flattening sends every index tuple (i,j,k,l) to its own column, so it is a permutation-like 0/1 matrix of rank N⁴. Tensor rank (and border rank) is at least any flattening rank, and the naive algorithm uses N⁴ terms. So R(T_N) = N⁴: 16 for 2×2, worse than the Strassen chain's 14.

Two-level schemes can only tie the chain

Suppose level 1 computes s products ℓ(A)·m(B), and level 2 multiplies t combinations of those by forms c(C). The span of the level-1 products must contain the four forms (AB)_ij, so s ≥ R(⟨N,N,N⟩). By symmetry t ≥ R(⟨N,N,N⟩). The chain attains both bounds. At N = 3, beating 46 here would mean 3×3 matrix multiplication in fewer than 23 products, a famous open problem.

2 · The search for 13 multiplications

The next model is broader: any circuit whose gates have degree ≤ 3, where linear forms may mix entries of A, B and C, which allows Winograd-style commutative tricks. Constants and linear terms provably don't help, so the question reduces to (ABC)i = Σr γir vr Σk αrk ℓk mk with cost s + t. Levenberg–Marquardt fitting started from 60 random points per split:

The 6 + 7 near-hits are not algorithms. Under ridge continuation (minimise fit² + μ‖p‖² while μ → 0) the error falls only as the coefficients diverge, with residual roughly ∝ ‖p‖⁻³. That is the signature of an order-1 border decomposition: each product carries three ε^−1/3 factors, the 1/ε terms cancel, and an O(ε) error remains. The same thing happens for A·B alone. 6 commutative products approximate 2×2 A·B arbitrarily well, but an exact algorithm needs 7 (Winograd 1971). So the "13" is just a border-6 A·B followed by the ordinary 7 for ·C, which is not a 3-way effect.

Exact check over GF(2)

The same model is encoded as SAT (CaDiCaL) and matched coefficient by coefficient as formal polynomials over GF(2). Finite fields have no border phenomena, so each answer is exact for GF(2), but says nothing directly about ℝ.

3 · EML trees: x·y·z costs 25, not 33

Odrzywołek showed that eml(a,b) = exp(a) − ln(b) together with the constant 1 generates every elementary function. His direct search found x·y at RPN length 17. An exhaustive search written for this page reproduces his values (x−y 11, −x 15, 1/x 15, x² 17, x·y 17). The x·y tree is a log-domain accumulator: s = e − ln x, then s ← eml(ln s, y) = s − ln y, then exp(e − s). Each additional factor costs 8 tokens, giving 8n + 1 instead of 16(n−1) + 1 for nested trees.

Why this doesn't help the matrix question. The EML tree has no multiplications in it: it is exp(ln x + ln y + ln z), and the saving is just dropping the exp→ln round trip between the two factors. In EML, products are cheap and addition is expensive (x+y needs 19 tokens). Each entry of A·B·C is a sum of four triple products, so the EML cost is dominated by additions.

4 · When writes are expensive and reads are free

This is the setting of non-volatile and asymmetric memory. Reads cost nothing, every word written to memory costs one, and a small scratch area of S registers is not counted as memory. The N² outputs must be written, so W ≥ N². Here three-way fusion does pay off: the intermediate A·B never has to be stored. Every schedule below runs in your browser, instrumented and checked against a direct product.

• Column fusion: hold column l of B·C in N registers, then compute D[:,l] = A·(B·C)[:,l]. That gives W = N² with exactly the chain's 2N³ multiplications.
• With fewer registers: S output accumulators and recomputed B·C entries give N³ + N⁴/S multiplications, still at W = N². This beats the chain whenever ω > N²/S − N.
• e/ln domain: a triple product is one exp of a sum of three logs. The fully fused sum then costs N⁴ exps with a single register. At N = 2 that is 16, a tie with the chain, but with half the writes and no scratch. Logs are recomputed on every read, which is free here; storing them would cost 3N² writes.

5 · Repeated values: caching products

If the matrices contain only k distinct values, the same scalar products recur and can be cached. Quantised neural-network weights (4-bit means k = 16) and ternary weights are examples. A cache keyed by value pair is the same thing as grouping by the distributive law: Dil = Σa a · (Σj : Aij=a (BC)jl) needs one multiplication per distinct value in the row. Four schedules, each instrumented and checked against the exact product:

• chain, no cache: 2N³ multiplications.
• chain + product cache: memoise a·b, then t·c. The intermediate A·B has many distinct values, but each one meets only k values of C, so the second stage costs about N²k. In effect the cache performs the grouping automatically.
• pair cache + grouping, A·(BC): at most k² products for B·C, then one multiplication per distinct value in each row of A. That is ≈ k² + N²k with a table of only k² entries.
• fused triple table: cache each value triple a·b·c once. The output is then pure additions, and the multiplications number at most k² + k³ whatever N is. This is the genuinely three-way effect, but it costs N⁴ additions instead of ≈ N³.

Results at N = 32

What this shows

• Caching a product is the same as grouping equal values. For example, a·x + a·y = a·(x + y) needs one multiplication instead of two.
• The ordinary chain benefits more than expected. Each entry of A·B meets only k different values of C, so a cache brings the chain down to about k² + N²k multiplications. That needs no three-way method.
• The genuinely three-way effect is the fused triple table. Compute every a·b·c once, then build the result by additions alone. That is at most k² + k³ multiplications however large N is (80 at k = 4, against 65,536 for the plain chain at N = 32).
• The price is additions: N⁴ instead of about N³. Compared with grouping, the triple table saves multiplications only when k < N (k³ < N²k). It is worth it only where a multiplication costs far more than an addition.
• With all-distinct values, caching gains nothing, and the triple table is the worst option (N³ + N⁴ multiplications).
• In the write-expensive setting of section 4, the cache itself costs writes unless it fits in scratch. The k² pair table usually fits; the k³ triple table may not.

This idea is well established for two-matrix products: the Four Russians method for Boolean matrices, the Mailman algorithm for small alphabets, and lookup-table kernels for quantised LLMs (LUT-GEMM, T-MAC). The triple-table trade-off for three factors was not found written up anywhere.

6 · An activation after the product: skip the saturated outputs

Neural networks rarely use A·B·C bare; they apply an activation: Y = f(A·B·C + β). Where f is flat, the exact value doesn't matter. For a sign only the side of 0 counts. A ReLU gives 0 for every negative input. A hardtanh or sigmoid is saturated far from 0. So the plan is to compute everything cheaply at b bits, attach a strict error bound that turns each output into an interval, and settle every output where f is flat across its interval. Only the undecided outputs get a deeper, full-precision multiplication. Cost here counts bit operations: a p-bit × q-bit multiply costs p·q, and full precision is 16 bits.

What this shows

• Across the whole chain the cascade costs more than computing exactly (≈ 107–181% at N = 32). Refining one output of A·B·C needs a whole row of A·B or column of B·C, about N² work rather than one dot product. The undecided outputs land in almost every row and column. (The fewest such rows and columns are picked exactly, as a minimum vertex cover via König's theorem.) The strict bound also piles up error through B·C and is about 30× too cautious on this data.
• Compute B·C exactly once and run the cascade on A·(B·C): a refined output is then one dot product, and the bound stops compounding. At N = 32 with one 8-bit pass: sign 52%, hardtanh and saturated sigmoid 60%, ReLU with bias −4 65%, plain ReLU 81% of the exact cost. Every decision is guaranteed correct. The floor is 33%, the exact B·C.
• ReLU saves least: only negative outputs can be skipped, because positive ones need their exact value.
• One well-chosen level beats several: a 4-bit pass decides almost nothing and only adds cost. Precision has to grow with N: 8 bits is best at N = 16–32, 10 bits at N = 64.

7 · Picking terms at random

Instead of lowering precision, estimate each output from a random sample of its terms and decide where it falls relative to a threshold (a sign or ReLU after the product). For three factors there's an appealing twist. Each output is a double sum, D_il = Σ_jΣ_k A_ijB_jkC_kl, so you can sample (j, k) pairs directly at 2 multiplications each, without ever building a row or column of the intermediate. Cost here is multiplications; the exact chain is 2N³.

What this shows

• Sampling terms mostly loses. A strict bound (Hoeffding–Serfling) never settles a single output: it has to allow for the largest possible term, so it is always too wide. An empirical bound saves only 1–17% on dense data, with a few wrong decisions. On sparse data it looks cheap but is badly fooled, with hundreds of wrong decisions: the first samples are often all zeros, so the estimate looks falsely certain.
• Why: to decide an output, the error must be smaller than its distance to the threshold. When outputs crowd around the threshold, partial samples can't get there. Seeing every term roughly (low precision) beats seeing some terms exactly (sampling).
• The exception is a random walk through all three matrices. For nonnegative data, walk i → j → k → l with j ∝ A_ij·(B·C·1)_j, k ∝ B_jk·(C·1)_k, l ∝ C_kl. The walk lands on output l with probability exactly D_il / Σ_l D_il, after a setup of only 5N² multiplications and without forming A·B or B·C. The visit counts give a guaranteed bound. On sparse, heavy-tailed data with a top-1% threshold, that needs only about a quarter of the multiplications, with zero wrong decisions. This is Cohen–Lewis sampling, extended to a chain of three.
• The catch: walks trade multiplications for random draws (about 393,000 draws at 4,096 walks per row, against 65,536 multiplications for the exact chain). They pay only when a multiplication is much more expensive than drawing a random index, the data is nonnegative and heavy-tailed, and the question is which few outputs are large.

Takeaways

The common thread is the intermediate product. Every gain came from treating A·B (or B·C) deliberately rather than as a black box:

• Writes: don't store it at all. Fusing halves the writes (section 4).
• Caching: cache products of the original entries. The intermediate's entries are sums, so they repeat far less than the inputs (section 5).
• Activations: compute the middle product exactly and save the cheap tricks for the last multiplication. Treating the chain as one block made things worse (section 6).
• Random walks: the one random method that works passes through the whole chain and never forms the intermediate at all. It weights each step by the rest of the chain instead (section 7).

Surprises

1. Applying the shortcut to the whole chain cost more than doing nothing (107–181%). In a single product, fixing one uncertain output costs one dot product (N work). In a chain it costs a whole row of A·B or column of B·C (N² work), so the smallest unit of correction is a row or column, not an entry.
2. The guaranteed bound was about 30× too cautious. With the bound shrunk to ≈ 0.03×, the whole-chain cascade reaches 50% with no wrong decisions on this data, beating the guaranteed hybrid (52%). Which strategy wins depends on whether you need a guarantee.
3. More precision levels made it worse. A 4-bit first pass settled almost nothing and only added cost; one well-chosen level was best.
4. ReLU, the most common activation, benefits least, because only its negative outputs can be skipped.
5. A classic graph result turned up naturally: choosing the fewest rows and columns to recompute is a minimum vertex cover in a bipartite graph, solved exactly by König's theorem.
6. Caching helps the plain chain more than expected. Each entry of A·B meets only k values of C, so a cache brings even the textbook chain to about k² + N²k multiplications.
7. Random sampling mostly loses to low precision, and its empirical bounds are dangerously overconfident on sparse data. To decide a threshold, seeing every term roughly beats seeing some terms exactly.

For practice

• Order the work so the product feeding the nonlinearity comes last. Keep the middle product exact and in fast memory, and apply early exit, sparsity or caching only to that last product.
• Fusing three-factor products pays off in memory traffic, not in arithmetic. That's why GPU libraries fuse back-to-back GEMMs and FlashAttention fuses attention.
• With repeated values, cache before the sums. A full table of value triples only pays off when k < N and multiplications cost far more than additions.
• Measure how loose your bounds are on real data before choosing between guaranteed and predictive skipping. Here the gap was 30×.
• Reach for random sampling only to find the few largest outputs of nonnegative, heavy-tailed products, and then with walks through the chain, not with per-term sampling.

Possibly new, and what to do next

Most results here rediscover known theory or practice (see the table below). Four may be new:

• six commutative products approximating 2×2 A·B arbitrarily well (numerical only, no explicit formula yet);
• for three factors, the activation shortcut belongs on the last product, not the whole chain;
• the fused triple table: at most k² + k³ multiplications for any N, at the price of N⁴ additions;
• guaranteed threshold decisions from random walks through all three factors, weighted by the rest of the chain.

Next steps:

• derive the explicit 6-product formula, which settles whether it is real;
• replace the worst-case bound with a probabilistic one that has a stated failure rate, to close most of the 30× gap;
• rerun the caching and activation experiments on trained network weights instead of random matrices;
• settle the remaining GF(2) splits (6+7, 7+6, 5+8): plain SAT gave no answer in 11 hours, so they need symmetry breaking or a smarter encoding.

Prior work, and how these results compare

Three-way matrix products have been studied from several directions. Most of what this page finds either restates known theory for the chain tensor or rediscovers established practice. The table puts each result next to the closest published one.

What is proven, what is numerical