---
title: Derivative-Free Levenberg–Marquardt
url: https://www.emergentmind.com/topics/derivative-free-levenberg-marquardt-algorithm
type: topic
---

# Derivative-Free Levenberg–Marquardt

Searching arXiv for the cited papers to ground the article in the current literature.
Derivative-free Levenberg–Marquardt algorithms are variants of the Levenberg–Marquardt method for nonlinear least-squares and inverse problems that avoid direct use of analytic Jacobians or classical derivatives. In the literature represented here, the term covers several distinct constructions: replacement of unavailable Fréchet derivatives by Bouligand subderivatives for non-smooth forward maps [1902.10596]; stochastic Jacobian approximation by orthogonal spherical smoothing for nonlinear least squares [2407.12542]; Jacobian-free schemes based on Broyden rank-one updates with Armijo line search for inverse problems [2203.04138]; derivative-free enhancements using geodesic acceleration, uphill acceptance, and Broyden updates [1201.5885]; and compressed-sensing-based sparse Jacobian reconstruction for sparse nonlinear least-squares problems [2507.06772]. Across these variants, the common objective is to preserve the regularized Gauss–Newton structure of Levenberg–Marquardt while relaxing the requirement of explicit derivative information.

## 1. Concept and scope

The classical Levenberg–Marquardt method solves nonlinear least-squares problems by forming a damped normal equation based on the Jacobian. For residual map $r(x)\in\mathbb{R}^m$ and objective
\[
f(x)=\frac12\|r(x)\|^2,
\]
the gradient is $\nabla f(x)=J(x)^\top r(x)$, with $J(x)\in\mathbb{R}^{m\times n}$ the residual Jacobian [2407.12542]. In inverse problems, an analogous role is played by the derivative of a forward operator $F:D(F)\subset U\to Y$, often between Hilbert spaces [1902.10596].

The derivative-free label does not denote a single algorithmic template. In one usage, it means “derivative-free relative to classical derivatives”: the Bouligand–Levenberg–Marquardt method replaces non-existing Fréchet derivatives by computable Bouligand elements tailored to a non-smooth PDE operator [1902.10596]. In another usage, it means Jacobian approximation from function values alone, either through orthogonal spherical smoothing [2407.12542], Broyden updates [1201.5885, 2203.04138], or sparse recovery from interpolation constraints [2507.06772]. This suggests that the phrase is best interpreted structurally rather than literally: the method remains derivative-based at the model level, but the derivative information is obtained indirectly or generalized.

A further unifying feature is damping. Even when the Jacobian is exact, Levenberg–Marquardt stabilizes the Gauss–Newton system by solving a regularized linear system. The same pattern persists in derivative-free variants, where the approximate or generalized Jacobian enters a system of the form
\[
(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)
\]
or its problem-specific analogue [2407.12542, 2507.06772]. In ill-posed inverse problems, this damping is explicitly interpreted as iterative regularization [1902.10596].

## 2. Generalized-derivative formulations for non-smooth inverse problems

A prominent nonclassical formulation appears in the Bouligand–Levenberg–Marquardt iteration for ill-posed inverse problems with non-smooth forward operators [1902.10596]. The setting uses real Hilbert spaces $U$ and $Y$, a forward operator $F:D(F)\subset U\to Y$, exact data $y^\dagger=F(u^\dagger)$, and noisy data $y^\delta$ satisfying
\[
\|y^\delta-y^\dagger\|_Y\le \delta.
\]
The inverse problem is to recover $u$ from $F(u)=y^\delta$ [1902.10596].

The motivating PDE is the non-smooth semilinear elliptic problem
\[
-\Delta y+y^+=u \quad \text{in }\Omega,\qquad y=0\quad \text{on }\partial\Omega,
\]
with $U=Y=L^2(\Omega)$ and $y^+(x):=\max(y(x),0)$, where $\Omega\subset\mathbb{R}^d$ for $d=2,3$ is bounded with Lipschitz boundary [1902.10596]. In this case the forward map fails to be Gâteaux differentiable at $u$ if the zero level set $\{x:F(u)(x)=0\}$ has positive measure; precisely, Gâteaux differentiability holds at $u$ if and only if $|\{F(u)=0\}|=0$ [1902.10596].

