Derivative-Free Levenberg–Marquardt

• Derivative-Free Levenberg–Marquardt algorithms are variants that bypass analytical Jacobians using surrogate models like Bouligand subderivatives, stochastic approximations, and secant updates.
• They maintain the regularized Gauss–Newton structure with damping strategies and acceptance rules to robustly solve nonlinear least-squares and inverse problems.
• Distinct methods—including orthogonal spherical smoothing and compressed sensing for sparse problems—offer practical improvements backed by theoretical convergence guarantees.

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 (Clason et al., 2019); stochastic Jacobian approximation by orthogonal spherical smoothing for nonlinear least squares (Chen et al., 2024); Jacobian-free schemes based on Broyden rank-one updates with Armijo line search for inverse problems (Piro et al., 2022); derivative-free enhancements using geodesic acceleration, uphill acceptance, and Broyden updates (Transtrum et al., 2012); and compressed-sensing-based sparse Jacobian reconstruction for sparse nonlinear least-squares problems (Feng et al., 9 Jul 2025). 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 (Chen et al., 2024). 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 (Clason et al., 2019).

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 (Clason et al., 2019). In another usage, it means Jacobian approximation from function values alone, either through orthogonal spherical smoothing (Chen et al., 2024), Broyden updates (Transtrum et al., 2012, Piro et al., 2022), or sparse recovery from interpolation constraints (Feng et al., 9 Jul 2025). 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 (Chen et al., 2024, Feng et al., 9 Jul 2025). In ill-posed inverse problems, this damping is explicitly interpreted as iterative regularization (Clason et al., 2019).

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 (Clason et al., 2019). 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 $f(x)=\frac12\|r(x)\|^2,$0 satisfying

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

The inverse problem is to recover $f(x)=\frac12\|r(x)\|^2,$2 from $f(x)=\frac12\|r(x)\|^2,$3 (Clason et al., 2019).

The motivating PDE is the non-smooth semilinear elliptic problem

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

with $f(x)=\frac12\|r(x)\|^2,$5 and $f(x)=\frac12\|r(x)\|^2,$6, where $f(x)=\frac12\|r(x)\|^2,$7 for $f(x)=\frac12\|r(x)\|^2,$8 is bounded with Lipschitz boundary (Clason et al., 2019). In this case the forward map fails to be Gâteaux differentiable at $f(x)=\frac12\|r(x)\|^2,$9 if the zero level set $\nabla f(x)=J(x)^\top r(x)$0 has positive measure; precisely, Gâteaux differentiability holds at $\nabla f(x)=J(x)^\top r(x)$1 if and only if $\nabla f(x)=J(x)^\top r(x)$2 (Clason et al., 2019).

The substitute for the missing derivative is the Bouligand subdifferential

$\nabla f(x)=J(x)^\top r(x)$3

where $\nabla f(x)=J(x)^\top r(x)$4 is the set of Gâteaux differentiability points (Clason et al., 2019). For the PDE application, an explicit Bouligand element $\nabla f(x)=J(x)^\top r(x)$5 is defined by solving the linearized PDE

$\nabla f(x)=J(x)^\top r(x)$6

with $\nabla f(x)=J(x)^\top r(x)$7 (Clason et al., 2019). The active-set indicator $\nabla f(x)=J(x)^\top r(x)$8 incorporates the non-smoothness directly.

The corresponding Levenberg–Marquardt step replaces the classical derivative by $\nabla f(x)=J(x)^\top r(x)$9 and solves

$J(x)\in\mathbb{R}^{m\times n}$0

with a geometrically decreasing a priori parameter rule

$J(x)\in\mathbb{R}^{m\times n}$1

and Morozov discrepancy stopping

$J(x)\in\mathbb{R}^{m\times n}$2

where $J(x)\in\mathbb{R}^{m\times n}$3 is minimal (Clason et al., 2019).

The analysis uses a generalized tangential cone condition,

$J(x)\in\mathbb{R}^{m\times n}$4

with $J(x)\in\mathbb{R}^{m\times n}$5, together with transfer operators $J(x)\in\mathbb{R}^{m\times n}$6 satisfying

$J(x)\in\mathbb{R}^{m\times n}$7

and compactness of $J(x)\in\mathbb{R}^{m\times n}$8 (Clason et al., 2019). Under these assumptions, the iteration is a regularization method, terminates with stopping index

$J(x)\in\mathbb{R}^{m\times n}$9

and satisfies

$F:D(F)\subset U\to Y$0

under the stated nullspace condition (Clason et al., 2019).

In the PDE application, the generalized tangential cone condition and transfer bounds are quantified by

$F:D(F)\subset U\to Y$1

so sufficiently small measure of the near-zero set enables the analysis (Clason et al., 2019). Numerical experiments reported an advantage over Bouligand–Landweber iteration, including logarithmic stopping-index growth and substantially lower total CPU time for small noise (Clason et al., 2019).

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 $F:D(F)\subset U\to Y$2 and the objective is again

$F:D(F)\subset U\to Y$3

under a Lipschitz continuous Jacobian assumption

$F:D(F)\subset U\to Y$4

