The Derived Optimizer: Six Questions and Failure on Real Architectures

1. Motivation

The starting point is the learnable diagonal variant of the induced-metric optimizer with curvature correction. The update rule for the per-coordinate log-scale $s_i$ that this optimizer used at the start of the program was

$$\dot s_i = \mu\bigl[\xi e^{s_i} g_i^2 - \beta_c \tanh(\hat H_i/\tau) e^{s_i} g_i^2\bigr] - \mu\lambda s_i,$$

followed by mean-centering ($\sum_i s_i = 0$) and clipping ($|s_i| \le c$). Here $g_i = \partial_i L$ is the gradient component, $\hat H_i$ is a cheap estimate of the diagonal Hessian entry, and the symbols $\mu, \xi, \beta_c, \tau, \lambda, c$ are hyperparameters. Adding the global denominator $r = 1/(1 + |\hat v|)$ from the smooth-clipping mechanism, the parameter update is $\theta_i \leftarrow \theta_i - \eta\,r\,e^{s_i}\,\hat m_i$ with $\hat m$ the bias-corrected momentum. The complete hyperparameter list is $\{\eta, \mu, \lambda, \xi, \beta_c, \tau, c, \beta, \beta_m, \varepsilon_{\text{sec}}\}$, ten in total. The motivation for the diagonal preconditioner class itself comes from the hierarchy of curvature corrections: the diagonal learnable metric is the lowest level of the hierarchy with enough freedom to act as a Newton preconditioner in principle.

Four structural problems with the heuristic rule were visible by inspection.

First, the metric explodes. The driving force $\xi e^{s_i} g_i^2$ grows exponentially in $s_i$ while the regulariser $\lambda s_i$ grows linearly. The equilibrium of the two requires $\xi g^2 / \lambda < 1/e$ (the existence condition for the Lambert $W$ function). In practice $\xi g^2 / \lambda$ is in the range $10^3$ to $10^4$, so the equilibrium does not exist and only clipping prevents divergence. The clip is doing structural work, not safety work.

Second, the curvature correction is anti-Newton. Without the $\tanh$ term, $s_i$ relaxes to $s_i^* \propto \xi g_i^2 / \lambda$. On a quadratic loss $L = \tfrac{1}{2}\sum h_i \theta_i^2$, where $g_i = h_i \theta_i$, this gives $s_i^* \propto h_i^2 \theta_i^2$. Directions with large curvature receive large metric, which is the opposite of Newton preconditioning. The Pearson correlation between the learned $s_i$ and $H_{ii}$, measured on a 10-dimensional quadratic with random condition number 100, comes out to $+0.71$. Newton would give $-1.0$.

Third, the role of the denominator $r$ is not understood. It clearly damps something when gradients are large, but its quantitative effect on the per-coordinate dynamics and its interaction with the per-coordinate metric are unclear.

Fourth, ten hyperparameters is a lot. Several of them interact in non-obvious ways. The optimal value of $\xi$ varies by four to six orders of magnitude across tasks.

The program below addresses each problem in turn by asking one question with a clean answer. The structure of this entry is to ask six questions and derive the corresponding piece of the algorithm from first principles, then assemble the pieces into a final algorithm with three free hyperparameters, then examine why the algorithm nonetheless fails on real neural networks.

2. Question 1: why does the metric explode?

The metric scale is $M_i = e^{s_i}$. The quantity we want to control is the denominator trace $v = \xi \sum_j M_j g_j^2$, which is linear in $M$,

$$\frac{\partial v}{\partial M_i} = \xi g_i^2.$$

This is bounded and well-behaved. The optimizer, however, updates $s = \log M$, not $M$ itself, and the chain rule introduces

$$\frac{\partial v}{\partial s_i} = \xi g_i^2 \cdot e^{s_i}.$$

The factor $e^{s_i}$ is the entire source of instability. Updating in the log parameterisation, the derivative we use accidentally inherits the exponential factor that comes from the chain rule, even though the underlying quantity $v$ in the $M$-parameterisation has no exponential dependence.

2.1 The fix

Update $s$ using the gradient with respect to $M$, not with respect to $s$. This is mirror descent in $M$-space with the log barrier, or equivalently the natural gradient for the metric parameterisation. (Natural gradient here means: rescale the ordinary gradient by the inverse of the Fisher metric on the parameter manifold, which in this case is just division by $M$ to undo the exponential Jacobian.) The corrected update is

$$\dot s_i = \mu[\xi g_i^2 - \lambda s_i].$$

The $e^{s_i}$ factor is gone. In code this is one line.

2.2 Stability proof

For fixed gradients, the equilibrium $s_i^* = \xi g_i^2 / \lambda$ always exists. To prove the equilibrium is globally asymptotically stable, take the Lyapunov candidate $V(s) = \tfrac{1}{2} \|s - s^*\|^2$. A Lyapunov function is a scalar function of the state that decreases monotonically along the dynamics, with its unique minimum at the equilibrium; the existence of such a function is the standard certificate for global asymptotic stability. Differentiate $V$ along the dynamics,

$$\dot V = (s - s^*)^\top \dot s = (s - s^*)^\top \mu[\xi g^2 - \lambda s] = -\mu \lambda (s - s^*)^\top (s - s^*) = -2 \mu \lambda V.$$

The substitution uses $\xi g^2 = \lambda s^*$ from the equilibrium definition. So $V$ decreases exponentially at rate $2 \mu \lambda$, regardless of initial conditions and regardless of $\xi$ or $\lambda$ individually. Global asymptotic stability follows without any linearisation.

2.3 Mean-centered version

In the actual optimizer, the gauge symmetry $s \mapsto s + c\,\mathbf{1}$ is fixed by mean-centering: we enforce $\sum_i s_i = 0$. With the constraint added via a Lagrange multiplier $\mu_L$, the equilibrium conditions become

$$\xi g_i^2 - \lambda s_i = \mu_L \quad \text{for all } i, \qquad \sum_i s_i = 0.$$

Solving the system gives $\mu_L = \xi \overline{g^2}$ where $\overline{g^2} = \tfrac{1}{N}\sum_j g_j^2$, and

$$s_i^* = \frac{\xi}{\lambda}\,(g_i^2 - \overline{g^2}).$$

The constrained equilibrium always exists and is always stable (the Hessian of $V$ on the constraint surface remains positive-definite). The Lagrange multiplier is the per-component bias that the constraint introduces.

2.4 Consequence for clipping

Before the fix, clipping was structurally required because the unconstrained dynamics had no equilibrium. After the fix, clipping is optional. It becomes a soft bound on the condition number of the preconditioner ($\kappa_M \le e^{2c}$) rather than the only thing preventing divergence. In practice we keep a modest clip ($c \approx 3$) as a safeguard for noisy stochastic gradients.

3. The separation principle

The Q1 fix has a generalisation that is useful for the rest of the program. The general form of the update rule is

$$\dot s_i = \phi(s_i)\,\bigl[D_i - \lambda s_i\bigr],$$

