SBP-Factorized Dissipative Residual Mixing

Integrate the operator into the token-mixing part of a residual sequence layer or into the linear state update of an SSM. Let H have shape batch x sequence_length x hidden_width, let W = diag(w_0,...,w_{T-1}) be a positive diagonal token-energy matrix, and let B be a banded sequence difference matrix. Initialize B as a second-difference operator: (BH)_t = H_t - 2H_{t-1} + H_{t-2} for interior positions. Use one-sided rows at the first and last two positions; a fourth-order version can use the stencil [1,-4,6,-4,1]. Compute the operator without forming dense matrices: BH = banded_apply(B,H); Z = W * BH; QH = -a * inverse(W) * banded_apply(B.T,Z). Update the layer as H_out = H + eta * (F_theta(H) + QH), where F_theta is the original attention, MLP, or SSM branch, eta is the residual step size, and a is a nonnegative scalar. Set W to the identity initially; optionally use w_t = 1 + lambda for padding positions to suppress boundary artifacts. If B is learned, learn only its band coefficients and retain Q = -a W^{-1} B^T W B rather than learning Q directly. The mathematical guarantee is the nonpositive weighted quadratic form; any W or B adaptation must be checked numerically by evaluating tr(H^T W QH) / ||H||_F^2. First test a six-layer causal Transformer or diagonal SSM on WikiText-2 and a synthetic long-range copy task. Compare the baseline, unconstrained fourth-order smoothing, and factorized SBP mixing at matched FLOPs. Record loss at equal updates, gradient-norm spikes, weighted activation energy, validation perplexity, and performance at sequences four times longer than training. Success means fewer exploding-gradient events and improved long-context accuracy without an appreciable loss increase at the training length.