Matching-guided block Newton solver

Apply this to the nonlinear solve inside DEQ-style layers, equilibrium GNNs, or recurrent models trained with implicit differentiation. Define residuals F_theta(z,x) and obtain a structural incidence matrix A by tracing autodiff dependencies or by using the known variable list consumed by each equation MLP. At initialization, run a maximum bipartite matching on A. Build an undirected graph containing every matched edge and every alternating unmatched edge; compute connected components, then topologically order the components using edges that point from an upstream variable or equation block to a downstream block under the paper's partial-orientation rule. During each solver iteration, evaluate all residuals, then process blocks in that order. For block b, extract z_b and residual F_b, form J_b with autodiff Jacobian-vector products, solve (J_b+lambda I)Delta_b=-F_b using a small dense solve for blocks up to 64 variables or conjugate gradients for larger blocks, and update z_b<-z_b+eta Delta_b. After one sweep, optionally apply one global residual correction to account for ignored off-block Jacobian terms. Recompute the matching only if the architecture's sparsity changes; do not recompute it per batch. The paper provides the structural decomposition principle, while Jacobian values, damping lambda, and the stopping criterion are estimated computationally. Test first on synthetic sparse equilibrium systems and a 2-layer implicit GNN, comparing dense Newton, Broyden, and block Newton at equal residual tolerance. Measure wall-clock time, peak memory, solver iterations, residual norm, and backward-pass failures. The expected signal is 1.5–3x lower solve cost and fewer divergence events when the Jacobian is sparse and feedback clusters are small, with no increase in final task loss.