where $\phi$ is the dynamics modulator and $D$ is the target driver. At equilibrium (assuming $\phi$ does not vanish), the equilibrium satisfies $D_i = \lambda s_i^*$, so $s_i^*$ depends only on $D$ and not on $\phi$. This is the separation principle: the rate at which the metric converges (governed by $\phi$) and the place it converges to (governed by $D$) are independent design choices.

The original rule corresponds to $\phi = e^s$ (unstable) and $D = \xi g^2$ (anti-Newton). The Q1 fix sets $\phi = 1$ while leaving $D = \xi g^2$ unchanged.

3.1 Why $\phi = 1$ is the optimal choice for $\phi$

Two arguments single out $\phi = 1$ within the class of $\phi$ that preserve the same equilibrium.

A minimax convergence-rate argument. Non-uniform $\phi$ produces non-uniform per-component convergence rates. The worst-case rate dominates the overall convergence. Maintaining stability with non-uniform $\phi$ forces a reduction of $\lambda$ to accommodate the fastest component, which slows every other component. The choice that maximises the worst-case rate is uniform $\phi$, which after normalisation is $\phi = 1$.

A noise-robustness argument. Under stochastic gradients, the variance of $\dot s_i$ scales as $\phi(s_i)^2 \cdot \mathrm{Var}(D_i)$. Any $\phi > 1$ amplifies the noise injected into $s$. The choice that minimises noise amplification is $\phi = 1$.

The dynamics question is fully resolved by $\phi = 1$. The remaining question is what $D$ should be.

4. Question 2: what should the metric target?

On a local quadratic model $L = \tfrac{1}{2} \theta^\top \mathrm{diag}(h)\,\theta$ with $h_i > 0$, the preconditioned gradient descent update is $\theta_i \leftarrow \theta_i - \alpha\,e^{s_i}\,g_i$. Substituting $g_i = h_i \theta_i$ gives $\theta_i \leftarrow (1 - \alpha e^{s_i} h_i) \theta_i$, so the per-coordinate contraction factor is

$$\rho_i = |1 - \alpha\,e^{s_i}\,h_i|.$$

The overall convergence rate is the maximum $\rho_i$ across coordinates. The objective is to choose $s$ minimising this maximum.

4.1 Theorem (Newton minimax preconditioner)

Among diagonal preconditioners $M = \mathrm{diag}(e^{s_i})$ with $\sum_i s_i = 0$,

(a) the minimax convergence rate $\min_s \max_i |1 - \alpha\,e^{s_i}\,h_i|$ is achieved uniquely at

$$s_i^* = -\log h_i + \overline{\log h},$$

where $\overline{\log h} = \tfrac{1}{N}\sum_j \log h_j$;

(b) at this optimum, all per-coordinate preconditioned eigenvalues $e^{s_i}\,h_i$ equal $h_{\mathrm{geo}}$, where $h_{\mathrm{geo}} = \exp(\overline{\log h})$ is the geometric mean of the curvatures (the contraction factors $|1 - \alpha\,e^{s_i}\,h_i|$ are then all equal to $|1 - \alpha h_{\mathrm{geo}}|$);

(c) with the additional optimal choice of learning rate $\alpha^* = 1/h_{\mathrm{geo}}$, the rate is zero (one-step convergence).

The term "minimax" here means the optimisation $\min_s \max_i$: we choose $s$ to make the worst-case contraction factor (over coordinates) as small as possible. The classical minimax solution equalises the worst-case quantity across all components, which is what produces the geometric-mean structure.

4.2 Proof

Let $p_i = \alpha\,e^{s_i}\,h_i$. Mean-centering enforces $\sum_i s_i = 0$, equivalently $\prod_i e^{s_i} = 1$, so

$$\prod_i p_i = \alpha^N \prod_i h_i$$

is independent of $s$. The problem reduces to minimising $\max_i |1 - p_i|$ subject to fixed product of $p_i$. The argument is an endpoint-equalisation: suppose the maximum of $|1 - p_i|$ over $i$ is attained both at the largest $p_i$ (call it $p^+$) and the smallest ($p^-$). If $p^- \to 0$ further, $|1 - p^-|$ grows; if $p^+ \to \infty$ further, $|1 - p^+|$ grows. So any spread that pushes $p^-$ down or $p^+$ up (while preserving the geometric mean by the product constraint) strictly increases $\max|1 - p_i|$. Conversely, contracting the spread toward the geometric mean weakly decreases the maximum. The optimum is therefore attained at $p_i = (\alpha^N \prod h_i)^{1/N} = \alpha h_{\mathrm{geo}}$ for all $i$, giving $e^{s_i} = h_{\mathrm{geo}}/h_i$ and (after taking logs and applying mean-centering) $s_i^* = -\log h_i + \overline{\log h}$. This proves (a). Part (b) is immediate from the definition of $p_i$ at the optimum. For (c), $\alpha^* = 1/h_{\mathrm{geo}}$ makes $p_i = 1$ for every $i$, giving $\rho_i = 0$.

(The naive "strict convexity of $\log p \mapsto |1 - e^{\log p}|$" version of this argument is false: the function $x \mapsto |1 - e^x|$ is concave for $x < 0$ and convex for $x > 0$, with a kink at $x = 0$. The endpoint argument above sidesteps the kink.)

4.3 Interpretation

The optimal metric is Newton preconditioning in log-space, with mean-centering as a gauge fix. Every preconditioned eigenvalue equals $h_{\mathrm{geo}}$. The condition number becomes 1. The landscape becomes isotropic to the optimiser. The mean-centering convention is what removes the overall scaling freedom and makes the answer unique; without it, any $s_i = -\log h_i + c$ for arbitrary $c$ would work, and the natural choice is $c = \overline{\log h}$.

4.4 Numerical comparison

On a 10-dimensional quadratic with condition number 100, run for 20{,}000 steps, the comparison between the original anti-Newton rule, the heuristic $\tanh$ rule that the curvature-aware variant added, and the Newton minimax target is striking.

Taking the unpreconditioned baseline at $4.7 \times 10^{-10}$ (Section 5.3):

Rule	Final loss	Ratio vs no preconditioning
$D_i = \xi g_i^2$ (anti-Newton)	$0.021$	$\approx 10^{7}$ times worse
$\tanh$ heuristic	$0.13$	$\approx 10^{8}$ times worse
Newton target $D_i = -\lambda \log h_i$	$1.3 \times 10^{-42}$	$\approx 10^{33}$ times better

The original rule is not merely suboptimal: it is actively worse than no preconditioning. The anti-Newton structure inflates the condition number from 100 to $1.6 \times 10^7$. The full numerical comparison, including loss-vs-step trajectories, is in notebooks/newton-minimax-comparison.nb.

4.5 The resulting metric update

By the separation principle, choose $D_i = -\lambda \log \max(\hat H_{ii}, \varepsilon)$, so the equilibrium satisfies $s_i^* = -\log H_{ii} + \mathrm{const}$ (after mean-centering). The discrete-time form of the update is

