GLM Energy-Preserving Feature Block

Integrate this as a replacement for a residual CNN block or as one step of a continuous-depth state-space model. First maintain four feature groups at each spatial location: T with C channels, J with dC channels, psi with dC channels, and phi with C channels; concatenate them into u. Parameterize each spatial operator with a convolution W_i and enforce skew-adjointness by A_i=(W_i-W_i^T)/2, where W_i^T is the adjoint convolution under the selected padding and boundary convention. Construct small fixed or learned symmetric field-coupling matrices B_i implementing the c_0 gradient/divergence and c_h curl couplings, then compute K u=-sum_i(B_i tensor A_i)u. Use the Cayley update by solving (I-hK/2)v=(I+hK/2)u with several GMRES iterations, or use an exact dense/block solve when the channel dimension is small. For damping, apply v <- v-h lambda Rv with diagonal R that selectively damps psi and phi. The paper determines the symmetric coupling and energy invariant; learned operators, h, lambda, and R are empirical choices. Measure norm drift ||v||^2/||u||^2 and residual-solve error. The first experiment should compare a 32-layer CIFAR-10 CNN against standard residual, antisymmetric, and orthogonal/Cayley blocks at matched parameters and FLOPs, testing h in {0.25,0.5,1,2}. Record loss descent, gradient norms, accuracy, norm drift, and NaN rate. Success means stable training at h=1 or 2, energy drift below 1e-4 per reversible block, and equal-or-better accuracy than the baselines.