---
title: Differentiable Levenberg-Marquardt Optimization
url: https://www.emergentmind.com/topics/fully-differentiable-levenberg-marquardt-optimization
type: topic
---

# Differentiable Levenberg-Marquardt Optimization

Searching arXiv for the specified papers and closely related work on differentiable Levenberg–Marquardt optimization.
Search query: title:"Intraoperative 2D/3D Registration via Spherical Similarity Learning and Inference-Time Differentiable Levenberg-Marquardt Optimization"
Fully differentiable Levenberg–Marquardt optimization denotes an implementation of the Levenberg–Marquardt (LM) method in which the full update path—residual construction, Jacobian extraction, curvature formation, linear solve, damping manipulation, and parameter update—remains inside an automatic-differentiation computation graph. In the cited literature, this formulation appears both as an inference-time pose-refinement mechanism for intraoperative 2D/3D registration and as a differentiable optimization layer for deep neural network training based on subsampled Gauss–Newton and natural-gradient ideas [2509.06890], [1906.02353]. Across these settings, the shared objective is a nonlinear least-squares problem of the form $C(x)=\tfrac12\|r(x)\|^2$ or $J(\theta)=\tfrac12\sum_i \|r_i(\theta)\|^2$, with LM supplying a damped Gauss–Newton step that interpolates between Gauss–Newton and gradient descent.

## 1. Problem formulation and residual structure

The central optimization problem is nonlinear least squares. In the inference-time registration formulation, a parametric pose vector $x\in\mathbb R^6$ is refined by minimizing

$$
C(x)=\frac12\|r(x)\|^2,
$$

where $r(x)\in\mathbb R^m$ is a residual vector encoding a per-pixel or per-feature discrepancy derived from the learned spherical similarity [2509.06890]. In practice, if $\phi_m(x)$ and $\phi_f$ are two $H\times W\times D$ feature fields, the residual at pixel $(i,j)$ is defined as

$$
r_{ij}(x)=\Phi_m[i,j]^T\cdot \Phi_f[i,j]-1,
$$

and stacking all $H\cdot W$ residuals yields an $m=H\cdot W$ dimensional residual vector.

In the training-oriented formulation for neural networks, the residuals are indexed by data point. Given a dataset $\{(x_i,y_i)\}_{i=1}^N$ and a model $\hat y_i(\theta)=F(x_i;\theta)$, one defines

$$
r_i(\theta)=F(x_i;\theta)-y_i\in\mathbb R^m,
$$

and the least-squares cost

$$
J(\theta)=\frac12\sum_{i=1}^N \|r_i(\theta)\|^2
=\frac12\sum_{i=1}^N r_i(\theta)^T r_i(\theta).
$$

The exact gradient is

$$
\nabla J(\theta)=\sum_{i=1}^N \frac{\partial r_i(\theta)}{\partial \theta}^T r_i(\theta)
=\sum_{i=1}^N J_i(\theta)^T r_i(\theta),
$$

where $J_i(\theta)=\frac{\partial r_i(\theta)}{\partial \theta}\in\mathbb R^{m\times n}$ [1906.02353].

These two formulations differ in application domain and dimensional regime, but they share the same structural decomposition: residuals, Jacobians, approximate second-order curvature, and a damped linear system. A plausible implication is that “fully differentiable LM” is best understood as a design pattern for embedding classical nonlinear least-squares updates into learned systems rather than as a single domain-specific algorithm.

## 2. LM update as damped Gauss–Newton

At iterate $x_k$, the registration formulation linearizes the residual as

$$
r(x_k+\delta)\approx r(x_k)+J\delta,
$$

with $J=\partial r/\partial x|_{x_k}\in\mathbb R^{m\times n}$ and $n=6$. The resulting quadratic model is

$$
\frac12\|r+J\delta\|^2+\frac{\lambda}{2}\delta^T D\delta,
$$

