LMM Subgradient Method for Composite Optimization

• The method is a preconditioned first-order approach for composite optimization, using LM-style regularization to counteract slow subgradient rates.
• It achieves local linear convergence under mild regularity and alignment conditions, with convergence driven solely by the conditioning of the convex function.
• The algorithm shows practical efficiency in ill-conditioned or overparameterized scenarios, such as nonnegative least squares, matrix sensing, and tensor factorization.

The Levenberg–Morrison–Marquardt subgradient method is a preconditioned first-order method for composite optimization problems of the form

$\min_{x\in\mathbb R^d} f(x):=h(F(x)),$

where $F:\mathbb R^d\to\mathbb R^m$ is $C^1$–smooth and $h:\mathbb R^m\to\mathbb R$ is convex and possibly nonsmooth. Its defining feature is an LM-style preconditioner,

$$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$

inserted between a subgradient of $h$ and the parameter update. In the formulation introduced in "Preconditioned subgradient method for composite optimization: overparameterization and fast convergence," the method is designed to overcome the slow sublinear behavior of the vanilla subgradient method in ill-conditioned or overparameterized composite models, and it is proved to converge linearly under mild regularity conditions at a rate determined solely by the convex function (Díaz et al., 15 Sep 2025).

1. Optimization model and notation

The method is posed for composite minimization

$\min_{x\in\mathbb R^d} f(x):=h(F(x)),$

with $F:\mathbb R^d\to\mathbb R^m$ of class $C^1$ and $h$ convex, possibly nonsmooth. The convex subdifferential of $F:\mathbb R^d\to\mathbb R^m$0 at $F:\mathbb R^d\to\mathbb R^m$1 is denoted $F:\mathbb R^d\to\mathbb R^m$2. The Jacobian is written $F:\mathbb R^d\to\mathbb R^m$3, with adjoint $F:\mathbb R^d\to\mathbb R^m$4. For a given point $F:\mathbb R^d\to\mathbb R^m$5, $F:\mathbb R^d\to\mathbb R^m$6 denotes the orthogonal projector onto the range of $F:\mathbb R^d\to\mathbb R^m$7, and

$F:\mathbb R^d\to\mathbb R^m$8

The analysis fixes a solution $F:\mathbb R^d\to\mathbb R^m$9 and writes $C^1$0. The regularity assumptions are organized around three components. Assumption A requires $C^1$1 to be $C^1$2–Lipschitz:

$C^1$3

Assumption B concerns the geometry of the image of $C^1$4 near $C^1$5. Under strong alignment, the rank of $C^1$6 is locally constant and equal to $C^1$7, so $C^1$8 is a $C^1$9–manifold and the tangent-manifold error satisfies

$h:\mathbb R^m\to\mathbb R$0

Under weak alignment, appropriate singular directions of $h:\mathbb R^m\to\mathbb R$1 capture the dominant part of $h:\mathbb R^m\to\mathbb R$2, while a lower singular-value bound scales with $h:\mathbb R^m\to\mathbb R$3.

The convex term $h:\mathbb R^m\to\mathbb R$4 is handled through two alternative hypotheses. Assumption $h:\mathbb R^m\to\mathbb R$5 covers nonsmooth sharp objectives: $h:\mathbb R^m\to\mathbb R$6 is convex on $h:\mathbb R^m\to\mathbb R$7, has a unique minimizer $h:\mathbb R^m\to\mathbb R$8, is $h:\mathbb R^m\to\mathbb R$9-sharp,

$$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$0

and satisfies restricted Lipschitz conditions involving $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$1 and $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$2. Assumption $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$3 covers smooth objectives with quadratic growth: $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$4 has $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$5–quadratic growth, restricted $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$6–smoothness, and analogous compatibility bounds with the projected and damped geometry (Díaz et al., 15 Sep 2025).

2. Algorithmic update and LM preconditioning

Algorithm 1 in the cited work takes positive sequences $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$7 and $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$8 and, at each iteration, chooses a subgradient $$(\nabla F(x)^\top \nabla F(x)+\lambda I)^{-1}\nabla F(x)^\top,$$9, solves

$h$0

and updates

$h$1

Equivalently,

$h$2

In the image space, the step is approximately

$h$3

with a second-order Taylor error of order $h$4. More precisely, the update in $h$5-space satisfies

$h$6

and the remainder is controlled by

$h$7

The role of the preconditioner is central. The term $h$8 is an ordinary subgradient of $h$9. The premultiplier

$\min_{x\in\mathbb R^d} f(x):=h(F(x)),$0