The substitute for the missing derivative is the Bouligand subdifferential
\[
\partial_B F(u)
:= \left\{ G \in \mathcal{L}(U,Y) \;\middle|\;
\exists\, \{u_n\}\subset D\ \text{with}\ u_n \to u \text{ in }U,\ F'(u_n)\ \text{exists, and}\ \forall h\in U,\ F'(u_n)h \to Gh \text{ in }Y\right\},
\]
where $D$ is the set of Gâteaux differentiability points [1902.10596]. For the PDE application, an explicit Bouligand element $G_u\in\partial_BF(u)$ is defined by solving the linearized PDE
\[
-\Delta \zeta + 1_{\{y_u>0\}}\,\zeta = h \quad \text{in }\Omega,\qquad \zeta=0 \quad \text{on }\partial\Omega,
\]
with $y_u=F(u)$ [1902.10596]. The active-set indicator $1_{\{y_u>0\}}$ incorporates the non-smoothness directly.

The corresponding Levenberg–Marquardt step replaces the classical derivative by $G_k\in\partial_BF(u_k)$ and solves
\[
\big(G_k^* G_k + \alpha_k I\big)\, s_k = -\, G_k^*\big(F(u_k) - y^\delta\big),\qquad
u_{k+1}=u_k+s_k,
\]
with a geometrically decreasing a priori parameter rule
\[
\alpha_k=\alpha_0 r^k,\qquad \alpha_0>0,\quad r\in(0,1)
\]
and Morozov discrepancy stopping
\[
\|F(u_{k_*})-y^\delta\|_Y\le \tau\delta,\qquad \tau>1,
\]
where $k_*$ is minimal [1902.10596].

The analysis uses a generalized tangential cone condition,
\[
\|F(\hat u)-F(u)-G_u(\hat u-u)\|_Y \le \eta(\rho)\,\|F(\hat u)-F(u)\|_Y,
\]
with $\eta(\rho_0)<1$, together with transfer operators $Q(u_1,u_2)$ satisfying
\[
G_{u_1}=Q(u_1,u_2)\,G_{u_2},\qquad \|I-Q(u,\hat u)\|_{\mathcal{L}(Y)}\le \kappa(\rho),
\]
and compactness of $G_{u^\dagger}$ [1902.10596]. Under these assumptions, the iteration is a regularization method, terminates with stopping index
\[
N_\delta=\mathcal{O}(1+|\log\delta|),
\]
and satisfies
\[
u^\delta_{N_\delta}\to u^\dagger \quad \text{as } \delta\to 0
\]
under the stated nullspace condition [1902.10596].

In the PDE application, the generalized tangential cone condition and transfer bounds are quantified by
\[
\eta(\rho)=C_\Omega\left|\{\,|\bar y|\le C_F\rho\,\}\right|^{1/14},\qquad
\kappa(\rho)=C_*\left|\{\,|\bar y|\le C_F\rho\,\}\right|^{1/3},
\]
so sufficiently small measure of the near-zero set enables the analysis [1902.10596]. Numerical experiments reported an advantage over Bouligand–Landweber iteration, including logarithmic stopping-index growth and substantially lower total CPU time for small noise [1902.10596].

## 3. Zeroth-order Jacobian modeling for nonlinear least squares

A second line of work constructs Jacobian models from residual evaluations. In the orthogonal spherical smoothing approach, the residual map is $r(x)=(r_1(x),\dots,r_m(x))$ and the objective is again
\[
f(x)=\frac12\|r(x)\|^2=\frac12\sum_{i=1}^m r_i(x)^2,
\]
under a Lipschitz continuous Jacobian assumption
\[
\|J(x)-J(y)\|\le \kappa_{lj}\|x-y\|,\qquad \forall x,y\in\mathbb{R}^n
\]
and a Lipschitz residual assumption
\[
\|r(x)-r(y)\| \le \kappa_{lr}\|x-y\|,\qquad \forall x,y\in\mathbb{R}^n
\]
[2407.12542].

