Projected Contact Momentum Optimizer

Integrate this at the optimizer level after computing a minibatch loss and gradient. Store parameters $\theta$, momentum $p$, and two temporary copies $(\theta_a,p_a)$ and $(\theta_b,p_b)$. Use a positive-definite metric $M(\theta)$, initially block-diagonal or diagonal-plus-low-rank, so computing $M(\theta)^{-1}p$ is affordable. Set $\rho=\exp(-\gamma h/2)$ and apply the first exact damping half-step to both momenta: $p_a\leftarrow\rho p_a$ and $p_b\leftarrow\rho p_b$. For copy a, use kick-drift-kick: compute $g_a=\nabla L(\theta_a)$, update $p_a\leftarrow p_a-(h/2)g_a$, solve $M(\theta_a)v_a=p_a$, update $\theta_a\leftarrow\theta_a+h v_a$, recompute $g_a$, and apply the second kick. For copy b, use the reverse drift-kick-drift ordering: solve $M(\theta_b)v_b=p_b$, update $\theta_b\leftarrow\theta_b+(h/2)v_b$, compute the gradient, update $p_b\leftarrow p_b-hg_b$, solve again, and finish the remaining half drift. Average the endpoints with $\theta=(\theta_a+\theta_b)/2$ and $p=(p_a+p_b)/2$, overwrite both copies with this average, and apply the second damping half-step $p\leftarrow\rho p$. Use the paper's projection and $O(h^3)$ symmetry principle directly; estimate metric-vector products using conjugate gradients, Cholesky for small blocks, or a Woodbury solve for low-rank updates. The first experiment should train a two-layer MLP on MNIST or CIFAR-10, comparing AdamW, standard momentum, and this optimizer at equal gradient evaluations. Test several step sizes and friction values, including a regime where Euler damping $p\leftarrow(1-\gamma h)p$ becomes unstable or negative. Measure training loss versus optimizer steps, gradient evaluations, wall-clock cost, exploding-run frequency, and validation accuracy. The key signal is lower loss at equal steps and substantially fewer divergent runs at large $h$; ablate projection, reversed split ordering, and exact exponential damping separately.