normalizes directions by the local Jacobian of $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$1, thereby undoing ill-conditioning or rank-deficiency of $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$2. When $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$3 and $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$4 has full column rank, $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$5, the projector onto $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$6. The eigenvalues of $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$7 lie in $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$8, and the eigenvalues of $\min_{x\in\mathbb R^d} f(x):=h(F(x)),$9 are at most $F:\mathbb R^d\to\mathbb R^m$0 (Díaz et al., 15 Sep 2025).

3. Alignment assumptions and local geometry

The convergence theory is built on an interaction between the local geometry of $F:\mathbb R^d\to\mathbb R^m$1 and the growth properties of $F:\mathbb R^d\to\mathbb R^m$2. Under strong alignment, the image of $F:\mathbb R^d\to\mathbb R^m$3 behaves locally as a $F:\mathbb R^d\to\mathbb R^m$4–manifold of constant rank. In that regime, the projected error onto the normal space is quadratic in $F:\mathbb R^d\to\mathbb R^m$5, which makes the geometry close to exact parameterization. Under weak alignment, the framework allows overparameterization: even if rank-deficiency or degeneracy is present, there exist singular directions of $F:\mathbb R^d\to\mathbb R^m$6 that capture most of the signal in $F:\mathbb R^d\to\mathbb R^m$7, while the residual component can be made arbitrarily small relative to $F:\mathbb R^d\to\mathbb R^m$8.

The assumptions on $F:\mathbb R^d\to\mathbb R^m$9 are correspondingly restricted to the image geometry. In the nonsmooth case, sharpness and restricted Lipschitzness are not imposed on the whole ambient space $C^1$0 in an unrestricted way; instead, they are formulated through the projected component $C^1$1 and the damped residual $C^1$2. In the smooth case, the analogous hypotheses are quadratic growth and restricted smoothness. This suggests that the method isolates the conditioning of the convex outer function from the possible ill-conditioning of the smooth parametrization map.

The paper emphasizes two regimes. In an exact-rank regime, the method reduces to a Riemannian subgradient on the image-manifold. In overparameterized or ill-conditioned regimes, the damping parameter $C^1$3 ensures invertibility and stability. The resulting theory is explicitly local, centered at a fixed solution $C^1$4 and its image $C^1$5 (Díaz et al., 15 Sep 2025).

4. Hyperparameter configurations and convergence guarantees

Two hyperparameter configurations are stated. The Polyak-type configuration requires access to $C^1$6 and $C^1$7 and uses

$C^1$8

with $C^1$9 between $h$0 and a multiple thereof. The geometric-decay configuration avoids knowledge of $h$1: for nonsmooth $h$2, one sets $h$3 and $h$4; for smooth $h$5, one sets $h$6 and $h$7, with $h$8.

The principal theorem in the nonsmooth weak-alignment setting states local linear convergence under Assumptions A, B(weak), and $h$9. If $F:\mathbb R^d\to\mathbb R^m$00 satisfies

$F:\mathbb R^d\to\mathbb R^m$01

and Algorithm 1 is run with the Polyak configuration, with $F:\mathbb R^d\to\mathbb R^m$02, then for all $F:\mathbb R^d\to\mathbb R^m$03,

$F:\mathbb R^d\to\mathbb R^m$04

The proof proceeds by comparing the true iterate with the one-step linearized iterate

$F:\mathbb R^d\to\mathbb R^m$05

Taylor’s theorem and the Lipschitzness of $F:\mathbb R^d\to\mathbb R^m$06 control the deviation $F:\mathbb R^d\to\mathbb R^m$07, while restricted sharpness and restricted Lipschitzness yield a contraction for the linearized step. Under strong alignment, the exponent improves from $F:\mathbb R^d\to\mathbb R^m$08 to $F:\mathbb R^d\to\mathbb R^m$09.

The smooth case is described as analogous, replacing sharpness by quadratic growth and replacing the ratio $F:\mathbb R^d\to\mathbb R^m$10 by the corresponding $F:\mathbb R^d\to\mathbb R^m$11 dependence. The central conclusion is that the local linear rate is governed solely by the conditioning of the convex function $F:\mathbb R^d\to\mathbb R^m$12, rather than by the ill-conditioning or overparameterization of $F:\mathbb R^d\to\mathbb R^m$13 (Díaz et al., 15 Sep 2025).

5. Relation to classical subgradient methods and nearby LM schemes

The comparison with the classical subgradient method is explicit. A standard update takes the form

$F:\mathbb R^d\to\mathbb R^m$14

Even with Polyak stepsizes,

$F:\mathbb R^d\to\mathbb R^m$15