For each residual component, the spherically smoothed version is
\[
r_{s,i}(x)=\mathbb{E}_{v\sim U(\mathbb{B})}[r_i(x+\gamma v)],
\]
with gradient identity
\[
\nabla r_{s,i}(x)=\frac{n}{\gamma}\,\mathbb{E}_{u\sim U(\mathbb{S})}[r_i(x+\gamma u)\,u]
\]
[2407.12542]. To reduce variance, the method samples $b$ orthonormal directions $u_1,\dots,u_b$ uniformly from the Stiefel manifold
\[
St(n,b)=\{U\in\mathbb{R}^{n\times b}:U^\top U=I_b\},\qquad 1\le b\le n,
\]
and defines
\[
\nabla \tilde r_i(x)=\frac{n}{b}\sum_{j=1}^b \frac{r_i(x+\gamma u_j)-r_i(x)}{\gamma}\,u_j.
\]
Stacking these rows yields the approximate Jacobian $\tilde J(x)$, and the gradient model is
\[
\nabla \tilde f(x)=\tilde J(x)^\top r(x)
\]
[2407.12542]. Orthogonality gives
\[
\mathbb{E}[u_j u_j^\top]=\frac1n I,\qquad u_j^\top u_k=0\ \text{for }j\neq k,
\]
which removes cross terms in variance calculations [2407.12542].

The resulting derivative-free Levenberg–Marquardt step uses $J_{m_k}:=\tilde J(x_k)$ and solves
\[
(J_{m_k}^\top J_{m_k}+\lambda_k I)d=-J_{m_k}^\top r(x_k),\qquad
\lambda_k=\theta_k\|J_{m_k}^\top r(x_k)\|.
\]
Step acceptance is determined by
\[
\rho_k=\frac{\|r(x_k)\|^2-\|r(x_k+d_k)\|^2}{\|r(x_k)\|^2-\|r(x_k)+J_{m_k}d_k\|^2},
\]
with threshold parameters $0<p_0<p_1<p_2$, scaling factors $a_1>1>a_2>0$ satisfying $a_1a_2<1$, and lower bound $\theta_{\min}>0$ for the damping update [2407.12542]. The stopping condition is
\[
\|J_{m_k}^\top r(x_k)\|\le \epsilon_0.
\]

A central analytical notion is probabilistically first-order accurate gradient models. With event
\[
A_k=\left\{\|\nabla f_{M_k}(X_k)-\nabla f(X_k)\|\le \xi_1\|\nabla f(X_k)\|+\xi_2\|D_{k-1}\|\right\},
\]
the requirement is
\[
\mathbb{P}(A_k\mid \mathcal{F}_{k-1})\ge \alpha
\]
[2407.12542]. When $b\ge n/2$ and $\gamma_k=\|d_{k-1}\|$, the gradient models are $\alpha$-probabilistically $(\xi_1,\xi_2)$-first-order accurate, with explicit constants stated in the paper [2407.12542]. The high-probability global complexity bound is of order $O(\epsilon^{-2})$ for attaining a first-order stationary point [2407.12542].

A different zeroth-order construction appears for sparse nonlinear least-squares problems, where the Jacobian is unavailable or expensive but each row is assumed $s$-sparse [2507.06772]. At iteration $k$, one samples $p<n$ directions $v_k^1,\dots,v_k^p$, forms interpolation points $x_k+\sigma_k v_k^j$, and imposes linear interpolation constraints
\[
A(v_k)g_k^{(i)}=\sigma_k^{-1} y_k^{(i)}
\]
for each residual component [2507.06772]. The gradient row is then recovered by $\ell_1$ minimization,
\[
\min_{g\in\mathbb{R}^n}\ \|g\|_1 \quad \text{s.t.}\quad \|A(v_k)g-\sigma_k^{-1}y_k^{(i)}\|_2\le \xi_k,
\]
and the Jacobian model $m_k$ is assembled row-wise [2507.06772]. Sampling matrices are drawn from Gaussian, Bernoulli, or Bernoulli-like distributions so that Restricted Isometry Property conditions hold with high probability when $p$ is large enough [2507.06772]. The resulting model satisfies a probabilistic first-order accuracy bound of the form
\[
\|m_k-J(x_k)\|_F\le \frac{\varepsilon}{\theta_k}
\]
with probability at least $\beta$, and the associated derivative-free Levenberg–Marquardt algorithm converges globally almost surely in the sense that
\[
\mathbb{P}\!\left(\lim_{k\to\infty}\|J(X_k)^\top F(X_k)\|_2=0\right)=1
\]
under the stated assumptions [2507.06772].

