fix(ops): build multiply_2d in O(n) via analytic quarter-square; guard SVD

multiply_2d was infeasible to BUILD at high breakpoint counts. It lowered to
piecewise_linear_2d's least-squares solve, whose `numpy.linalg.svd(...,
full_matrices=True)` materializes a discarded (n1*n2, n1*n2) left-singular
matrix U. At 257 breakpoints/axis that U is ~35 GB of float64 — a single
graph build OOM'd a 30 GB machine.

Two layered fixes:

1. Analytic product fast path (the named-op fix). A product x*y is exactly
   bilinear, so multiply_2d no longer routes through the generic solve. On
   the common unit grid it builds the piecewise-linear interpolant in closed
   form via the quarter-square identity x*y = ((x+y)^2 - (x-y)^2)/4. The
   interpolant of t^2 on a uniform grid of spacing h has a constant slope
   change of 2h at every interior breakpoint, so each square is an
   equal-weight ReLU staircase: one MLP sublayer, ~2*(2n-2) neurons, no
   solve. The generic piecewise_linear_2d path is untouched for true 2D
   functions. Build at n=257 over [-1500,1500]: OOM -> 230 ms / 0.44 MB.
   Noise is unchanged (multiply_uniform_pm10 0.06250 -> 0.06249, within GPU
   FP jitter); the interpolant is the same diagonal-aligned product.

2. Conditional full_matrices guard in piecewise_linear_2d. The fit reads the
   nullspace basis Vh[rank:].T, which needs every right-singular vector up to
   the column dimension. When A_v is TALL (uniform/collapsed grids) the
   reduced SVD already returns all of them, so full_matrices=False is
   bit-identical and skips the giant U. When A_v is WIDE (non-uniform grids,
   where the sum/diff hyperplane families don't collapse) the reduced SVD
   drops the nullspace rows, so full_matrices=True is required there — an
   earlier unconditional flip was the source of the historical 19.4 vs 7.78
   diff_trig regression. The guard keeps the non-uniform path at 7.78 while
   making the tall path buildable.

Tests: add a build-cost regression and a multi-chunk regression to
test_multiply_2d.py; add scripts/validate_multiply_2d_build_cost.py. D7:
re-measured noise, refreshed footers/JSON/markdown, and updated
numerical_noise_findings.md (the stale 19.4 piecewise_linear_2d finding is
resolved; added a quarter-square note).

Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>

Author: physicsrob

docs/numerical_noise.md

@@ -1,7 +1,7 @@
 # Numerical noise reference
 
 *Generated by `scripts/measure_op_noise.py`; do not edit by hand.*
-*Measured at commit `c234ba4` on `2026-04-29`.*
+*Measured at commit `fe133f9` on `2026-06-02`.*
 
 Every piecewise-linear op in `torchwright/ops/` is measured on one or more
 named input distributions characterising its expected input ranges.
@@ -26,13 +26,13 @@ For commentary on the numbers — which are expected, which warrant investigatio
 | `exp` | `torchwright.ops.arithmetic_ops` | 0.02797 | 0.0001924 | `exp_pm5` |
 | `floor_int` | `torchwright.ops.arithmetic_ops` | 0.9993 | 0.953 | `floor_near_boundary_10` |
 | `linear_bin_index` | `torchwright.ops.arithmetic_ops` | 1 | 0.9854 | `linear_bin_index_tex_col_16` |
-| `log` | `torchwright.ops.arithmetic_ops` | 0.003026 | 0.003721 | `log_6decades_wide` |
+| `log` | `torchwright.ops.arithmetic_ops` | 0.003023 | 0.003705 | `log_6decades_wide` |
 | `log_abs` | `torchwright.ops.arithmetic_ops` | 0.000699 | 0.1428 | `log_abs_3decades_pm100` |
-| `low_rank_2d` | `torchwright.ops.arithmetic_ops` | 0.0651 | 0.3314 | `atan_cross_dot_nonuniform` |
+| `low_rank_2d` | `torchwright.ops.arithmetic_ops` | 0.0651 | 0.3315 | `atan_cross_dot_nonuniform` |
 | `max` | `torchwright.ops.arithmetic_ops` | 3.815e-06 | 1.953e-06 | `minmax_uniform_pm50` |
 | `min` | `torchwright.ops.arithmetic_ops` | 3.815e-06 | 1.953e-06 | `minmax_uniform_pm50` |
 | `mod_const` | `torchwright.ops.arithmetic_ops` | 0 | 0 | `mod_integers_0_100_by7` |
-| `multiply_2d` | `torchwright.ops.arithmetic_ops` | 0.0625 | 2.15 | `multiply_uniform_pm10` |
+| `multiply_2d` | `torchwright.ops.arithmetic_ops` | 0.06249 | 2.15 | `multiply_uniform_pm10` |
 | `multiply_integers` | `torchwright.ops.arithmetic_ops` | 0 | 0 | `multiply_integers_0_10` |
 | `piecewise_linear` | `torchwright.ops.arithmetic_ops` | 0.25 | 231.1 | `parabola_0_10_step1` |
 | `piecewise_linear_2d` | `torchwright.ops.arithmetic_ops` | 7.778 | 0.6024 | `diff_trig_nonuniform` |
@@ -333,12 +333,12 @@ log(x) over [0.01, 100] — 4 decades. With per-decade sectioning (default), eac
 | Metric | Value |
 | --- | --- |
 | Samples | 4096 |
-| Max abs error | 0.000247 |
-| Mean abs error | 0.0001126 |
-| p99 abs error | 0.0001694 |
-| Max rel error | 0.003721 (over 4096 samples with |ref| ≥ 1e-6) |
-| Mean rel error | 4.257e-05 |
-| p99 rel error | 0.0002646 |
+| Max abs error | 0.0002475 |
+| Mean abs error | 0.0001119 |
+| p99 abs error | 0.0001681 |
+| Max rel error | 0.003705 (over 4096 samples with |ref| ≥ 1e-6) |
+| Mean rel error | 4.237e-05 |
+| p99 rel error | 0.0002638 |
 | Worst input | `x=0.01797` |
 
 ### `log_6decades_wide`
@@ -348,12 +348,12 @@ log(x) over [0.01, 30000] — 6 decades. Stresses the sectioning path: a single-
 | Metric | Value |
 | --- | --- |
 | Samples | 4096 |
-| Max abs error | 0.003026 |
+| Max abs error | 0.003023 |
 | Mean abs error | 0.001523 |
-| p99 abs error | 0.002798 |
-| Max rel error | 0.001149 (over 4096 samples with |ref| ≥ 1e-6) |
-| Mean rel error | 0.0001697 |
-| p99 rel error | 0.0003993 |
+| p99 abs error | 0.002797 |
+| Max rel error | 0.001148 (over 4096 samples with |ref| ≥ 1e-6) |
+| Mean rel error | 0.0001696 |
+| p99 rel error | 0.0003989 |
 | Worst input | `x=2.707e+04` |
 
 ## `log_abs`
@@ -389,10 +389,10 @@ atan2(cross/dot) over [-2, 2] × [0.5, 10] on a non-uniform grid — rank-3 SVD
 | Max abs error | 0.0651 |
 | Mean abs error | 0.01586 |
 | p99 abs error | 0.05744 |
-| Max rel error | 0.3314 (over 4096 samples with |ref| ≥ 1e-6) |
+| Max rel error | 0.3315 (over 4096 samples with |ref| ≥ 1e-6) |
 | Mean rel error | 0.08032 |
 | p99 rel error | 0.2172 |
-| Worst input | `a=2`, `b=3.365` |
+| Worst input | `a=-2`, `b=3.365` |
 
 ## `max`
 
@@ -453,7 +453,7 @@ Integer inputs on [0, 100] with divisor=7. `mod_const` is x − divisor × `ther
 
 ## `multiply_2d`
 
