Puiseux Arclength Continuation for Implicit Layers

Integrate this at the implicit-state solve of an equilibrium layer, for example z=phi_theta(z,x,t) with residual h(z,t)=z-phi_theta(z,x,t); t interpolates a stable initialization or regularized problem at t=1 toward the target problem at small t. Maintain the last four accepted pairs (t_j,z_j). Before each step, estimate a local exponent by fitting log(||z_j-z_{j-1}||+epsilon)=b+alpha log(t_j+epsilon) with least squares over the recent differences. Clamp alpha to [0.1,4], reject unstable fits, and retain up to three distinct exponents alpha_i. Solve the displayed coefficient system using recent states and finite-difference tangents dot z_j=(z_j-z_{j-1})/(t_j-t_{j-1}); with more history than equations, use regularized least squares. Predict at t_new=t_j-delta_t using z_pred=a_0+sum_i a_i t_new^(alpha_i). Correct jointly in (z,t) with Newton or damped least squares on h(z,t)=0 plus the pseudo-arclength equation q_z^T(z-z_pred)+q_t(t-t_new)=0. Compute h_z and h_t with automatic differentiation, but estimate exponents and coefficients empirically; the paper contributes the fractional-power model and derivative scaling rather than neural-specific hyperparameters. Use step rejection when Newton fails, and adapt delta_t by doubling after success or halving after failure. First test a small 2- or 3-layer DEQ on MNIST or CIFAR-10, comparing this method with linear prediction plus fixed-parameter Newton over a homotopy from t=1 to 10^-4. Record failed solves, Newton iterations, accepted step sizes, wall-clock time, final residual, and task accuracy. A positive result is at least 1.5x larger average continuation steps near rank deficiency and fewer failed corrections at equal residual tolerance.