$$s_{i,t+1} = (1 - \mu\lambda)\,s_{i,t} - \mu\lambda\,\log \max(\hat H_{ii,t}, \varepsilon),$$

which is an exponential moving average (EMA) of $-\log H_{ii}$ with smoothing factor $\beta_s = \mu \lambda$. The floor $\varepsilon$ avoids the logarithm of zero or negative estimates of the diagonal Hessian. Three of the original hyperparameters ($\beta_c, \tau, \varepsilon_{\text{sec}}$) are absorbed or removed by this construction.

5. Question 3: how to estimate the Hessian diagonal?

Two candidate estimators: the secant approximation, and Hutchinson's stochastic estimator.

5.1 Secant estimator

Compare gradients at two consecutive iterates,

$$\hat h_i^{\text{sec}} = \frac{g_i^{(t)} - g_i^{(t-1)}}{\theta_i^{(t)} - \theta_i^{(t-1)}}.$$

This is free in the sense that it reuses gradients already computed. On a diagonal Hessian, it is exact. On a non-diagonal Hessian, expand $g_i^{(t)} = g_i^{(t-1)} + \sum_j H_{ij}\,\Delta\theta_j$ to get

$$\hat h_i^{\text{sec}} = H_{ii} + \sum_{j \ne i} H_{ij}\,\frac{\Delta\theta_j}{\Delta\theta_i}.$$

The cross terms are systematic bias, not noise. The bias does not shrink as the step size shrinks because the ratios $\Delta\theta_j / \Delta\theta_i$ are determined by the optimisation trajectory, not by the step size. On a 10-dimensional quadratic with modest off-diagonal structure, the initial maximum relative error across coordinates exceeds 1100% in numerical tests, and the smallest-curvature coordinate exceeds 730%.

5.2 Hutchinson estimator

Draw a Rademacher vector $v \in \{-1, +1\}^N$ (each component independently $\pm 1$ with equal probability). Compute the Hessian-vector product $Hv$ via reverse-mode automatic differentiation (one extra backward pass per step). Then set

$$\hat h_i = v_i\,(Hv)_i.$$

Take expectation over $v$: $\mathbb{E}[v_i v_j] = \delta_{ij}$ because the components are independent and have unit variance. So

$$\mathbb{E}[\hat h_i] = \mathbb{E}[v_i \sum_j H_{ij} v_j] = \sum_j H_{ij}\,\delta_{ij} = H_{ii}.$$

The estimator is exactly unbiased. The variance is

$$\mathrm{Var}(\hat h_i) = \mathbb{E}[\hat h_i^2] - H_{ii}^2 = \sum_{j \ne i} H_{ij}^2.$$

The Rademacher distribution is the variance-minimising choice within the class of zero-mean unit-variance distributions for this estimator; this is why it is preferred over Gaussian $v$ (which would give a larger variance constant).

5.3 Numerical comparison

On a 10-dimensional quadratic with condition number 100, run for 20{,}000 steps:

Estimator	Final loss	Ratio vs no preconditioning
No preconditioning	$4.7 \times 10^{-10}$	baseline
Secant + clip	$2.6 \times 10^{-12}$	$100$ times better
Hutchinson + clip	$1.3 \times 10^{-84}$	$10^{74}$ times better
Perfect Newton	$1.3 \times 10^{-84}$	$10^{74}$ times better

Hutchinson with EMA smoothing exactly matches perfect Newton. The factor-of-two per-step cost (one extra Hessian-vector product, comparable to one extra backward pass) buys a $10^{72}$-times improvement in final accuracy on this problem.

5.4 Robustness theorem

Suppose every $\hat h_i$ is within multiplicative factor $(1 \pm \delta)$ of $H_{ii}$, that is, $|\hat h_i / H_{ii} - 1| \le \delta$ for some $\delta \in [0, 1)$. Then the preconditioned condition number satisfies

$$\kappa_{\mathrm{prec}} \le \frac{1 + \delta}{1 - \delta}.$$

This bound is independent of the original condition number $\kappa(H)$ and independent of the dimension $N$.

The proof is short. The preconditioned spectrum has $e^{s_i^*} H_{ii} = h_{\mathrm{geo}}\,(H_{ii}/\hat H_{ii})$. The ratio $H_{ii}/\hat H_{ii} \in [1/(1+\delta), 1/(1-\delta)]$, so the preconditioned eigenvalues lie in $[h_{\mathrm{geo}}/(1+\delta), h_{\mathrm{geo}}/(1-\delta)]$, and the ratio is $(1+\delta)/(1-\delta)$. The full proof, including extension from multiplicative error to log-space error, is in Appendix B.

Numerical implications. Even with 50% estimation error ($\delta = 0.5$), $\kappa_{\mathrm{prec}} \le 3$. Even with 90% error ($\delta = 0.9$), $\kappa_{\mathrm{prec}} \le 19$, which on a problem with $\kappa(H) = 10^4$ is still a 530-times improvement.

The mechanism is the log-space parameterisation. A factor-of-$k$ multiplicative error in $\hat h_i$ shifts $s_i$ by $\log k$, which is small for moderate $k$. On a $\kappa = 10^4$ problem the target range of $s_i$ values is $\log(10^4) \approx 9$ nats, so a factor-of-10 error (about 2.3 nats) is a 25% perturbation of the relevant scale.

The secant estimator's problem is that its initial relative error exceeds $\delta = 1$ by an order of magnitude, placing it outside the domain where the robustness theorem applies. Hutchinson stays inside the domain from step 1.

6. Question 4: does the joint $(\theta, s)$ system converge?

Per-component theorems analyse $s$ assuming $\theta$ is frozen. In reality the two coordinates evolve simultaneously, and coupled dynamical systems can oscillate, cycle, or diverge even when each subsystem is stable in isolation. We need a joint convergence result.

6.1 On quadratics, the coupling vanishes

For $L = \tfrac{1}{2}\sum_i h_i \theta_i^2$, the Hessian diagonal is the constant $h_i$. The Newton target $s_i^* = -\log h_i + \overline{\log h}$ does not depend on $\theta$. The metric converges to a fixed point, not a moving one, and the $s$ dynamics are decoupled from the $\theta$ dynamics.

6.2 Joint Lyapunov function

Combine the energy in $\theta$ with the energy in $s$ relative to the Newton target,

$$V(\theta, s) = L(\theta) + \tfrac{\alpha}{2}\|s - s^*\|^2,$$

for some $\alpha > 0$ to be chosen. The time derivative along the joint dynamics is

$$\dot V = -\eta\,r \sum_i e^{s_i}\,h_i^2\,\theta_i^2 - \alpha\,\mu\,\lambda\,\|s - s^*\|^2.$$

Both terms are non-positive. The energy is strictly decreasing except at the joint equilibrium $\theta = 0$, $s = s^*$. This is global asymptotic stability of the joint system on quadratics, valid for any positive ratio $\mu/\eta$. No timescale separation is required.

