From 4cdbe30773d4b21efb00cb83a9244200158f2c14 Mon Sep 17 00:00:00 2001 From: godofecht Date: Fri, 21 Aug 2026 12:49:29 +0100 Subject: [PATCH] fix: centre multinomial softmax parameters to the sum-to-zero gauge lin433_multinomial_fit now subtracts the class mean from every feature column and from the intercept vector before returning, so the model it reports is the sum-to-zero representative scikit-learn reports. The means are accumulated in f64 before the f32 store. The change is small on purpose. theta starts at zero and the softmax cross-entropy gradient sums to zero across classes, so every LBFGS iterate already lies in the sum-to-zero subspace in exact arithmetic. What the centring removes is f32 accumulation drift out of that subspace: measured at 1.1e-05 on digits and 1.0e-06 on iris before, 1.4e-07 and 7.5e-08 after. It also makes the gauge a property of the returned model rather than of the starting point. Predictions do not move. On the canonical benchmark both LogisticRegression accuracies are identical (0.933333337 iris, 0.975000024 digits) and 160 of the 165 RESULT/DETAIL records are byte-identical; the five that move are LogisticRegression coefficient records at the 1e-06 level. On a direct before/after comparison over the test sets, zero of 30 iris and zero of 360 digits predicted labels change, and the largest softmax probability movement is 8.3e-07. tests/test_multinomial_gauge.flow pins the gauge, the invariance and the recentring. Verified to bite by deleting the centring block and rerunning. Co-Authored-By: Claude Opus 5 (1M context) --- lib/scikit/linear.flow | 45 ++++++ tests/test_multinomial_gauge.flow | 249 ++++++++++++++++++++++++++++++ 2 files changed, 294 insertions(+) create mode 100644 tests/test_multinomial_gauge.flow diff --git a/lib/scikit/linear.flow b/lib/scikit/linear.flow index feda9e5..cb3c7cd 100644 --- a/lib/scikit/linear.flow +++ b/lib/scikit/linear.flow @@ -1166,6 +1166,51 @@ function lin433_multinomial_fit(X: Matrix, y: ptr, classes: ptr, nc: i all_biases[c] = theta[base + n] } + # Issue #470: pin the softmax gauge to sum-to-zero across classes. + # + # Softmax scores are invariant to adding the same constant to every class, + # so (W, b) and (W + 1 c^T, b + d 1) describe the same model. The + # parameters are identified only up to that shift. scikit-learn reports the + # sum-to-zero representative, so anything comparing coefficients against it + # needs Flow to report the same one. + # + # This loop is close to a no-op in practice, and deliberately so. theta + # starts at zero, and the gradient of the softmax cross-entropy sums to + # zero across classes for every feature and for the intercept, because + # sum_c (p_ic - t_ic) = 0. The L2 term is linear in theta and so preserves + # that. Every LBFGS direction is a linear combination of gradients and of + # earlier steps, so in exact arithmetic each iterate stays in the + # sum-to-zero subspace it started in. What is left is f32 accumulation + # drift: measured at 1.1e-05 on digits and 1.0e-06 on iris, against + # coefficients up to 1.83 and 2.49. Centring removes that drift and makes + # the gauge a property of the returned model rather than of the starting + # point, which is what a warm start or a different initialisation would + # otherwise break. + # + # The means are accumulated in f64 before the f32 store, for the reason + # PR #443 moved its line-search objective to f64. + if nc > 0 { + let inv_nc: f64 = 1.0 / (nc as f64) + for j in 0 to n { + let mut w_total: f64 = 0.0 + for c in 0 to nc { + w_total = w_total + (all_weights[c * n + j]) as f64 + } + let w_mean: f64 = w_total * inv_nc + for c in 0 to nc { + all_weights[c * n + j] = (((all_weights[c * n + j]) as f64) - w_mean) as f32 + } + } + let mut b_total: f64 = 0.0 + for c in 0 to nc { + b_total = b_total + (all_biases[c]) as f64 + } + let b_mean: f64 = b_total * inv_nc + for c in 0 to nc { + all_biases[c] = (((all_biases[c]) as f64) - b_mean) as f32 + } + } + free(X_aug as ptr) free(tgt as ptr) free(Z as ptr) diff --git a/tests/test_multinomial_gauge.flow b/tests/test_multinomial_gauge.flow new file mode 100644 index 0000000..dd9a51c --- /dev/null +++ b/tests/test_multinomial_gauge.flow @@ -0,0 +1,249 @@ +# Issue #470: the multinomial softmax gauge. +# +# Softmax scores are invariant to adding the same constant to every class, so +# (W, b) and (W + 1 c^T, b + d 1) are the same model written two ways. The +# parameters are identified only up to that shift. scikit-learn reports the +# sum-to-zero representative; lin433_multinomial_fit now does too. +# +# Three checks: +# 1. A fitted multinomial model sums to zero across classes, in every feature +# column and in the intercept vector. The bound is a few float32 ulps of +# the largest coefficient, which is all the f32 store can leave behind +# once the means are removed in f64. +# 2. The invariance itself: shifting a fitted model by an arbitrary per-class +# constant leaves every softmax probability and every predicted label +# unchanged. This is what makes check 1 free to impose. +# 3. Removing the class mean from the shifted model returns it to the fitted +# parameters, so centring picks out one representative of the family and +# always the same one. +# +# The fixture is 8 classes over 32 features, large enough that the LBFGS run +# accumulates visible f32 drift out of the zero-sum subspace. Without the +# centring the fit lands about an order of magnitude further off the gauge, and +# check 1 fails: verified by deleting the centring block and rerunning. + +import "lib/scikit/scikit.flow" + +extern { + function printf(fmt: string, ...) -> i32 +} + +function mgz470_abs(x: f32) -> f32 { + if x < 0.0 { return -x } + return x +} + +# Softmax of one score row in place. exp() lives here because the transpiler +# miscompiles it inside a nested while loop (compiler bug #421). +function mgz470_softmax_row(row: ptr, nc: i32) -> void { + let mut max_val: f32 = row[0] + for c in 1 to nc { + if row[c] > max_val { max_val = row[c] } + } + let mut total: f32 = 0.0 + for c in 0 to nc { + let e: f32 = exp((row[c] - max_val) as f64) as f32 + row[c] = e + total = total + e + } + if total > 0.0 { + let inv: f32 = 1.0 / total + for c in 0 to nc { row[c] = row[c] * inv } + } +} + +# Softmax probabilities of an explicit (W, b) against X, row-major rows x nc. +function mgz470_probs(W: ptr, b: ptr, nc: i32, nf: i32, X: Matrix, out: ptr) -> void { + for i in 0 to X.rows { + for c in 0 to nc { + let mut z: f32 = b[c] + for j in 0 to nf { + z = z + W[c * nf + j] * matrix_at(X, i, j) + } + out[i * nc + c] = z + } + mgz470_softmax_row(out + i * nc, nc) + } +} + +function main() -> i32 { + let nc: i32 = 8 + let nf: i32 = 32 + let m: i32 = 400 + + let X: Matrix = matrix_new(m, nf) + let y: ptr = array_new_f32(m) + for i in 0 to m { + let cls: i32 = i % nc + y[i] = cls as f32 + for j in 0 to nf { + let mut v: f32 = ((((i * 31 + j * 7 + cls * 13) % 29) as f32) / 14.0) - 1.0 + if j % nc == cls { v = v + 1.5 } + matrix_set(X, i, j, v) + } + } + + let model: LogisticRegression = logistic_regression_fit(X, y, nc, 200, 0.1, penalty_l2_from_c(1.0, m)) + if model.n_classes != nc { + println("FAIL: multinomial fit did not find eight classes") + return 1 + } + + # ---- 1. the fitted parameters sum to zero across classes ---- + let mut max_coef: f32 = 0.0 + for k in 0 to nc * nf { + let a: f32 = mgz470_abs(model.weights[k]) + if a > max_coef { max_coef = a } + } + let mut max_intercept: f32 = 0.0 + for c in 0 to nc { + let a: f32 = mgz470_abs(model.biases[c]) + if a > max_intercept { max_intercept = a } + } + + # Rounding one already-centred f64 value to f32 costs at most half an ulp, + # and nc of them are summed, so the residual cannot exceed + # nc * ulp(max) / 2, with 2^-23 the float32 mantissa step. The bound below + # is twice that, which leaves the check room for summation order and still + # bites: without the centring this fixture reports 1.74e-06 against a + # coefficient bound of 1.10e-06 and 6.24e-07 against an intercept bound of + # 6.92e-08. + let coef_bound: f32 = (nc as f32) * max_coef * 0.00000011920929 + let intercept_bound: f32 = (nc as f32) * max_intercept * 0.00000011920929 + + let mut worst_coef_sum: f32 = 0.0 + for j in 0 to nf { + let mut s: f64 = 0.0 + for c in 0 to nc { + s = s + (model.weights[c * nf + j]) as f64 + } + let a: f32 = mgz470_abs(s as f32) + if a > worst_coef_sum { worst_coef_sum = a } + } + let mut bs: f64 = 0.0 + for c in 0 to nc { + bs = bs + (model.biases[c]) as f64 + } + let intercept_sum: f32 = mgz470_abs(bs as f32) + + printf(" max coefficient %.7f, worst class sum %.9g, bound %.9g\n", max_coef, worst_coef_sum, coef_bound) + printf(" max intercept %.7f, class sum %.9g, bound %.9g\n", max_intercept, intercept_sum, intercept_bound) + + if worst_coef_sum > coef_bound { + println("FAIL: fitted coefficients are not centred across classes") + return 1 + } + if intercept_sum > intercept_bound { + println("FAIL: fitted intercepts are not centred across classes") + return 1 + } + println(" OK: the fitted multinomial parameters sum to zero across classes") + + # ---- 2. an arbitrary per-class shift changes nothing observable ---- + let shifted_w: ptr = array_new_f32(nc * nf) + let shifted_b: ptr = array_new_f32(nc) + for j in 0 to nf { + let shift: f32 = 0.25 + 0.125 * ((j % 5) as f32) + for c in 0 to nc { + shifted_w[c * nf + j] = model.weights[c * nf + j] + shift + } + } + for c in 0 to nc { + shifted_b[c] = model.biases[c] + 0.75 + } + + let probs_fit: ptr = array_new_f32(m * nc) + let probs_shifted: ptr = array_new_f32(m * nc) + mgz470_probs(model.weights, model.biases, nc, nf, X, probs_fit) + mgz470_probs(shifted_w, shifted_b, nc, nf, X, probs_shifted) + + let mut worst_prob: f32 = 0.0 + for k in 0 to m * nc { + let d: f32 = mgz470_abs(probs_fit[k] - probs_shifted[k]) + if d > worst_prob { worst_prob = d } + } + + let pred_fit: ptr = logistic_predict(model, X) + let shifted_model: LogisticRegression = LogisticRegression { + weights: shifted_w, + bias: shifted_b[0], + biases: shifted_b, + n_features: nf, + classes: model.classes, + n_classes: nc, + penalty: model.penalty, + fitted: true, + n_iter: model.n_iter + } + let pred_shifted: ptr = logistic_predict(shifted_model, X) + let mut flips: i32 = 0 + for i in 0 to m { + if pred_fit[i] != pred_shifted[i] { flips = flips + 1 } + } + + printf(" shifted model: worst probability gap %.9g, label flips %d\n", worst_prob, flips) + if worst_prob > 0.00002 { + println("FAIL: a per-class shift moved the softmax probabilities") + return 1 + } + if flips != 0 { + println("FAIL: a per-class shift moved a predicted label") + return 1 + } + println(" OK: a per-class shift leaves probabilities and labels unchanged") + + # ---- 3. centring the shifted model returns the fitted parameters ---- + let inv_nc: f64 = 1.0 / (nc as f64) + for j in 0 to nf { + let mut s: f64 = 0.0 + for c in 0 to nc { + s = s + (shifted_w[c * nf + j]) as f64 + } + let mean: f64 = s * inv_nc + for c in 0 to nc { + shifted_w[c * nf + j] = (((shifted_w[c * nf + j]) as f64) - mean) as f32 + } + } + let mut sb: f64 = 0.0 + for c in 0 to nc { + sb = sb + (shifted_b[c]) as f64 + } + let mean_b: f64 = sb * inv_nc + for c in 0 to nc { + shifted_b[c] = (((shifted_b[c]) as f64) - mean_b) as f32 + } + + let mut worst_w: f32 = 0.0 + for k in 0 to nc * nf { + let d: f32 = mgz470_abs(shifted_w[k] - model.weights[k]) + if d > worst_w { worst_w = d } + } + let mut worst_b: f32 = 0.0 + for c in 0 to nc { + let d: f32 = mgz470_abs(shifted_b[c] - model.biases[c]) + if d > worst_b { worst_b = d } + } + printf(" recentred: worst coefficient gap %.9g, worst intercept gap %.9g\n", worst_w, worst_b) + if worst_w > 0.00002 { + println("FAIL: centring the shifted model did not restore the fitted coefficients") + return 1 + } + if worst_b > 0.00002 { + println("FAIL: centring the shifted model did not restore the fitted intercepts") + return 1 + } + println(" OK: centring recovers the fitted parameters from a shifted model") + + array_free_f32(probs_fit) + array_free_f32(probs_shifted) + array_free_f32(pred_fit) + array_free_f32(pred_shifted) + array_free_f32(shifted_w) + array_free_f32(shifted_b) + array_free_f32(y) + matrix_free(X) + logistic_regression_free(model) + + println("All issue #470 multinomial gauge tests passed!") + return 0 +}