where $D$ is a positive-definite damping matrix, often chosen as $\operatorname{diag}(J^T J)$ or the identity. Setting the derivative with respect to $\delta$ to zero yields the normal equations

$$
(J^T J+\lambda D)\delta=-J^T r,
$$

and hence the LM step

$$
\delta=-(J^T J+\lambda D)^{-1}J^T r,
$$

followed by $x_{k+1}=x_k+\delta$ [2509.06890].

The training-oriented formulation presents the same idea in parameter space. The Gauss–Newton matrix is the positive semi-definite approximation

$$
G_N(\theta)\approx \sum_{i=1}^N J_i(\theta)^T J_i(\theta),
$$

and the damped LM update solves

$$
(G_N(\theta)+\lambda I)\Delta\theta=-\nabla J(\theta),
$$

so that

$$
\Delta\theta=-(G_N(\theta)+\lambda I)^{-1}\nabla J(\theta).
$$

In both cases, the damping parameter $\lambda$ is the mechanism that makes LM interpolate between Gauss–Newton for small $\lambda$ and gradient descent for large $\lambda$ [1906.02353].

This algebra clarifies the role of differentiability. The requirement is not merely that the objective be differentiable; rather, the linearized solve itself must be implemented so that gradients can propagate through the update map. That distinction is central to the cited work.

## 3. End-to-end differentiability in automatic-differentiation frameworks

The inference-time registration implementation specifies that all key operations are realized with auto-differentiable tensor operations: residual construction by forward propagation through the CNN+Transformer, spherical exponential map, and inner-product; Jacobian computation by back-propagating the residuals with respect to the six pose parameters; formation of $A=J^T J+\lambda D$ and $g=J^T r$ by differentiable matrix multiplications; solution of $A\delta=-g$ through a differentiable linear solve such as Cholesky factorization or triangular solve; and pose update by differentiable addition, with retraction or group composition if one works directly on the manifold [2509.06890]. Because every tensor operation is tracked by PyTorch or TensorFlow, the full LM loop is end-to-end differentiable.

The training-oriented treatment makes the same point in a higher-dimensional regime, while emphasizing that one never forms $J$ explicitly in frameworks like PyTorch or TensorFlow. Instead, two primitives suffice: for $v\in\mathbb R^n$, one computes $Jv$ by forward-mode AD or “perturb-and-clone” tricks; for $u\in\mathbb R^m$, one computes $J^T u$ by reverse-mode AD. These primitives are then used to accumulate the gradient, form $(J J^T)w$ products, solve for $\Delta\theta$ by Sherman–Morrison–Woodbury or a small linear solve, and differentiate through the full LM step as a differentiable layer [1906.02353].

A common misconception is that a differentiable second-order method necessarily requires explicit Jacobian materialization. The training-oriented formulation directly contradicts that view: explicit storage of the full Jacobian is not required, and the LM step can still be differentiated through linear-algebra operations already supported by deep-learning frameworks.

## 4. Inference-time differentiable LM for 2D/3D registration

In the intraoperative 2D/3D registration setting, the optimization problem is embedded within a broader geometric learning pipeline. The cited method addresses the limitation that existing Euclidean approximations distort manifold structure and slow convergence. To mitigate this, it explores similarity learning in non-Euclidean spherical feature spaces, extracts feature embeddings using a CNN-Transformer encoder, projects them into spherical space, and approximates their geodesic distances with Riemannian distances in the bi-invariant SO(4) space. The stated purpose is a more expressive and geometrically consistent deep similarity metric that enhances the ability to distinguish subtle pose differences [2509.06890].

Within that pipeline, inference replaces gradient descent with fully differentiable Levenberg–Marquardt optimization. The high-level loop is: compute the residual vector $r(x)$, evaluate the cost $C=\tfrac12\|r\|^2$, compute the Jacobian $J=\partial r/\partial x$, form $A=J^T J+\lambda D$ and $g=J^T r$, solve $A\delta=-g$, evaluate a trial pose $x_{\mathrm{trial}}=x+\delta$, recompute the residual and trial cost, accept or reject the step based on cost decrease, adapt $\lambda$, and stop when $\|\delta\|<\epsilon$ or when the maximum number of iterations is reached [2509.06890].

