Positive Bidiagonal Diffusion Layer

Integrate this as a structured transition layer inside an SSM, graph-neural-network propagation block, or residual diffusion module. The input is a feature matrix X in R^{n x d}, where n is the ordered spatial/state dimension and d is the channel dimension; apply the same transition to every channel. Store unconstrained logits alpha[k,i] and beta[k,i], and set a[k,i]=sigmoid(alpha[k,i]) and b[k,i]=sigmoid(beta[k,i]). For each forward pass, execute the factors sequentially rather than forming T: initialize H=X; for k=1,...,p, compute H_i <- (1-a[k,i])H_i + a[k,i]H_{i-1}; then for k=1,...,q in the prescribed reverse product order, compute H_i <- (1-b[k,i])H_i + b[k,i]H_{i+1}. Use zero, reflecting, or periodic boundary handling, but keep boundary rows normalized. Optionally add a residual gate y=(1-g)X+gH with g=sigmoid(gamma). Each factor costs O(nd), so the full layer costs O((p+q)nd), versus O(n(p+q) d) for explicitly applying a band matrix and O(n^2d) for an unrestricted matrix. The mathematical quantities computed exactly are the bidiagonal updates, nonnegativity, row-stochasticity, and effective bandwidth; no spectral estimator is required. Initialize alpha and beta near logit(0.25) so propagation is initially local and stable. First test on length-256 copy/reconstruction and sequential-MNIST tasks using a linear SSM or 1D residual baseline with an unconstrained banded transition of the same nominal bandwidth. Compare validation accuracy, gradient norm growth, parameter count, and wall-clock time. Success means no exploding activations over 100+ recurrent steps, lower memory than the dense baseline, and equal or better accuracy at matched parameter count; a stronger result is faster loss descent at equal FLOPs.