and a Lipschitz residual assumption

$F:D(F)\subset U\to Y$5

(Chen et al., 2024).

For each residual component, the spherically smoothed version is

$F:D(F)\subset U\to Y$6

with gradient identity

$F:D(F)\subset U\to Y$7

(Chen et al., 2024). To reduce variance, the method samples $F:D(F)\subset U\to Y$8 orthonormal directions $F:D(F)\subset U\to Y$9 uniformly from the Stiefel manifold

$(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$0

and defines

$(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$1

Stacking these rows yields the approximate Jacobian $(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$2, and the gradient model is

$(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$3

(Chen et al., 2024). Orthogonality gives

$(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$4

which removes cross terms in variance calculations (Chen et al., 2024).

The resulting derivative-free Levenberg–Marquardt step uses $(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$5 and solves

$(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$6

Step acceptance is determined by

$(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$7

with threshold parameters $(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$8, scaling factors $(J_k^\top J_k+\lambda_k I)d_k=-J_k^\top r(x_k)$9 satisfying $U$0, and lower bound $U$1 for the damping update (Chen et al., 2024). The stopping condition is

$U$2

A central analytical notion is probabilistically first-order accurate gradient models. With event

$U$3

the requirement is

$U$4

(Chen et al., 2024). When $U$5 and $U$6, the gradient models are $U$7-probabilistically $U$8-first-order accurate, with explicit constants stated in the paper (Chen et al., 2024). The high-probability global complexity bound is of order $U$9 for attaining a first-order stationary point (Chen et al., 2024).

A different zeroth-order construction appears for sparse nonlinear least-squares problems, where the Jacobian is unavailable or expensive but each row is assumed $Y$0-sparse (Feng et al., 9 Jul 2025). At iteration $Y$1, one samples $Y$2 directions $Y$3, forms interpolation points $Y$4, and imposes linear interpolation constraints

$Y$5

for each residual component (Feng et al., 9 Jul 2025). The gradient row is then recovered by $Y$6 minimization,

$Y$7

and the Jacobian model $Y$8 is assembled row-wise (Feng et al., 9 Jul 2025). Sampling matrices are drawn from Gaussian, Bernoulli, or Bernoulli-like distributions so that Restricted Isometry Property conditions hold with high probability when $Y$9 is large enough (Feng et al., 9 Jul 2025). The resulting model satisfies a probabilistic first-order accuracy bound of the form

$F:D(F)\subset U\to Y$0

with probability at least $F:D(F)\subset U\to Y$1, and the associated derivative-free Levenberg–Marquardt algorithm converges globally almost surely in the sense that

$F:D(F)\subset U\to Y$2

under the stated assumptions (Feng et al., 9 Jul 2025).

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 $F:D(F)\subset U\to Y$3, possibly weighted by $F:D(F)\subset U\to Y$4, and the method solves the damped system

$F:D(F)\subset U\to Y$5

where $F:D(F)\subset U\to Y$6 is an approximate Jacobian updated by Broyden’s rank-one formula (Piro et al., 2022): $F:D(F)\subset U\to Y$7 with

$F:D(F)\subset U\to Y$8

The method uses only residual evaluations, no derivatives (Piro et al., 2022).

Global progress is controlled by Armijo backtracking. With gradient approximation $F:D(F)\subset U\to Y$9, the line search seeks $y^\dagger=F(u^\dagger)$0 such that

$y^\dagger=F(u^\dagger)$1

where

$y^\dagger=F(u^\dagger)$2

(Piro et al., 2022). The paper emphasizes applications where analytic derivatives are impractical because the residual itself may contain a nested optimization problem or black-box computations (Piro et al., 2022). 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 (Piro et al., 2022).

The earlier study on improvements to Levenberg–Marquardt likewise advocates Broyden updates to reduce Jacobian evaluations after one initial finite-difference Jacobian (Transtrum et al., 2012). After an accepted step $y^\dagger=F(u^\dagger)$3 and residual change $y^\dagger=F(u^\dagger)$4, the update is

$y^\dagger=F(u^\dagger)$5

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 (Transtrum et al., 2012).

These secant-based approaches preserve the least-squares structure more directly than generic derivative-free optimizers. The inverse-problem formulation in (Piro et al., 2022) explicitly contrasts Broyden with BFGS: $y^\dagger=F(u^\dagger)$6 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,

$y^\dagger=F(u^\dagger)$7

with discrepancy-principle stopping and no additional line search or trust-region mechanism required in the analysis (Clason et al., 2019). In the orthogonal spherical smoothing method, damping is tied to the model gradient norm,

$y^\dagger=F(u^\dagger)$8

and the parameter $y^\dagger=F(u^\dagger)$9 is increased or decreased according to the actual-to-predicted reduction ratio $f(x)=\frac12\|r(x)\|^2,$00 and threshold tests involving $f(x)=\frac12\|r(x)\|^2,$01 (Chen et al., 2024). In the sparse-Jacobian reconstruction method, the same structure appears with model $f(x)=\frac12\|r(x)\|^2,$02: $f(x)=\frac12\|r(x)\|^2,$03 and acceptance ratio

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

(Feng et al., 9 Jul 2025).

The geodesic-acceleration variant enriches the standard Levenberg–Marquardt step by a second-order correction (Transtrum et al., 2012). Writing

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

the correction is

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

with directional second derivative $f(x)=\frac12\|r(x)\|^2,$07 estimated by finite differences along $f(x)=\frac12\|r(x)\|^2,$08 (Transtrum et al., 2012). The truncation safeguard

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

controls the step, with $f(x)=\frac12\|r(x)\|^2,$10 reported as a good default and smaller values improving robustness on difficult cases (Transtrum et al., 2012).

The same paper also studies controlled uphill acceptance. If $f(x)=\frac12\|r(x)\|^2,$11 is the cosine between the current and previous accepted Levenberg–Marquardt velocity parts,

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

then uphill steps may be accepted when

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

or, more conservatively,

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

with $f(x)=\frac12\|r(x)\|^2,$15 (Transtrum et al., 2012). 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 (Clason et al., 2019); others operate in a trust-region-like regime with predicted reduction ratios (Chen et al., 2024, Feng et al., 9 Jul 2025); still others combine Levenberg–Marquardt linearization with line search and Armijo sufficient decrease (Piro et al., 2022). 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 $f(x)=\frac12\|r(x)\|^2,$16 residual evaluations at perturbed points $f(x)=\frac12\|r(x)\|^2,$17 plus possibly one at $f(x)=\frac12\|r(x)\|^2,$18, followed by linear algebra for an $f(x)=\frac12\|r(x)\|^2,$19 positive definite system (Chen et al., 2024). The sparse-Jacobian method requires evaluation of $f(x)=\frac12\|r(x)\|^2,$20 at $f(x)=\frac12\|r(x)\|^2,$21 shifted points, solution of $f(x)=\frac12\|r(x)\|^2,$22 independent $f(x)=\frac12\|r(x)\|^2,$23 problems, and then a sparse Levenberg–Marquardt solve; when $f(x)=\frac12\|r(x)\|^2,$24, this can be substantially cheaper than finite-difference Jacobians requiring $f(x)=\frac12\|r(x)\|^2,$25 evaluations (Feng et al., 9 Jul 2025). Broyden-based schemes need one initial Jacobian or initial secant seed and then only residual evaluations per iteration, with occasional refreshes to counter drift (Transtrum et al., 2012, Piro et al., 2022).

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 (Clason et al., 2019). 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 (Piro et al., 2022). Orthogonal spherical smoothing addresses general nonlinear least squares with probabilistic first-order accuracy and high-probability complexity guarantees (Chen et al., 2024). 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 (Feng et al., 9 Jul 2025). The geodesic-acceleration paper studies broad nonlinear least-squares benchmarks, including MINPACK-2, NIST datasets, and large real-world models (Transtrum et al., 2012).

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 (Clason et al., 2019). 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 $f(x)=\frac12\|r(x)\|^2,$26 of problems at tolerance $f(x)=\frac12\|r(x)\|^2,$27 and around $f(x)=\frac12\|r(x)\|^2,$28 at tolerance $f(x)=\frac12\|r(x)\|^2,$29, whereas DFLM-FD solved fewer (Chen et al., 2024). In the geodesic-acceleration study, acceleration improved speed and robustness, with speedups up to $f(x)=\frac12\|r(x)\|^2,$30 fewer Jacobian evaluations on some problems and typical improvements of $f(x)=\frac12\|r(x)\|^2,$31–$f(x)=\frac12\|r(x)\|^2,$32; bold acceptance sometimes reduced iteration counts by up to $f(x)=\frac12\|r(x)\|^2,$33 but could reduce robustness if used alone (Transtrum et al., 2012). In sparse least squares, DFLM-SNLS is reported to outperform MATLAB baselines in function evaluations, particularly at higher dimension and higher sparsity (Feng et al., 9 Jul 2025).

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 (Chen et al., 2024, Feng et al., 9 Jul 2025). The Bouligand method is derivative-free only in the classical smooth sense, not in the sense of dispensing with linearized models altogether (Clason et al., 2019). Conversely, Broyden-based methods are Jacobian-free after initialization, but they still maintain an evolving Jacobian approximation rather than abandoning derivative structure (Transtrum et al., 2012, Piro et al., 2022).

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 (Clason et al., 2019). The convergence result is expressed in the regularization sense: $f(x)=\frac12\|r(x)\|^2,$34 with stopping index complexity

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

(Clason et al., 2019).

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 $f(x)=\frac12\|r(x)\|^2,$36 to reach a first-order stationary point, with an explicit tail bound on the probability of failure (Chen et al., 2024). 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: $f(x)=\frac12\|r(x)\|^2,$37 (Feng et al., 9 Jul 2025).

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 (Piro et al., 2022). 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 (Transtrum et al., 2012).

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 (Clason et al., 2019), a probabilistic zeroth-order algorithm for nonlinear least squares (Chen et al., 2024, Feng et al., 9 Jul 2025), or a secant-updated Jacobian-free practical solver for black-box inverse problems (Transtrum et al., 2012, Piro et al., 2022). 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.