The damping rule follows the standard LM heuristic. If $C(x_k+\delta)<C(x_k)$, the step is accepted and $\lambda\leftarrow \lambda/\nu$; otherwise, the step is rejected and $\lambda\leftarrow \lambda\cdot \nu$. This is the mechanism by which the method moves toward Gauss–Newton on successful steps and damps more aggressively on unsuccessful ones.

The reported empirical behavior is specific: the differentiable LM converges in fewer than 10 iterations to a sub-millimeter pose error, whereas SGD or Adam needed tens of iterations and still achieved worse final loss. The overall runtime per registration is reported as $\approx 6$ s with LM versus $>10$ s for a same-step-size gradient-based solver. The paper further states that its patient-specific and patient-agnostic experiments show state-of-the-art sub-millimeter accuracy with real-time-compatible runtimes [2509.06890].

## 5. Subsampled Gauss–Newton, stochastic variants, and low-rank inversion

In deep-network training, the dimensional imbalance is different: the number of parameters $n$ can be enormous, and full curvature formation is impractical. The cited work therefore introduces practical LM variants of Gauss–Newton and natural-gradient methods based on subsampling [1906.02353]. Two nested mini-batches are used at iteration $t$: $S_1\subset\{1,\dots,N\}$ of size $N_1$ for the gradient estimate and $S_2\subset S_1$ of size $N_2$ for the curvature estimate. The resulting approximations are

$$
g_{\rm sub}(\theta)=\frac1{|S_1|}\sum_{i\in S_1} J_i(\theta)^T r_i(\theta),
$$

and

$$
G_{N,\rm sub}(\theta)=\frac1{|S_2|}\sum_{i\in S_2} J_i(\theta)^T J_i(\theta).
$$

The fully or semi-stochastic LM step remains

$$
\Delta\theta=-\bigl(G_{N,\rm sub}(\theta)+\lambda I\bigr)^{-1} g_{\rm sub}(\theta).
$$

Damping adaptation is expressed through the ratio

$$
\rho=
\frac{J(\theta)-J(\theta+\Delta\theta)}
{\tfrac12\,\Delta\theta^T(G_{N,\rm sub})\Delta\theta+g_{\rm sub}^T\Delta\theta}.
$$

If $\rho$ is small, $\lambda$ is increased, described as “shrink trust-region”; if $\rho$ is large, $\lambda$ is decreased, described as “enlarge trust-region” [1906.02353].

When $n\gg m|S_2|$, the matrix $G_{N,\rm sub}=\tfrac1{|S_2|}J^T J$ is low-rank, and the Sherman–Morrison–Woodbury formula is used to avoid inverting an $n\times n$ matrix. Writing $J\in\mathbb R^{(|S_2|m)\times n}$, one obtains

$$
(G_{N,\rm sub}+\lambda I)^{-1}
=
\frac1\lambda
\left(
I-\frac1{|S_2|}J^T\left(\lambda I+\frac1{|S_2|}JJ^T\right)^{-1}J
\right),
$$

and therefore

$$
\Delta\theta=
-\frac1\lambda
\left[
I-\frac1{|S_2|}J^T D^{-1}J
\right]g_{\rm sub},
\qquad
D=\lambda I+\frac1{|S_2|}JJ^T.
$$

This formulation is not merely an algebraic convenience. It is the mechanism by which LM becomes computationally viable in parameter regimes where direct second-order methods would otherwise be prohibitive.

## 6. Convergence, complexity, and conceptual interpretation