## 4. Jacobian-free secant updates and quasi-Newton variants

A different derivative-free tradition keeps the Levenberg–Marquardt linear algebra but replaces repeated Jacobian computation by secant updates. In the Jacobian-free deterministic method for inverse problems, the residual vector is $r(x)\in\mathbb{R}^m$, possibly weighted by $W\succ0$, and the method solves the damped system
\[
(B_k^\top W B_k+\lambda_k \operatorname{diag}(B_k^\top W B_k))p_k=-B_k^\top W r_k,
\]
where $B_k\approx J_k$ is an approximate Jacobian updated by Broyden’s rank-one formula [2203.04138]:
\[
B_k=B_{k-1}+\frac{(y_k-B_{k-1}s_k)s_k^\top}{s_k^\top s_k},
\]
with
\[
s_k=x_k-x_{k-1},\qquad y_k=r(x_k)-r(x_{k-1}).
\]
The method uses only residual evaluations, no derivatives [2203.04138].

Global progress is controlled by Armijo backtracking. With gradient approximation $g_k=B_k^\top W r_k$, the line search seeks $\alpha_k$ such that
\[
f_W(x_k+\alpha_k p_k)\le f_W(x_k)+c_1\alpha_k g_k^\top p_k,
\]
where
\[
f_W(x)=\frac12 r(x)^\top W r(x)
\]
[2203.04138]. The paper emphasizes applications where analytic derivatives are impractical because the residual itself may contain a nested optimization problem or black-box computations [2203.04138]. Examples include thermodynamic model calibration through a thermodynamic code treated as a black box and calibration of a fission gas diffusion model through a COMSOL-based multiphysics code [2203.04138].

The earlier study on improvements to Levenberg–Marquardt likewise advocates Broyden updates to reduce Jacobian evaluations after one initial finite-difference Jacobian [1201.5885]. After an accepted step $\Delta\theta_i=\theta_i-\theta_{i-1}$ and residual change $\Delta r_i=r_i-r_{i-1}$, the update is
\[
J_i=J_{i-1}+\left[\frac{\Delta r_i-J_{i-1}\Delta\theta_i}{|\Delta\theta_i|^2}\right]\Delta\theta_i^\top.
\]
The recommendation is to refresh the Jacobian by a new finite-difference evaluation after one or two consecutive rejected steps, since repeated rank-one updates can drift [1201.5885].

These secant-based approaches preserve the least-squares structure more directly than generic derivative-free optimizers. The inverse-problem formulation in [2203.04138] explicitly contrasts Broyden with BFGS: $J$ is generally rectangular, making a Jacobian update more natural than a symmetric Hessian approximation. A plausible implication is that derivative-free Levenberg–Marquardt occupies an intermediate position between model-based derivative-free optimization and structured Gauss–Newton methods.

## 5. Step computation, damping, and acceptance mechanisms

Despite differences in Jacobian modeling, derivative-free Levenberg–Marquardt methods share a small set of algorithmic primitives: a damped linearized step, an acceptance test, and a rule for updating the damping or trust parameter.