linear convergence requires $F:\mathbb R^d\to\mathbb R^m$16 itself to be sharp and Lipschitz on $F:\mathbb R^d\to\mathbb R^m$17. For composite objectives $F:\mathbb R^d\to\mathbb R^m$18, this requirement is restrictive because the composition may be ill-conditioned or rank-deficient near a solution. In that setting, Gauss–Newton subgradient with $F:\mathbb R^d\to\mathbb R^m$19 and the unpreconditioned subgradient can stagnate at sublinear rate or diverge, whereas LMM uses the damped LM preconditioner to recover linear rates controlled by $F:\mathbb R^d\to\mathbb R^m$20 (Díaz et al., 15 Sep 2025).

A related but distinct line of work appears in "A Levenberg-Marquardt Method for Nonsmooth Regularized Least Squares," which studies problems of the form

$F:\mathbb R^d\to\mathbb R^m$21

with $F:\mathbb R^d\to\mathbb R^m$22 proper, lower-semicontinuous, and possibly nonsmooth and nonconvex. There, the LM step is obtained by minimizing a regularized linear least-squares model plus a model of $F:\mathbb R^d\to\mathbb R^m$23, approximately solved by a first-order method such as proximal gradient, and both a regularization variant and a trust-region variant are shown to have worst-case complexity $F:\mathbb R^d\to\mathbb R^m$24 for reducing a stationarity measure below $F:\mathbb R^d\to\mathbb R^m$25 (Aravkin et al., 2023). This suggests a methodological neighborhood rather than an identity: both frameworks use LM regularization around Jacobian-linearized models, but the LMM subgradient method is formulated as a preconditioned subgradient method for $F:\mathbb R^d\to\mathbb R^m$26, while the nonsmooth regularized least-squares method targets $F:\mathbb R^d\to\mathbb R^m$27 and analyzes model-based outer iterations with inner proximal solves.

6. Applications and computational behavior

The LMM method is applied to several classes of problems. For squared-variable reformulations of nonnegative least squares, with $F:\mathbb R^d\to\mathbb R^m$28 and $F:\mathbb R^d\to\mathbb R^m$29 or $F:\mathbb R^d\to\mathbb R^m$30, weak alignment for $F:\mathbb R^d\to\mathbb R^m$31 and restricted sharpness or Lipschitzness for $F:\mathbb R^d\to\mathbb R^m$32 can be verified via $F:\mathbb R^d\to\mathbb R^m$33. In that setting, LMM converges linearly at rate $F:\mathbb R^d\to\mathbb R^m$34 in practice and outperforms both Gauss–Newton-subgradient and the unpreconditioned subgradient.

For low-rank matrix sensing, the method is analyzed for positive-semidefinite factors $F:\mathbb R^d\to\mathbb R^m$35 and asymmetric factors $F:\mathbb R^d\to\mathbb R^m$36. Under RIP, one obtains strong or weak alignment for $F:\mathbb R^d\to\mathbb R^m$37 and restricted smoothness or growth for $F:\mathbb R^d\to\mathbb R^m$38 or for $F:\mathbb R^d\to\mathbb R^m$39 losses. The resulting local linear convergence is described as independent of the condition number of $F:\mathbb R^d\to\mathbb R^m$40 and governed solely by $F:\mathbb R^d\to\mathbb R^m$41–RIP constants. For asymmetric factors, local weak alignment together with balanced spectral initialization again yields linear rates for $F:\mathbb R^d\to\mathbb R^m$42 and $F:\mathbb R^d\to\mathbb R^m$43 losses, including robustness to gross outliers.

For canonical-polyadic tensor factorization, the paper considers symmetric CP,

$F:\mathbb R^d\to\mathbb R^m$44

and asymmetric CP,

$F:\mathbb R^d\to\mathbb R^m$45

Local strong alignment is established near a full-rank CP factorization, and $F:\mathbb R^d\to\mathbb R^m$46–quadratic growth together with Lipschitz properties are verified for $F:\mathbb R^d\to\mathbb R^m$47. Consequently, LMM enjoys local linear convergence with $F:\mathbb R^d\to\mathbb R^m$48 loss. The reported numerical experiments on nonnegative least squares, matrix sensing, and tensor factorization confirm consistent linear convergence to exact solutions, improved performance over unpreconditioned subgradient and Gauss–Newton subgradient in ill-conditioned or overparameterized regimes, and robustness to gross outliers when using the $F:\mathbb R^d\to\mathbb R^m$49 loss (Díaz et al., 15 Sep 2025).

A nearby empirical pattern appears in the nonsmooth regularized least-squares LM literature: for a group-lasso basis-pursuit denoise problem, a nonlinear support vector machine, and a FitzHugh–Nagumo neuron-firing inverse problem, the LM variants required fewer outer iterations than adaptive proximal gradient and a quasi-Newton trust-region method, although at the expense of many more inner proximal-gradient steps (Aravkin et al., 2023). A plausible implication is that LM-style curvature information can substantially reduce outer-iteration counts, while leaving the efficiency of the inner linearized composite solve as a critical implementation bottleneck.