Beckmann Flow Boundary Regularizer

Integrate this at the loss level for a model whose input contains n binary features, bits, or a small discrete latent code. For each minibatch point x, evaluate the scalar logit f_theta(x), all one-bit neighbors f_theta(x^{\oplus i}), and therefore b(x)=Lf_theta(x). Maintain a small vector field V_i(x) for the sampled cube points. Solve the convex approximation min_V (1/|B|) sum_x sqrt(sum_i V_i(x)^2+epsilon) subject to div V(x)=b(x). A practical MVP uses 10-30 projected primal-dual iterations: initialize V=0 and dual variable u=0; update u <- u + eta_u(div(V)-b); update V <- V - eta_V grad_V[<u,div(V)> + lambda_V sum_x sqrt(||V(x)||_2^2+epsilon)]; then optionally project each local vector onto the corresponding L2 ball. Cache neighbor indices so the divergence is a sparse gather/scatter operation. During early experiments, use a penalty on constraint residual plus the group norm, rather than requiring exact feasibility: R_B=(1/|B|)sum_x sqrt(||V(x)||_2^2+epsilon)+rho(div(V)-b)^2. Add lambda_B R_B to cross-entropy. The mathematical object computed exactly is the discrete Laplacian target b and the flow objective; the finite-batch flow and primal-dual iterate are empirical approximations. First test on a 10-20 dimensional synthetic parity-plus-noise task and a tabular binary UCI dataset, comparing cross-entropy plus ordinary edge-sensitivity against cross-entropy plus Beckmann regularization at matched FLOPs. Measure test accuracy, prediction changes under one-bit and multi-bit perturbations, flow constraint residual, and loss curvature. Success is lower multi-bit adversarial sensitivity or improved test accuracy at equal one-bit sensitivity, without a large training-time increase.