In the Bouligand formulation, the regularization parameter is prescribed geometrically,
\[
\alpha_k=\alpha_0 r^k,\qquad r\in(0,1),
\]
with discrepancy-principle stopping and no additional line search or trust-region mechanism required in the analysis [1902.10596]. In the orthogonal spherical smoothing method, damping is tied to the model gradient norm,
\[
\lambda_k=\theta_k\|J_{m_k}^\top r(x_k)\|,
\]
and the parameter $\theta_k$ is increased or decreased according to the actual-to-predicted reduction ratio $\rho_k$ and threshold tests involving $p_0,p_1,p_2,a_1,a_2,\theta_{\min}$ [2407.12542]. In the sparse-Jacobian reconstruction method, the same structure appears with model $m_k$:
\[
(m_k^\top m_k+\lambda_k I)d_k=-m_k^\top F_k,\qquad \lambda_k=\theta_k\|m_k^\top F_k\|_2,
\]
and acceptance ratio
\[
\rho_k=\frac{\|F_k\|_2^2-\|F(x_k+d_k)\|_2^2}{\|F_k\|_2^2-\|F_k+m_kd_k\|_2^2}
\]
[2507.06772].

The geodesic-acceleration variant enriches the standard Levenberg–Marquardt step by a second-order correction [1201.5885]. Writing
\[
\delta\theta_1=-(J^\top J+\lambda D^\top D)^{-1}J^\top r,
\]
the correction is
\[
\delta\theta_2\approx -\frac12 (J^\top J+\lambda D^\top D)^{-1}J^\top r'',
\]
with directional second derivative $r''$ estimated by finite differences along $\delta\theta_1$ [1201.5885]. The truncation safeguard
\[
\frac{2|\delta\theta_2|}{|\delta\theta_1|}\le \alpha
\]
controls the step, with $\alpha\approx 0.75$ reported as a good default and smaller values improving robustness on difficult cases [1201.5885].

The same paper also studies controlled uphill acceptance. If $\beta_i$ is the cosine between the current and previous accepted Levenberg–Marquardt velocity parts,
\[
\beta_i=\cos(\delta\theta_1^{new},\delta\theta_1^{old}),
\]
then uphill steps may be accepted when
\[
(1-\beta_i)^b C_{i+1}\le C_i
\]
or, more conservatively,
\[
(1-\beta_i)^b C_{i+1}\le \min(C_1,\dots,C_i),
\]
with $b\in\{1,2\}$ [1201.5885]. This mechanism is intended for narrow valleys in which downhill-only acceptance can force tiny steps.

The diversity of acceptance rules is one of the main distinctions among derivative-free Levenberg–Marquardt methods. Some formulations treat damping primarily as regularization [1902.10596]; others operate in a trust-region-like regime with predicted reduction ratios [2407.12542, 2507.06772]; still others combine Levenberg–Marquardt linearization with line search and Armijo sufficient decrease [2203.04138]. This suggests that the derivative-free aspect concerns model construction more than globalization strategy.

## 6. Computational behavior, applications, and comparisons

The computational profile of derivative-free Levenberg–Marquardt depends strongly on how the Jacobian surrogate is built. Orthogonal spherical smoothing requires $b$ residual evaluations at perturbed points $x_k+\gamma_k u_j$ plus possibly one at $x_k$, followed by linear algebra for an $n\times n$ positive definite system [2407.12542]. The sparse-Jacobian method requires evaluation of $F$ at $p$ shifted points, solution of $m$ independent $\ell_1$ problems, and then a sparse Levenberg–Marquardt solve; when $p\approx Cs\log(n/s)$, this can be substantially cheaper than finite-difference Jacobians requiring $O(n)$ evaluations [2507.06772]. Broyden-based schemes need one initial Jacobian or initial secant seed and then only residual evaluations per iteration, with occasional refreshes to counter drift [1201.5885, 2203.04138].

The literature emphasizes different application domains. The Bouligand method is developed for ill-posed inverse problems and demonstrated on an inverse source problem for a non-smooth semilinear elliptic PDE [1902.10596]. The Jacobian-free deterministic method is designed for inverse problems involving black-box simulators and nested optimization, including thermodynamic calibration and nuclear fuel performance modeling [2203.04138]. Orthogonal spherical smoothing addresses general nonlinear least squares with probabilistic first-order accuracy and high-probability complexity guarantees [2407.12542]. The sparse-Jacobian method targets sparse nonlinear least-squares problems in which underlying Jacobian sparsity can be exploited both in model construction and linear algebra [2507.06772]. The geodesic-acceleration paper studies broad nonlinear least-squares benchmarks, including MINPACK-2, NIST datasets, and large real-world models [1201.5885].