We verified this numerically: zero monotonicity violations across 2{,}000 random initialisations over 20{,}000 steps each. Details in notebooks/joint-lyapunov-verification.nb.

6.3 Beyond quadratics

At any non-degenerate local minimum $\theta^*$ of a smooth (non-quadratic) loss, the joint Jacobian of the update is block-triangular,

$$J(\theta^*, s^*) = \begin{pmatrix} A & 0 \\ B & C \end{pmatrix}.$$

The block-zero in the upper right is because $\ell = 0$ at $\theta^*$, killing all cross-terms between the $\theta$ block and the $s$ block. Eigenvalues of a block-triangular matrix are the eigenvalues of the diagonal blocks. Both blocks contract under standard learning-rate conditions (the $\theta$ block is the usual contraction of preconditioned gradient descent, the $s$ block contracts at rate $1 - \mu\lambda$). The joint system therefore inherits local asymptotic stability from each subsystem, without requiring that the metric converge faster than the parameters.

For stochastic gradients with diminishing step sizes, Borkar's two-timescale stochastic-approximation theorem applies (Borkar, 1997), giving almost-sure convergence. With constant step sizes the analysis gives convergence to an $O(\eta + \mu)$ neighbourhood of the joint equilibrium. These guarantees are comparable in strength to the convergence results available for standard adaptive optimisers.

7. Question 5: what about the embedding weight $\xi$?

The smooth-clipping denominator $r = 1/(1 + |\hat v|)$ with $v = \xi \sum_i e^{s_i}\,g_i^2$ remains in the optimizer after Q1-Q4 have done their work. With Newton-targeted $s$, does the denominator still help?

7.1 The trace simplifies at equilibrium

At Newton equilibrium $e^{s_i^*} = h_{\mathrm{geo}}/h_i$, the denominator trace is

$$v = \xi \sum_i \frac{h_{\mathrm{geo}}}{h_i}\,(h_i \theta_i)^2 = \xi\,h_{\mathrm{geo}}\,\sum_i h_i\,\theta_i^2 = 2\,\xi\,h_{\mathrm{geo}}\,L.$$

The trace is proportional to the loss. This was verified symbolically (Mathematica, $N = 4$) and numerically (Python, $N = 20$, $\kappa = 1000$).

7.2 The denominator vanishes near the minimum

The effective per-coordinate convergence rate becomes

$$\eta_{\mathrm{eff}} = \frac{\eta\,h_{\mathrm{geo}}}{1 + 2\,\xi\,h_{\mathrm{geo}}\,L}.$$

As $L \to 0$, $\eta_{\mathrm{eff}} \to \eta\,h_{\mathrm{geo}}$. The asymptotic convergence rate is completely independent of $\xi$.

7.3 Numerical evidence

Take $N = 20$, $\kappa = 1000$, $\eta = 0.8/h_{\mathrm{geo}}$.

$\xi$	$L$ at step 50	$L$ at step 100
0	$2.9 \times 10^{-3}$	$1.2 \times 10^{-3}$
0.001	$4.9 \times 10^{-1}$	$2.4 \times 10^{-4}$
0.1	$77$	$6.0$
1	$188$	$175$

$\xi = 0$ converges fastest. Increasing $\xi$ progressively slows convergence by damping the Newton step. The "safety" argument that $\xi > 0$ prevents instability is redundant: the Newton metric already constrains every preconditioned eigenvalue to $h_{\mathrm{geo}}$, so $\eta$ can be set directly against $h_{\mathrm{geo}}$ without needing a separate safety mechanism. With stochastic gradients $\xi > 0$ is actively harmful: gradient noise contributes permanent damping that does not vanish as $L \to 0$.

7.4 Conclusion

Set $\xi = 0$. This eliminates the denominator, the denominator EMA, and the hyperparameters $\beta$ and $\xi$. The parameter update simplifies to

$$\theta_i \leftarrow \theta_i - \eta\,e^{s_i}\,\hat m_i.$$

8. Question 6: which embedding function?

With $\xi = 0$, the embedding is moot. If $\xi > 0$ is retained as a transient safety net before the metric converges, the choice of embedding function $f(L)$ matters.

8.1 Non-interference condition

For the Newton preconditioning to function at convergence, the damping must vanish: $r_t \to 1$ as $L \to 0$. With general embedding $\phi = (\theta, \sqrt{\xi}\,f(L))$, the optimizer descends $L$ (not $f(L)$) under the pullback metric, so the gradient in the embedded coordinates is $\sqrt{\xi}\,f'(L)\,\nabla L$ and the Sherman-Morrison denominator becomes

$$r_t = \frac{f'(L)}{1 + \xi\,f'(L)^2 \cdot 2\,h_{\mathrm{geo}}\,L},$$

