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 +}