Several comparative claims recur. In the non-smooth inverse-problem setting, Bouligand–Levenberg–Marquardt is reported to outperform Bouligand–Landweber in stopping index and total CPU time, especially for small noise [1902.10596]. In the smoothing-based framework, orthogonal spherical smoothing is presented as a low-variance alternative to coordinate-wise forward differences, with performance profiles showing that DFLM-OSS variants solved roughly $94\%$ of problems at tolerance $10^{-3}$ and around $90\%$ at tolerance $10^{-5}$, whereas DFLM-FD solved fewer [2407.12542]. In the geodesic-acceleration study, acceleration improved speed and robustness, with speedups up to $70\times$ fewer Jacobian evaluations on some problems and typical improvements of $2$–$10\times$; bold acceptance sometimes reduced iteration counts by up to $\sim30\times$ but could reduce robustness if used alone [1201.5885]. In sparse least squares, DFLM-SNLS is reported to outperform MATLAB baselines in function evaluations, particularly at higher dimension and higher sparsity [2507.06772].

A persistent misconception is that derivative-free Levenberg–Marquardt means the same thing as generic zeroth-order optimization. The papers surveyed here do not support that interpretation. Even the most function-evaluation-based variants still exploit residual structure, linearized least-squares models, and damped normal equations [2407.12542, 2507.06772]. The Bouligand method is derivative-free only in the classical smooth sense, not in the sense of dispensing with linearized models altogether [1902.10596]. Conversely, Broyden-based methods are Jacobian-free after initialization, but they still maintain an evolving Jacobian approximation rather than abandoning derivative structure [1201.5885, 2203.04138].

## 7. Theoretical guarantees and open methodological distinctions

Theoretical guarantees differ sharply across formulations. For non-smooth ill-posed inverse problems, the Bouligand method is proved to be an iterative regularization scheme under the generalized tangential cone condition, transfer-operator bounds, compactness, geometric parameter decay, and discrepancy stopping [1902.10596]. The convergence result is expressed in the regularization sense:
\[
u^\delta_{N_\delta}\to u^\dagger \quad \text{as } \delta\to 0,
\]
with stopping index complexity
\[
N_\delta=\mathcal{O}(1+|\log\delta|)
\]
[1902.10596].

For orthogonal spherical smoothing, the guarantee is probabilistic first-order model accuracy and high-probability global complexity. Under the stated smoothness and boundedness assumptions, the algorithm attains worst-case global complexity of order $O(\epsilon^{-2})$ to reach a first-order stationary point, with an explicit tail bound on the probability of failure [2407.12542]. For sparse nonlinear least squares, the analysis uses compressed sensing, RIP-based sparse recovery, and a random-walk submartingale argument to show almost-sure convergence of the first-order optimality measure:
\[
\mathbb{P}\!\left(\lim_{k\to\infty}\|J(X_k)^\top F(X_k)\|_2=0\right)=1
\]
[2507.06772].

By contrast, the Broyden-based inverse-problem method presents global convergence in terms of standard assumptions combined with Armijo sufficient decrease and increasingly accurate Jacobian approximations in the limit [2203.04138]. The geodesic-acceleration paper is primarily an algorithmic and empirical study of improvements to Levenberg–Marquardt rather than a full derivative-free global complexity theory [1201.5885].

These distinctions matter because “derivative-free Levenberg–Marquardt” can refer to at least three analytically different objects: a regularization method for non-smooth inverse operators [1902.10596], a probabilistic zeroth-order algorithm for nonlinear least squares [2407.12542, 2507.06772], or a secant-updated Jacobian-free practical solver for black-box inverse problems [1201.5885, 2203.04138]. A plausible implication is that no single convergence framework subsumes all such methods without significant abstraction. The topic is therefore best understood as a family of Levenberg–Marquardt adaptations united by damping and least-squares structure, but separated by the kind of derivative information they replace and by the analytical machinery used to justify the replacement.

Source: https://www.emergentmind.com/topics/derivative-free-levenberg-marquardt-algorithm