The two cited settings provide complementary characterizations of performance. In the registration setting, the per-iteration cost of a forward pass and Jacobian extraction is $O(m\cdot n)$ with $m\approx H\cdot W$ and $n=6$; forming $A=J^T J$ costs $O(n^2\cdot m)$, but with $n=6$ this is negligible; and solving the $n\times n$ system is $O(n^3)=O(6^3)=\text{constant}$. By contrast, a naive gradient-descent update requires computing $\nabla_x C=J^T r$ at cost $O(m\cdot n)$ and a tunable step size, but empirically needs many more iterations to reach the same residual norm [2509.06890].

In the training setting, the complexity accounting is broader. With $n$ parameters, output dimension $m$, and mini-batch sizes $N_1,N_2$, the cited cost breakdown is: forward plus backprop for the gradient, $O(N_1\cdot n)$; building low-rank Jacobian blocks, $O(N_2\cdot n)$ calls to $Jv$ or $J^T v$; forming $D\in\mathbb R^{(mN_2)\times(mN_2)}$, $O(m^2N_2^2)$; inverting $D$, $O((mN_2)^3)$; and the SMW update, $O(n\cdot mN_2+(mN_2)^2)$. The total is summarized as

$$
O\!\bigl(N_1n + N_2n + (mN_2)^2 + (mN_2)^3 + n\,mN_2\bigr).
$$

Memory is dominated by activation storage $O(N_1 n)$, low-rank Jacobian factors, and the matrix $D$ of size $(mN_2)^2$ [1906.02353].

The convergence statements are likewise domain-specific. For the semi-stochastic training variant, under the assumptions that the mini-batch Gauss–Newton blocks $B_t$ are uniformly bounded, $\|B_t\|\le \beta$, and the full Hessian is Lipschitz-bounded, $\|\nabla^2 J(\theta)\|\le L$, a semi-stochastic LM method with full gradient, subsampled $G_N$, and step acceptance only when $\rho\ge \eta$ converges to a stationary point, namely $\|\nabla J(\theta_t)\|\to 0$. The proof is stated to follow classical trust-region arguments with LM–TR equivalence and control of the $\rho$-ratio [1906.02353].

Taken together, these results suggest a bifurcated interpretation of fully differentiable LM. In small-state inference problems such as se(3) pose refinement, its value lies in fast local convergence with differentiable coupling to learned residuals. In large-parameter training problems, its value lies in making second-order structure accessible through subsampling, low-rank algebra, and AD-compatible matrix-free operations.

## 7. Scope, misconceptions, and research significance

The cited literature places fully differentiable LM at the intersection of classical nonlinear least squares, manifold-aware geometric learning, and differentiable programming. In the registration case, the method is embedded in a system that explicitly seeks to better capture and fit complex manifold structure by moving from Euclidean approximations to spherical feature spaces and bi-invariant SO(4) geometry [2509.06890]. In the training case, the method is embedded in large-scale neural optimization, where practical LM variants are derived from subsampled Gauss–Newton and Fisher-information ideas and realized via SMW and automatic differentiation [1906.02353].

Several clarifications follow directly from these formulations. First, “fully differentiable” does not mean that the optimizer ceases to be iterative; the registration implementation explicitly runs an LM loop with residual recomputation, trial updates, acceptance or rejection, and convergence checks. Second, “second-order” does not imply use of the exact Hessian; both papers rely on Gauss–Newton structure, with damping added for robustness. Third, differentiability does not eliminate the need for damping policy; in both settings, adaptation of $\lambda$ remains central to the algorithmic behavior.

A plausible implication is that the main conceptual contribution of fully differentiable LM is not a new objective class but a new systems-level integration: classical LM updates are recast as differentiable computational modules that can be inserted into learned pipelines without severing gradient flow. In the available applications, this integration is used either to accelerate inference-time registration or to construct scalable LM-style updates for deep networks.

Source: https://www.emergentmind.com/topics/fully-differentiable-levenberg-marquardt-optimization