-Single-sublayer product via 2D piecewise-linear lookup. Analytical worst-cell absolute bound is `step1*step2/4 = 0.25` for this configuration.
+Single-sublayer product built analytically via the quarter-square identity `a*b = ((a+b)^2 - (a-b)^2)/4` (O(n) neurons, no least-squares solve). Worst-cell absolute bound is `step1*step2/4 = 0.25` for this configuration.
 
 ### `multiply_uniform_pm10`
 
@@ -462,7 +462,7 @@ a·b over [-10, 10]² with step=1 — generic uniform product grid; analytical b
 | Metric | Value |
 | --- | --- |
 | Samples | 4096 |
-| Max abs error | 0.0625 |
+| Max abs error | 0.06249 |
 | Mean abs error | 0.0202 |
 | p99 abs error | 0.05847 |
 | Max rel error | 2.15 (over 4096 samples with |ref| ≥ 1e-6) |

docs/numerical_noise_findings.md

@@ -15,28 +15,60 @@ observations and removing findings a fix has invalidated. See the
 
 ### `piecewise_linear_2d` oscillates on non-uniform grids
 
-`diff_trig_nonuniform` hits **19.4 absolute error** on a product whose
-reference value is `-24.3` — a ~80% relative deviation.
-`diff_vel_nonuniform` hits **13.6 absolute error**. The op's own docstring
+`diff_trig_nonuniform` hits **7.78 absolute error** on a product whose
+reference value is `-24.3` — a ~32% relative deviation.
+`diff_vel_nonuniform` hits **5.44 absolute error**. The op's own docstring
 already flags this: the constrained least-squares fit used on non-uniform
 breakpoint grids can oscillate inside cell interiors. **Preferred
 replacement:** `low_rank_2d`, which measures 0.065 absolute / 33% relative
 (rank-3 SVD, with a compile-time `σ_{K+1}` bound).
 
