Jacobian-Free Secant MPC for Learned Dynamics

Integrate this at the planning layer of a model-based RL agent, immediately after the learned dynamics model $f_\theta(x,u)$ and before action selection; do not alter the neural model or backpropagate through planning iterations. Given the current state $x_0$, a warm-start control sequence $u[0:H-1]$, and reference states/actions, repeat the following loop for one to three iterations. First roll out $f_\theta$ to obtain $\bar{x}[0:H]$. For every stage $j$, construct $A[j]$ and $B[j]$ by coordinate telescoping. Start with $z_{prev}=0$ and $y_{prev}=f_\theta(z_{prev})$; for each coordinate $q$ of $z=(\bar{x}[j],u[j])$, change only coordinate $q$ by $\delta_q=\operatorname{sign}(z_q)\max(|z_q|,10^{-3})$, evaluate $y_{new}=f_\theta(z_{new})$, and set column $q$ to $(y_{new}-y_{prev})/\delta_q$. Then set $z_{prev}=z_{new}$ and $y_{prev}=y_{new}$. Split the columns into $A[j]$ and $B[j]$. Solve the banded quadratic program with constraints $x[j+1]=A[j]x[j]+B[j]u[j]$, action bounds, and optional state bounds. Apply only the first action, shift the solution for warm-starting, and re-roll out the nonlinear model. The secant matrices are computed from function evaluations; no neural Jacobian is required. Estimate local contraction empirically with $\hat{\rho}=\|u_{new}-u_{old}\|/\max(\|u_{old}-u_{older}\|,10^{-8})$ and stop when the sequence change is below tolerance or $\hat{\rho}\geq 1$. Use Riccati recursion or a sparse QP solver rather than a dense generic SQP solver. First test on a 4D or 8D learned cartpole or pendulum world model with horizon 32, comparing against shooting MPC with autodiff Jacobians and CEM at equal model evaluations. Record wall-clock planning time, constraint violations, return, and robustness to model noise. Success means at least 1.5x lower planning latency at matched return, fewer gradient-related numerical failures, and stable performance when the horizon is doubled.