where the $f'(L)$ in the numerator is the multiplicative correction from differentiating the embedded coordinate, and $\|\ell_{\mathrm{embed}}\|^2 = \xi f'(L)^2 \cdot \|\nabla L\|^2 = \xi f'(L)^2 \cdot 2\,h_{\mathrm{geo}}\,L$ on the equilibrium-preconditioned quadratic. Among power-law embeddings $f(L) = L^p$, the limit $\lim_{L \to 0} r_t = 1$ requires $p = 1$. More generally, $f'(0)$ finite and non-zero forces $f(L) \sim cL + \mathrm{const}$ near $L = 0$. The plain loss embedding (Thomas Harvey's original choice) is the unique non-interfering choice within this family.

8.2 Why the log embedding fails

For $f = \log L$, we have $f'(L) = 1/L$, so as $L \to 0$ the denominator approaches $r_t \to 1/(2 \xi h_{\mathrm{geo}})$. This is a permanent step-size cap that depends on $\xi$ but not on the geometry of the problem. The cap decouples the learning rate from the curvature, which is the opposite of what Newton preconditioning is trying to achieve.

8.3 Transient stability

Before the metric converges, the preconditioned eigenvalues span $[h_{\min}, h_{\max}]$. Stability requires $\eta\,r_t\,\max(e^{s_i}\,h_i) < 2$. We compared the maximum stable $\eta$ at $\kappa = 100$ across the three embedding choices.

Embedding	Maximum stable $\eta$ at $\kappa = 100$
Plain ($f = L$)	$\sim 1.0$
Log ($f = \log L$)	$\sim 1.0$
Identity ($\xi = 0$)	$\sim 0.03$

The plain embedding is self-stabilising: large gradients automatically reduce $r_t$, permitting much larger $\eta$ during the metric transient. The identity embedding ($\xi = 0$) gives no safety net, and the stability window is 30 times smaller during the transient. So the practical recommendation is $\xi = 0$ at deployment (Section 7), with the plain embedding kept as a fallback if the optimisation hits a transient instability and needs the geometric safety net.

9. The final algorithm

The composition of the answers to the six questions gives a clean algorithm with three free hyperparameters (plus momentum if used).

Inputs: learning rate $\eta$, metric smoothing rate $\beta_s = \mu\lambda$, clip bound $c$. Optional: momentum coefficient $\beta_m$.

Each step:

(1) Compute the gradient $g_t = \nabla L(\theta_t)$.

(2) Estimate the Hessian diagonal via Hutchinson: draw a Rademacher vector $v \in \{-1, +1\}^N$, compute $Hv$ via reverse-mode autodiff, set $\hat h_{i,t} = v_i\,(Hv)_i$.

(3) Update the metric: $s_{i,t+1} = (1 - \beta_s)\,s_{i,t} - \beta_s\,\log \max(\hat h_{i,t}, \varepsilon)$.

(4) Mean-center and clip: $s_i \leftarrow s_i - \bar s$, then $s_i \leftarrow \mathrm{clip}(s_i, -c, +c)$.

(5) Update parameters: $\theta_{t+1} = \theta_t - \eta\,\mathrm{diag}(e^s)\,g_t$.

9.1 Before versus after

	Original heuristic	Derived
Hyperparameters	10	3 (+1 if using momentum)
Metric stability	Diverges without clipping	Globally asymptotically stable
Preconditioning quality	Anti-Newton (correlation $+0.71$)	Newton (correlation $-1.0$)
Hessian estimation	Secant (1100% initial bias)	Hutchinson (unbiased)
Global damping	Unclear role	Eliminated
Curvature correction	Heuristic $\tanh$	Derived from minimax theorem
Convergence guarantee	None	Lyapunov + Borkar

9.2 Closest published cousins

Adahessian (Yao et al., 2021) uses diagonal Hutchinson Hessian estimation with momentum. The key difference is parameterisation: Adahessian works in linear space (EMA of $h^2$, then square root, in the style of Adam), while the algorithm above works in log space (EMA directly on $\log h$). The log-space parameterisation is what enables the robustness theorem.

Sophia (Liu et al., 2024) uses a diagonal Hessian with clipping. Sophia clips the update ratio $m/h$ element-wise and updates the Hessian every $k = 10$ steps. The derived optimiser clips the log-preconditioner $s$ (bounding the condition number of the metric) and updates the metric every step as a stable ODE.

What is genuinely new in the derived optimizer:

(a) the log-space stable dynamical system for the preconditioner update, with a formal Lyapunov proof;

(b) the separation principle, decomposing the update into independent dynamics and target subproblems with separate optimality results;

(c) the robustness theorem explaining log-space tolerance to estimation error.

10. Convergence rate analysis

The derived optimizer has two phases.

Phase 1, the metric transient: the EMA learns the Newton target. Metric error decays as $(1 - \beta_s)^t$, reaching tolerance $\delta$ in $T_{\mathrm{metric}} \approx (1/\beta_s)\,\log(\|e_0\|/\delta)$ steps, where $e_0$ is the initial metric error. During this phase the preconditioned eigenvalues are not yet equal: the optimizer sees a partially-corrected canyon.

Phase 2, the metric equilibrium: all preconditioned eigenvalues are approximately $h_{\mathrm{geo}}$, the effective condition number is approximately 1, and convergence is fast.

Phase 1 dominates the total iteration count. With a fixed learning rate $\eta$ chosen safely for the worst (uncorrected) direction, the iteration complexity is $O(\sqrt\kappa)$, the same as plain momentum without preconditioning. The metric converges to a perfect Newton preconditioner, but the learning rate is set for the worst case throughout. The "perfect Newton" state is only reached at the end of Phase 1, and Phase 2 is short.

10.1 The $\sqrt L$ embedding restores $O(\log \kappa)$ during the transient

The Question 5 / Question 6 conclusion ($\xi = 0$ at the optimum, plain embedding when $\xi > 0$) is the right asymptotic choice but discards an adaptive step-size mechanism that helps during the transient. A practical refinement: use the $\sqrt L$ embedding ($f(L) = \sqrt L$) with $\xi > 0$ during early training and let $\xi$ decay toward $0$ as the optimizer approaches the minimum. Among power-law embeddings $f(L) = L^p$, only $p = 1/2$ makes the denominator dimensionally scale-independent under $L \to c L$. The scale-independent denominator absorbs the worst-case eigenvalue during the metric transient, giving $O(\log \kappa)$ iteration complexity instead of $O(\sqrt \kappa)$; much better than $\sqrt \kappa$ for large $\kappa$, though not strictly an exponential improvement (the comparison is $\log$ vs $\sqrt{}$).

A quantitative comparison: on a problem with $\kappa = 10^5$, the $\sqrt L$ embedding is 11 times faster than the plain embedding (with $\xi$ tuned in both cases). Full derivation and numerical comparison in notebooks/convergence-rates.nb.

The complete recipe is therefore: use the $\sqrt L$ embedding with $\xi > 0$ during the metric transient (Phase 1), and the conclusions of Section 7 (which favour $\xi \to 0$ at convergence) take over in Phase 2 as the denominator naturally approaches $1$. The two regimes compose continuously rather than contradicting each other.

11. Failure on real architectures

The derivation lives in a problem class: diagonal quadratics with deterministic gradients. Real neural-network Hessians do not lie in this class. Two failure modes were verified directly.

11.1 Failure mode 1: diagonal Newton is wrong for non-diagonal $H$

The Question 2 theorem assumes diagonal $H$. The Gauss-Newton Hessian of a dense layer $y = Wx$ has the Kronecker structure $H_{\mathrm{layer}} = A \otimes B$, where $A = X^\top X / |\mathcal B|$ is the per-layer input second moment and $B = (G_o)^\top G_o / |\mathcal B|$ is the per-layer output-gradient second moment (averaged over batch size $|\mathcal B|$, with $G_o$ the matrix of output-gradients per sample). The symbol $\otimes$ denotes the Kronecker tensor product. The off-diagonal mass of such a matrix is structured and large.

For the MNIST MLP dense2 layer ($64 \times 64$, approximately 4{,}000 parameters) measured after 100 training steps, the condition numbers of the preconditioned Hessian under three different preconditioner choices are:

Preconditioner	$\kappa(MH)$
Identity (no preconditioning)	$1.79 \times 10^9$
Diagonal Newton	$1.16 \times 10^9$
Kronecker ($A \otimes B$)	$15.5$

Diagonal Newton barely improves the condition number (the diagonal preconditioner cannot break up the off-diagonal block structure that dominates the spectrum). The Kronecker preconditioner, which respects the true factor structure of the Hessian, yields a condition number of $15.5$, which is not literally "near 1" in the conditioning sense but is near-isotropic compared with $10^9$ (an improvement of approximately eight orders of magnitude). The Kronecker upgrade that resolves this is the subject of the companion entry on pullback unification.

11.2 Failure mode 2: Hutchinson breaks on real architectures

The variance of the Hutchinson estimator scales as $\sum_{j \ne i} H_{ij}^2$. Normalise this by the diagonal entry to get the diagonal-dominance ratio

$$\rho_i := \frac{\sqrt{\sum_{j \ne i} H_{ij}^2}}{H_{ii}}.$$

The estimator's variance, relative to the magnitude of $H_{ii}$, is $\rho_i^2$. The robustness theorem of Section 5.4 requires multiplicative error bounded by a value $\delta < 1$, which translates to $\rho \lesssim 0.3$ for the error to stay in the theorem's domain over the timescales the EMA averages.

We measured $\rho$ on real architectures. The measurement uses the Hutchinson samples themselves, via the coefficient of variation of $\hat h_i$ across multiple draws: $\rho_i = \mathrm{std}(\hat h_i) / |\mathrm{mean}(\hat h_i)|$. The code is in scripts/rho-measurements.py.

MNIST MLP (55K parameters, three dense layers), Optuna-tuned, 200 epochs:

Layer	$\rho$ at epoch 0	$\rho$ at epoch 199	$P(\hat h_i < 0)$
Layer 1 kernel	74	41	38-39%
Layer 2 kernel	74	39	48-50%
Layer 3 kernel	186	11	44-49%

Miniature 4-layer transformer (843K parameters), Optuna-tuned, 200 epochs:

Layer	$\rho$
Self-attention layer 0, value kernel	35{,}402
Self-attention layer 0, output kernel	33{,}718
Dense layer 3 kernel	6{,}598
Embedding layer 1	855

For MNIST MLP, $\rho$ values of $11$ to $186$ are approximately one to two orders of magnitude outside the usable-threshold of about $0.3$. For the transformer attention layers, $\rho$ values of $10^3$ to $10^4$ are three to five orders of magnitude outside. The sign-error probability (the empirical fraction of $\hat h_i < 0$) is 38 to 50% on the dense layers. Across 200 training epochs, the metric never gets within 10 nats of the Newton target on either architecture.

The variance scaling has a Kronecker amplification effect. For a Kronecker-structured Hessian $H = A \otimes B$, where the factor matrices $A$ and $B$ have their own diagonal-dominance ratios $\rho_A$ and $\rho_B$, the joint diagonal-dominance ratio for the Kronecker product satisfies $\rho^2 = \rho_A^2 + \rho_B^2 + \rho_A^2 \rho_B^2$. Even modest non-diagonality of each Kronecker factor compounds multiplicatively. This explains the catastrophic $\rho$ values on transformer attention layers, where the factor matrices are themselves moderately non-diagonal.

11.3 The definitive evidence: full-batch ablation

To isolate whether the failure is caused by gradient noise or by the target itself, MNIST was trained with the entire training set as one batch (so the gradient noise from minibatch sampling is zero) using the exact Gauss-Newton diagonal (so the Hessian estimation noise is zero). Running 1000 steps at a fixed learning rate, sweeping the clip $c$:

Clip $c$	$\kappa(M)$	Accuracy at step 1000	Fraction of components clipped
0.5	2.7	97.00%	92%
1.0	7.4	97.07%	84%
2.0	55	96.92%	69%
4.0	2{,}981	95.61%	37%
8.0	80{,}349	94.89%	0%

At $c = 8$ the clip never activates: the full Newton target is expressed. This is the worst result. Accuracy decreases monotonically as preconditioning aggressiveness increases. For reference, Adam reaches 97.17% under the same conditions and SGD reaches 94.22%.

This is the definitive evidence. Minibatch noise is not the cause: full-batch results are within 0.3% of minibatch results step-for-step. The cause is the target itself: diagonal Newton preconditioning of a non-diagonal Hessian inflates the effective condition number, as predicted by Section 11.1. The derived optimiser works best in this regime when the clip forces minimal preconditioning ($c \approx 1$, $\kappa(M) \approx 7$), at which point the optimiser is essentially SGD with a mild rescaling and matches Adam. The full reproduction code is in scripts/fullbatch-ablation.py.

11.4 Standing in the 1000-trial sweep

The curvature-aware diagonal (sgd_learn_diag_curv) — the optimizer this entry derives and cleans up — is in the project's 1000-trial MNIST MLP Optuna sweep (results/mnist_mlp/). Its standing corroborates the failure analysis: a high best-run ceiling but a median well below the baselines.

Optimizer	Best	Median	Mean	Worst
sgd_learn_diag_curv	98.39%	78.04%	60.88%	5.36%
Adam	98.09%	94.73%	72.69%	9.80%
SGD	98.04%	78.19%	58.43%	2.57%
Muon	98.39%	93.73%	70.53%	9.29%

The ceiling ties Muon for the best single run, but the median ($78.04\%$) sits ${\approx}16$ points below Adam — the same high-ceiling / low-robustness signature the full-batch ablation (§11.3) traces to diagonal misspecification against Kronecker-structured Hessians.

12. What survives

The derived optimiser is mathematically clean. Every component has a theorem behind it. Ten hyperparameters reduce to three. Lyapunov stability is global. The robustness theorem is striking. The Lyapunov-stable preconditioner update with the natural-gradient correction is a foundational contribution: no prior optimizer frames the preconditioner update as a provably stable ODE.

The Question 2 minimax theorem is tautological in the formal sense, as discussed in the companion entry on the master minimax framework. An instance of the framework is tautological when the preconditioner class can exactly invert any element of the problem class: in that case the "minimax-optimal preconditioner" is just the inverse, and the theorem is an algorithmic re-derivation of a known optimizer (Newton, in this case). The same is true for the Kronecker-level upgrade, which reproduces KFAC (Martens and Grosse, 2015).

On real neural networks, the diagonal preconditioner class is structurally misspecified for the Kronecker structure of layer Hessians (Section 11.1), and the Hutchinson estimator's variance scales quadratically with each Kronecker factor's non-diagonality (Section 11.2). Both failure modes are mechanistically explained, not statistically noisy.

What survives as contribution: the Lyapunov / natural-gradient stable update for the preconditioner; the separation principle for dynamics and target; the robustness theorem explaining why log-space targets tolerate enormous estimation error within their domain; the unified pullback-metric view connecting SGD, Newton, Gauss-Newton, Fisher, and KFAC; and precise mechanistic explanations for why naive diagonal Newton and Hutchinson curvature estimation do not transfer to neural networks.

What does not survive: a "derived state-of-the-art optimizer from differential geometry" claim. On MNIST, the derived optimizer loses to SGD. Hutchinson is not a viable curvature estimator for real architectures. Diagonal is not the right structural level for neural-network preconditioning.

The framework's contribution at the diagonal level is foundational understanding rather than a new dominant optimizer. The natural next move, minimax at the Kronecker level, produces KFAC, which is itself competitive but not Adam-beating on benchmarks.

Appendix A: Lyapunov stability of the corrected metric update

This appendix gives the full proof of global asymptotic stability for the metric update $\dot s_i = \mu[\xi g_i^2 - \lambda s_i]$.

Assume the gradient components $g_i$ are constants (this is the appropriate setting for fixed-point analysis of the $s$ subsystem; the joint case is in Appendix D). The unique equilibrium is $s_i^* = \xi g_i^2 / \lambda$. Define the Lyapunov candidate

$$V(s) = \tfrac{1}{2} \|s - s^*\|^2 = \tfrac{1}{2} \sum_i (s_i - s_i^*)^2.$$

This is non-negative, equal to zero only at $s = s^*$, and radially unbounded. Differentiate along the dynamics:

$$\dot V = \sum_i (s_i - s_i^*)\,\dot s_i = \sum_i (s_i - s_i^*)\,\mu\,(\xi g_i^2 - \lambda s_i).$$

Substitute $\xi g_i^2 = \lambda s_i^*$ from the equilibrium definition:

$$\dot V = \sum_i (s_i - s_i^*)\,\mu\,(\lambda s_i^* - \lambda s_i) = -\mu \lambda \sum_i (s_i - s_i^*)^2 = -2 \mu \lambda V.$$

So $V(t) = V(0)\,e^{-2 \mu \lambda t}$, decaying exponentially at rate $2 \mu \lambda$. The decay is global: no linearisation was used and the result holds for all initial conditions. This is global asymptotic stability.

In discrete time, with stepsize $\Delta t$, the update is $s_{t+1} = (1 - \mu\lambda \Delta t) s_t + \mu \xi g^2 \Delta t$. Stability requires $|1 - \mu\lambda \Delta t| < 1$, which holds for $\Delta t < 2/(\mu\lambda)$. Standard implementations use $\Delta t = 1$ and absorb $\mu\lambda$ into a single rate $\beta_s = \mu\lambda$ with $\beta_s \in (0, 1)$.

Appendix B: Robustness theorem

Let $H$ be a diagonal $N \times N$ matrix with positive diagonal entries $H_{ii}$, and let $\hat H_{ii}$ be estimates with multiplicative error bounded by $\delta \in [0, 1)$:

$$\left|\frac{\hat H_{ii}}{H_{ii}} - 1\right| \le \delta \quad \text{for every } i.$$

With Newton-style preconditioner $M = \mathrm{diag}(\hat H_{ii}^{-1})$ multiplied by a normalising scalar $\alpha = (\prod_i \hat H_{ii})^{1/N}$ (which corresponds to mean-centering in $s$-space), the preconditioned matrix $MH$ has diagonal entries

$$(MH)_{ii} = \alpha\,\frac{H_{ii}}{\hat H_{ii}}.$$

By the multiplicative error bound, $H_{ii}/\hat H_{ii} \in [1/(1+\delta), 1/(1-\delta)]$. Therefore

$$\frac{\max_i (MH)_{ii}}{\min_i (MH)_{ii}} \le \frac{1/(1-\delta)}{1/(1+\delta)} = \frac{1+\delta}{1-\delta}.$$

This is the bound on the preconditioned condition number. Note that the bound is independent of $\kappa(H)$ and independent of $N$.

Numerical illustration: at $\delta = 0.5$, the bound gives $\kappa_{\mathrm{prec}} \le 3$; at $\delta = 0.9$, the bound gives $\kappa_{\mathrm{prec}} \le 19$. Either way, this is a dramatic reduction from $\kappa(H)$ even when the estimation error is large.

The bound is tight in the worst case: take $\hat H_{ii}$ to alternate between $H_{ii}(1+\delta)$ and $H_{ii}(1-\delta)$, and the bound is achieved with equality.

The reason the theorem is so strong is the log-space parameterisation. A multiplicative error in $\hat H_{ii}$ is an additive error in $\log \hat H_{ii}$, hence in $s_i$. The target spectrum spans $\log(\kappa(H))$ nats, so an additive error of $\log(1+\delta) \approx \delta$ nats is small relative to the target span. The condition-number bound is the exponential of twice the additive error in $s$, which is exactly $(1+\delta)/(1-\delta)$.

The theorem does not apply if the estimate has the wrong sign (multiplicative error greater than 1 in magnitude). On real architectures the wrong-sign rate is 40 to 50% (Section 11.2), placing the estimator outside the theorem's domain.

Appendix C: Hutchinson estimator details

The Hutchinson estimator for the diagonal of a square matrix $H$ uses a random vector $v$ to project the matrix into a scalar that, in expectation, recovers a chosen diagonal entry.

Algorithm. Draw $v_i \in \{-1, +1\}$ independently with equal probability (Rademacher distribution). Compute $w = Hv$ via one application of $H$. The Rademacher choice can be replaced by any zero-mean unit-variance distribution; the Rademacher distribution is variance-optimal within this class (see below). Set $\hat h_i = v_i\,w_i$ for each $i$.

Unbiasedness. $\mathbb{E}[v_i v_j] = \delta_{ij}$ because the components are independent. Therefore

$$\mathbb{E}[\hat h_i] = \mathbb{E}\Bigl[v_i \sum_j H_{ij}\,v_j\Bigr] = \sum_j H_{ij}\,\mathbb{E}[v_i v_j] = \sum_j H_{ij}\,\delta_{ij} = H_{ii}.$$

Variance.

$$\mathrm{Var}(\hat h_i) = \mathbb{E}[\hat h_i^2] - H_{ii}^2 = \mathbb{E}\Bigl[v_i^2\,\Bigl(\sum_j H_{ij}\,v_j\Bigr)^2\Bigr] - H_{ii}^2.$$

Expanding the square and using $v_i^2 = 1$ identically (the special feature of the Rademacher distribution),

$$\mathrm{Var}(\hat h_i) = \sum_{j, k} H_{ij}\,H_{ik}\,\mathbb{E}[v_j\,v_k] - H_{ii}^2 = \sum_j H_{ij}^2 - H_{ii}^2 = \sum_{j \ne i} H_{ij}^2.$$

The variance is the off-diagonal squared row sum of $H$, normalised by $H_{ii}^2$ it is exactly $\rho_i^2$.

Why Rademacher specifically. Among zero-mean unit-variance distributions, the Rademacher distribution achieves the lower bound on $\mathbb{E}[v_i^4]$ (which is 1 because $v_i^2 = 1$ identically). The variance calculation involves $\mathbb{E}[v_i^2 v_j^2]$ terms; these are minimised by Rademacher. A Gaussian $v$ would give a variance constant of 2 times the Rademacher value, which is the standard trade-off when choosing between the two.

Computational cost. One $Hv$ product is one extra reverse-mode autodiff pass per training step, approximately doubling the per-step cost relative to a vanilla gradient computation. For amortisation across multiple samples per step (the multi-Hutchinson estimator), the cost scales linearly in the number of samples.

Appendix D: Joint $(\theta, s)$ Lyapunov function

The joint Lyapunov function for the coupled $(\theta, s)$ system on a diagonal quadratic loss is

$$V(\theta, s) = L(\theta) + \tfrac{\alpha}{2}\,\|s - s^*\|^2$$

for any $\alpha > 0$. We show that $\dot V \le 0$ along the joint dynamics, with equality only at the joint equilibrium.

The joint dynamics. The parameter update on a quadratic with diagonal $H = \mathrm{diag}(h)$ is $\dot \theta_i = -\eta\,r\,e^{s_i}\,g_i = -\eta\,r\,e^{s_i}\,h_i\,\theta_i$. The metric update with the Newton-targeted driving force is $\dot s_i = \mu\,\lambda\,(s_i^* - s_i)$, where $s_i^* = -\log h_i + \overline{\log h}$ is constant on a fixed quadratic (it does not depend on $\theta$). This is the form used here. The earlier gradient-magnitude rule $\dot s_i = \mu[\xi g_i^2 - \lambda s_i]$ has a $\theta$-dependent target via $g_i = h_i \theta_i$, and the Lyapunov argument below would need modification to account for that coupling; on a quadratic, the Newton-targeted form has the convenient property of a fixed target, which is what makes the joint analysis clean.

The loss is $L(\theta) = \tfrac{1}{2} \sum_i h_i\,\theta_i^2$, so $\dot L = \sum_i h_i\,\theta_i\,\dot \theta_i = -\eta\,r \sum_i e^{s_i}\,h_i^2\,\theta_i^2$. This term is non-positive: $\eta, r > 0$, $e^{s_i} > 0$, $h_i^2 > 0$, $\theta_i^2 \ge 0$.

The metric error term has $\tfrac{1}{2}\,d/dt(\|s - s^*\|^2) = \sum_i (s_i - s_i^*)\,\dot s_i = -\mu\lambda \|s - s^*\|^2$ by direct substitution of $\dot s_i = \mu\lambda(s_i^* - s_i)$. This is also non-positive.

Therefore

$$\dot V = \dot L + \tfrac{\alpha}{2}\,\frac{d}{dt}\|s - s^*\|^2 = -\eta\,r \sum_i e^{s_i}\,h_i^2\,\theta_i^2 \;-\; \alpha\,\mu\,\lambda\,\|s - s^*\|^2.$$

Both terms are non-positive, and the sum is zero only when $\theta = 0$ and $s = s^*$. This is global asymptotic stability of the joint system.

The parameter $\alpha > 0$ is free in the sense that any positive choice gives a valid Lyapunov function. In particular, no relationship between $\mu$ and $\eta$ is required: the joint dynamics are stable for any ratio $\mu/\eta > 0$. This was a key concern at the start of the program (whether the slow-fast timescale ratio mattered), and the Lyapunov calculation shows that on quadratics it does not.

Beyond quadratics: at a non-degenerate local minimum, the joint Jacobian is block-triangular because $\ell = 0$ at the minimum kills the cross-coupling. Local asymptotic stability follows from each diagonal block being a contraction. With diminishing step sizes and additional assumptions (martingale-difference noise, Lipschitz dynamics, boundedness of the iterates, and the standard timescale-ratio condition $\eta_t/\mu_t \to 0$), Borkar's two-timescale stochastic-approximation theorem (Borkar, 1997) gives almost-sure convergence to the stable equilibrium set; with constant step sizes, the iterates converge to an $O(\eta + \mu)$ neighbourhood. The Borkar hypotheses are non-trivial to verify on real architectures and we do not claim to do so here; the role of the theorem is to suggest that the joint stability of the quadratic analysis extends to the stochastic case under conventional regularity, not to provide an end-to-end proof of convergence for deep network training.

Appendix E: Full-batch ablation methodology

The full-batch ablation in Section 11.3 is constructed to isolate the contribution of each potential cause of the diagonal-Newton failure on MNIST. We list the controls explicitly.

The training set is the full 60{,}000-sample MNIST training set, used as a single batch for every gradient evaluation. This eliminates minibatch sampling noise from the gradient.

The Hessian diagonal is computed exactly via the Gauss-Newton approximation, $H_{\mathrm{GN}} = J^\top J$ for a regression layer and the appropriate generalisation for softmax-cross-entropy. The diagonal entries are extracted analytically rather than estimated via Hutchinson. This eliminates Hessian-estimation noise.

The model is a three-layer MLP with ReLU activations, batch normalisation off, no dropout, fixed random seed. The model has approximately 55K parameters. The clip parameter $c$ is swept across $\{0.5, 1.0, 2.0, 4.0, 8.0\}$. All other hyperparameters are held fixed at their values from the best Optuna run for $c = 1$.

The accuracy reported is on the held-out 10{,}000-sample test set after 1000 gradient steps. Adam and SGD baselines use the same loss function and step count with their own tuned hyperparameters.

The result is that accuracy decreases monotonically with $c$ when both gradient noise and Hessian noise are removed. Since the only thing varying is the aggressiveness of the diagonal Newton preconditioning, the cause of the failure is the preconditioning itself, not the estimation. This is the experimental complement to the structural argument of Section 11.1 (Kronecker structure makes diagonal preconditioning misspecified).

Appendix F: Attached materials

Notebooks (notebooks/):

• notebooks/newton-minimax-comparison.nb: symbolic verification of the Newton minimax theorem and numerical comparison on the 10-dimensional quadratic. Generates the Section 4.4 table.
• notebooks/lyapunov-stability.nb: symbolic verification of the Appendix A Lyapunov calculation and numerical integration of the metric ODE from multiple initial conditions, confirming exponential decay.
• notebooks/hutchinson-robustness.nb: numerical verification of the robustness theorem for sweeps over $\delta$, plus the Section 5.3 estimator comparison table.
• notebooks/joint-lyapunov-verification.nb: numerical integration of the coupled $(\theta, s)$ dynamics from 2000 random initialisations on quadratics of dimension up to 50; confirms monotonic decrease of $V(\theta, s)$.
• notebooks/convergence-rates.nb: phase analysis of the $O(\sqrt\kappa)$ scaling, derivation of the $\sqrt L$ embedding $O(\log\kappa)$ result, and the 11-times speedup comparison on $\kappa = 100{,}000$.

Scripts (scripts/):

• scripts/rho-measurements.py: measures the diagonal-dominance ratio $\rho$ from Hutchinson samples on a given model and dataset, producing the tables in Section 11.2.
• scripts/fullbatch-ablation.py: reproduces the Section 11.3 MNIST full-batch sweep, including the clip-sweep table.
• scripts/make_figures.py: reproduces every figure shown in this entry (PDF and PNG written to figures/).

References

Borkar, V.S. (1997). Stochastic approximation with two time scales. Systems and Control Letters, 29(5), 291-294.

Hutchinson, M.F. (1989). A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines. Communications in Statistics: Simulation and Computation, 18(3), 1059-1076.

Liu, H., Tian, X., Chen, X., Tao, R., Liu, Y. (2024). Sophia: A scalable stochastic second-order optimizer for language model pre-training. International Conference on Learning Representations.

Martens, J., Grosse, R. (2015). Optimizing neural networks with Kronecker-factored approximate curvature. International Conference on Machine Learning, 2408-2417.

Yao, Z., Gholami, A., Shen, S., Mustafa, M., Keutzer, K., Mahoney, M.W. (2021). AdaHessian: An adaptive second order optimizer for machine learning. AAAI Conference on Artificial Intelligence, 35(12), 10665-10673.