-**Noise ~doubled after commit `3895411`** ("perf(ops): drop full_matrices in
-piecewise_linear_2d SVD"). That change moved `torch.linalg.svd` from
-`full_matrices=True` to `False` to save ~15 GB of unused U. For `m > n` the
-two SVD modes compute the same `Vh` in theory, but torch's reduced-SVD
-algorithm produces a slightly different orthonormal basis in the nullspace
-span (which is what the piecewise_linear_2d fit uses via `Vh[rank:]`).
-Different nullspace basis → different constrained least-squares solution →
-different oscillation pattern, same shape but higher amplitude at DIFF_BP's
-sparsely-sampled outer cells. The commit did not re-measure noise (D7
-violation caught by `test_numerical_noise_drift` and resolved on this
-branch). The underlying issue — piecewise_linear_2d is brittle on highly
+**The SVD's `full_matrices` flag is load-bearing here, and only for *tall*
+design matrices.** `piecewise_linear_2d` reads the nullspace basis
+`N_basis = Vh[rank:].T`, which needs every right-singular vector up to the
+column dimension `3+K`. The shape of `A_v` decides whether the two SVD
+modes agree:
+
+* **Tall** (`vertex rows n1·n2 >= 3+K`): this happens on uniform / collapsed
+  grids — e.g. the unit grid `multiply_2d` used to feed here, where the
+  sum/diff hyperplane families collapse to O(n) distinct lines. The reduced
+  SVD (`full_matrices=False`) already returns all `3+K` right-singular
+  vectors, so `Vh` — and therefore the solution — is identical to the full
+  SVD, while skipping the discarded `(n1·n2, n1·n2)` left-singular matrix U
+  (~35 GB of float64 at n1·n2 ≈ 66049 — a real build-time OOM).
+* **Wide** (`rows < 3+K`): this is the *non-uniform* case — DIFF_BP × TRIG_BP
+  has 117 vertex rows but 219 hyperplane columns. The reduced SVD returns
+  only the first 117 right-singular vectors, so the nullspace rows
+  `rank..(3+K)` are *missing*. `N_basis` then can't span the nullspace the
+  interior-sample regularization needs, the oscillation is no longer pinned,
+  and `diff_trig_nonuniform`'s max abs error jumps from 7.78 to **19.4**.
+
+So the op now guards the flag: `full_matrices=False` only when `A_v` is
+tall, `full_matrices=True` (the U it builds is only `(rows, rows)`, and the
+non-uniform grids that hit this path are small) otherwise. An earlier commit
+that flipped the flag unconditionally was the source of the historical 19.4
+regression; the conditional guard keeps the non-uniform path at its
+original 7.78 while making the tall path buildable at high breakpoint
+counts. The underlying issue — `piecewise_linear_2d` is brittle on highly
 non-uniform grids — is unchanged; callers should prefer `low_rank_2d` or a
 uniform grid.
 
+### `multiply_2d` builds its product grid analytically (quarter-square)
+
+`multiply_2d` no longer routes the product through `piecewise_linear_2d`'s
+generic least-squares solve. A product is exactly bilinear, so on the
+common unit grid the op now builds the piecewise-linear interpolant in
+closed form via the quarter-square identity
+`u·v = ((u+v)² − (u−v)²)/4`: the interpolant of `t²` on a uniform grid of
+spacing `h` has a *constant* slope change of `2h` at every interior
+breakpoint, so each square is an equal-weight ReLU staircase and the whole
+product is one ReLU bank of `~2·(2n−2)` neurons — O(n) build, no SVD, no
+O(n²) design matrix (the old path OOM'd at ~257 breakpoints/axis). This is
+the same single MLP sublayer and the same diagonal-aligned interpolant, so
+`multiply_uniform_pm10` is unchanged at ~0.0625 absolute (it moved by ~1e-5,
+within GPU FP run-to-run jitter). The near-zero high-relative-error
+pathology below is structural to the product and unaffected.
+
 ### `reciprocal` relative error is load-bearing at the high-x tail
 
 `reciprocal_03_200` reports 0.007 absolute error, which sounds small. But

docs/op_noise_data.json

@@ -1,7 +1,7 @@
 {
   "schema_version": 1,
-  "commit": "c234ba4",
-  "measured_at": "2026-04-29",
+  "commit": "fe133f9",
+  "measured_at": "2026-06-02",
   "ops": [
     {
       "name": "abs",
@@ -315,12 +315,12 @@
           "name": "log_4decades_001_100",
           "description": "log(x) over [0.01, 100] — 4 decades. With per-decade sectioning (default), each section's pre-cancellation magnitude is bounded by `section_factor=10`, so float32 ULP noise is ~1.2e-6 per section. The dominant residual error comes from the multiply_2d blending grid in routing.",
           "n_samples": 4096,
-          "max_abs_error": 0.000247002,
-          "mean_abs_error": 0.000112631,
-          "p99_abs_error": 0.000169402,
-          "max_rel_error": 0.00372105,
-          "mean_rel_error": 4.25736e-05,
-          "p99_rel_error": 0.000264624,
+          "max_abs_error": 0.000247478,
+          "mean_abs_error": 0.000111867,
+          "p99_abs_error": 0.000168088,
+          "max_rel_error": 0.00370548,
+          "mean_rel_error": 4.2368e-05,
+          "p99_rel_error": 0.000263829,
           "rel_valid_samples": 4096,
           "worst_input": {
             "x": 0.0179683
@@ -330,12 +330,12 @@
           "name": "log_6decades_wide",
           "description": "log(x) over [0.01, 30000] — 6 decades. Stresses the sectioning path: a single-section piecewise log over this range fails outright (cancellation floor `(x_max/x_min) · 2⁻²³ ≈ 0.36 absolute, observed worst ~1.0 abs). Sectioning drops the floor to ~5e-3 absolute, dominated by compare-cancellation noise at large `x` propagating through multiply_2d blending.",
           "n_samples": 4096,
-          "max_abs_error": 0.00302601,
-          "mean_abs_error": 0.00152333,
-          "p99_abs_error": 0.00279808,
-          "max_rel_error": 0.00114918,
-          "mean_rel_error": 0.000169692,
-          "p99_rel_error": 0.000399339,
+          "max_abs_error": 0.00302315,
+          "mean_abs_error": 0.00152277,
+          "p99_abs_error": 0.00279686,
+          "max_rel_error": 0.00114823,
+          "mean_rel_error": 0.000169627,
+          "p99_rel_error": 0.000398892,
           "rel_valid_samples": 4096,
           "worst_input": {
             "x": 27070.5
@@ -374,15 +374,15 @@
           "name": "atan_cross_dot_nonuniform",
           "description": "atan2(cross/dot) over [-2, 2] × [0.5, 10] on a non-uniform grid — rank-3 SVD via `low_rank_2d`.",
           "n_samples": 4096,
-          "max_abs_error": 0.0651017,
+          "max_abs_error": 0.0651015,
           "mean_abs_error": 0.0158553,
           "p99_abs_error": 0.0574396,
-          "max_rel_error": 0.33142,
-          "mean_rel_error": 0.0803217,
-          "p99_rel_error": 0.217151,
+          "max_rel_error": 0.331451,
+          "mean_rel_error": 0.0803216,
+          "p99_rel_error": 0.21716,
           "rel_valid_samples": 4096,
           "worst_input": {
-            "a": 2.0,
+            "a": -2.0,
             "b": 3.36508
           }
         }
@@ -459,17 +459,17 @@
     {
       "name": "multiply_2d",
       "module": "torchwright.ops.arithmetic_ops",
-      "notes": "Single-sublayer product via 2D piecewise-linear lookup. Analytical worst-cell absolute bound is `step1*step2/4 = 0.25` for this configuration.",
+      "notes": "Single-sublayer product built analytically via the quarter-square identity `a*b = ((a+b)^2 - (a-b)^2)/4` (O(n) neurons, no least-squares solve). Worst-cell absolute bound is `step1*step2/4 = 0.25` for this configuration.",
       "distributions": [
         {
           "name": "multiply_uniform_pm10",
           "description": "a·b over [-10, 10]² with step=1 — generic uniform product grid; analytical bound `step1*step2/4 = 0.25`.",
           "n_samples": 4096,
-          "max_abs_error": 0.0624998,
-          "mean_abs_error": 0.0201988,
-          "p99_abs_error": 0.0584654,
-          "max_rel_error": 2.15046,
-          "mean_rel_error": 0.00729335,
+          "max_abs_error": 0.0624883,
+          "mean_abs_error": 0.0201983,
+          "p99_abs_error": 0.0584726,
+          "max_rel_error": 2.14985,
+          "mean_rel_error": 0.00729313,
           "p99_rel_error": 0.0661755,
           "rel_valid_samples": 4096,
           "worst_input": {

scripts/measure_op_noise.py

@@ -812,9 +812,10 @@ def _target_ops() -> List[TargetOp]:
             reference_fn=lambda inputs: inputs["a"] * inputs["b"],
             distribution_names=("multiply_uniform_pm10",),
             notes=(
-                "Single-sublayer product via 2D piecewise-linear lookup. "
-                "Analytical worst-cell absolute bound is `step1*step2/4 = 0.25` "
-                "for this configuration."
+                "Single-sublayer product built analytically via the "
+                "quarter-square identity `a*b = ((a+b)^2 - (a-b)^2)/4` "
+                "(O(n) neurons, no least-squares solve). Worst-cell absolute "
+                "bound is `step1*step2/4 = 0.25` for this configuration."
             ),
         ),
         TargetOp(

scripts/validate_multiply_2d_build_cost.py

@@ -0,0 +1,93 @@
+"""Validate multiply_2d's build cost and accuracy at high breakpoint counts.
+
+Background.  ``multiply_2d`` used to lower to ``piecewise_linear_2d``'s
+generic least-squares solve, which built an O(n²)-row vertex design matrix
+and ran ``numpy.linalg.svd(..., full_matrices=True)`` whose left-singular
+matrix is ``(n², n²)``.  At ~257 breakpoints per axis that discarded U
+matrix is ~35 GB of float64 — a single graph build OOM'd a 30 GB machine.
+
+A product ``x·y`` is exactly bilinear (rank-1), so ``multiply_2d`` now
+builds its piecewise-linear interpolant analytically via the quarter-square
+identity ``x·y = ((x+y)² − (x−y)²)/4`` — O(n) neurons, no global solve.
+
+This script:
+
+* (i)  builds a 257×257 grid over a realistic operand range and confirms it
+       constructs in well under a second using a trivial amount of memory;
+* (ii) checks the max abs error of the product against exact ``x·y`` over a
+       dense sample, at magnitudes up to ~1500·1500, stays inside the
+       bilinear ``step²/4`` interpolation bound.
+
+Run with the workspace venv::
+
+    /data/torchdoom/.venv/bin/python -m scripts.validate_multiply_2d_build_cost
+"""
+
+from __future__ import annotations
+
+import time
+import tracemalloc
+
+import torch
+
+from torchwright.ops.arithmetic_ops import multiply_2d
+from torchwright.ops.inout_nodes import create_input
+
+
+def main() -> None:
+    max_abs = 1500.0
+    n_bp = 257
+    step = (2 * max_abs) / (n_bp - 1)  # 257 breakpoints per axis over [-1500, 1500]
+
+    # (i) Construction cost.
+    a = create_input("a", 1)
+    b = create_input("b", 1)
+    tracemalloc.start()
+    t0 = time.perf_counter()
+    node = multiply_2d(
+        a, b, max_abs1=max_abs, max_abs2=max_abs, step1=step, step2=step
+    )
+    elapsed = time.perf_counter() - t0
+    _cur, peak = tracemalloc.get_traced_memory()
+    tracemalloc.stop()
+
+    print(f"breakpoints per axis : {n_bp}")
+    print(f"operand range        : [-{max_abs}, {max_abs}], step {step:.4f}")
+    print(f"build time           : {elapsed * 1000:.1f} ms")
+    print(f"build alloc (traced) : {peak / 1e6:.2f} MB")
+    print()
+    print(
+        "old generic path at this n built a (n², n²) = "
+        f"({n_bp ** 2}, {n_bp ** 2}) float64 SVD U matrix "
+        f"≈ {n_bp ** 4 * 8 / 1e9:.1f} GB → OOM"
+    )
+    print()
+
+    assert elapsed < 1.0, f"build took {elapsed:.2f}s — expected well under a second"
+    assert peak < 200e6, f"build allocated {peak / 1e6:.1f} MB — expected a few MB"
+
+    # (ii) Accuracy vs exact product over a dense sample.
+    torch.manual_seed(0)
+    n_samples = 200_000
+    xs = (torch.rand(n_samples) * 2 - 1) * max_abs
+    ys = (torch.rand(n_samples) * 2 - 1) * max_abs
+    worst = 0.0
+    worst_at = None
+    # compute() takes one position at a time here; batch through n_pos.
+    inputs = {"a": xs.unsqueeze(1), "b": ys.unsqueeze(1)}
+    out = node.compute(n_pos=n_samples, input_values=inputs).flatten()
+    err = (out - xs * ys).abs()
+    worst = err.max().item()
+    worst_at = (xs[err.argmax()].item(), ys[err.argmax()].item())
+
+    bound = step * step / 4
+    print(f"max abs error        : {worst:.4f} at {worst_at}")
+    print(f"bilinear step²/4 bnd : {bound:.4f}")
+    assert worst < bound + 1.0, f"error {worst} exceeds bilinear bound {bound}"
+
+    print()
+    print("PASS: builds in O(n) time/memory and stays within the interpolation bound.")
+
+
+if __name__ == "__main__":
+    main()