S1-26: OUT_PROD_ID CPU kernel — the activation half of MUL_MAT_ID's backward

[dillon-blake]

Stacked on #20 (S1-27). With this, MUL_MAT_ID has a complete backward.

d_b[k,l,t] = sum over {i : i % ne_b1 == l} of  sum_j  as[k,j, ids[i,t]] * grad[j,i,t]

as may be quantized, and that is the whole difficulty rather than a corner case. In a LoRA MoE
graph the base expert stacks are frozen Q4_K. They take no gradient themselves, but the activations
flowing into them do — because an earlier layer's LoRA is upstream. So d(b) must propagate back
through a quantized weight: the same dequantize-a-row-at-a-time path as OUT_PROD, accumulating in
F32 per ADR-0002.

MODE_GRAD cannot check the quantized path, and it took a while to work out why

ggml's quantized matmul quantizes the activations on the fly to vec_dot_type before the dot
product. So the forward is a staircase in b: its finite difference at the FD step scale measures
quantization edges, not a derivative. The analytic gradient is the gradient of the intended smooth
function; the harness evaluates the actual quantized one. They disagree by construction — measured
MAA 0.028–0.071, with the kernel exact throughout.

Upstream reached the same conclusion and encoded it without a word of comment:

// test_mul_mat
if (!ggml_is_quantized(type_a)) {
    ggml_set_param(b);        // b is a param ONLY when the weight is not quantized
}

So the quantized cases are deliberately absent from MODE_GRAD, and the dequantize path is verified
a stronger way instead: OUT_PROD_ID with a quantized as, against OUT_PROD_ID with that same as
dequantized to F32. Bit-exact on q8_0, q4_K and q4_0.

Two test fixes, both from measuring rather than assuming

The token count is 32, not 5. mean_abs_asymm divides by (gn + ga) — not (|gn| + |ga|) — so
any output element whose true gradient is near zero sends that ratio to infinity. This op manufactures
such elements: d_as[k,j,e] sums only over the slots that routed to expert e, so with 5 tokens and 3
experts it is a sum of one or two random products. A single element with asymm ≈ 700 drags the mean
of 90 up to ≈ 8.

tokens	result
5	3/10 runs fail, MAA up to 7.9 — with the kernel exact throughout
32	0/30 runs fail

Condition the test; don't widen the tolerance until the flapping stops.

max_maa_err = 5e-3, measured on both sides:

	MAA
worst FD noise, 40 runs	4.5e-4
ignore grad entirely	9.21
ignore the ne_b1 broadcast	3.39 / 1.15
read grad row 0 always	4.41
route every slot to every expert	5.30

10× above the noise floor, 200–1800× below every real defect. Five mutations introduced, five
caught.

Verification

Both kernels checked against a naive double-precision reference across 15 configurations covering
both broadcast modes: relative error ~1e-7, i.e. float32 epsilon. That reference is what proved the FD
failures were the harness and not the arithmetic.

Note its one blind spot, stated honestly: it shares the kernels' assumption about the broadcast rule,
so it cannot check that. The rule was confirmed by reading the forward instead — row_mapping.i1 is
the slot index (the in-code comment calling it "selected expert index" is wrong), so the forward's
b column is slot % ne_b1.