PhD research · SIAM Journal on Optimization, 2026

R2N: fast steps for problems with kinks and jumps

A quasi-Newton method, with convergence guarantees and a Julia implementation, for optimization problems that combine a smooth data-fitting term with a nonsmooth regularizer such as a sparsity or low-rank penalty.

Paper
A proximal modified quasi-Newton method for nonsmooth regularized optimization
Youssef Diouane, Mohamed Laghdaf Habiboullah, Dominique Orban
SIAM Journal on Optimization, 36(2):534–563, 2026
In one minute

Many problems in signal processing, statistics and machine learning ask for a solution that fits the data and has a simple structure: few nonzero entries, low rank, or values in a given range. They are written as the minimization of \(f + h\), where \(f\) measures the fit and the regularizer \(h\) rewards structure. The regularizer is typically nonsmooth, sometimes even discontinuous, and both terms may be nonconvex. R2N computes each step from a quasi-Newton model of \(f\), a model of \(h\) and an adaptive penalty on long steps. It converges under weaker assumptions than competing methods, and it can directly solve problems, such as minimum-rank matrix completion, that trust-region methods currently cannot.

The problem

Suppose you want to recover a signal from a few noisy measurements. Infinitely many signals explain the data, but you know the true one is sparse: most of its entries are zero. You can encode this knowledge with a regularizer that counts nonzero entries (the ℓ0 “norm”), or with its convex cousin, the ℓ1 norm. The resulting objective is not smooth: it has kinks (ℓ1) or even jumps (ℓ0), so the classical tools of smooth optimization, which rely on derivatives, do not apply directly.

Plots of x squared over two, absolute value of x, and the l0 norm
Three functions of one variable. A smooth function has a tangent everywhere; the ℓ1 norm has a kink at zero; the ℓ0 “norm” jumps at zero. R2N handles all three kinds of regularizers.
Sparse signal recovery with and without a regularizer
A toy sparse-recovery problem: 10 spikes among 200 entries, observed through 80 random measurements. Without a regularizer, the reconstruction is spread over all entries; with an ℓ0 regularizer, the spikes are recovered. (Illustration computed for this page with a simple proximal-gradient method.)

What R2N does differently

Proximal-gradient methods handle \(f + h\) with a gradient step on \(f\) followed by a proximal step on \(h\), which pulls the point toward the desired structure. They are simple, but they only use first-order information and can be slow.

Trust-region methods also use curvature and are usually more efficient, but in the nonsmooth setting each step requires the proximal operator of \(h\) combined with the trust-region constraint. For important regularizers, such as the rank, this combined operator is not known in closed form.

R2N replaces the trust-region constraint by an adaptive quadratic penalty \(\tfrac{\sigma}{2}\|s\|^2\) on the step. At each iteration it minimizes the sum of a quasi-Newton model of \(f\), a model of \(h\) and this penalty. If the step works well, \(\sigma\) decreases and longer steps are allowed; otherwise \(\sigma\) increases. The subproblem can be solved by a variant of the proximal-gradient method that only needs the proximal operator of \(h\) itself. The name reflects that R2N extends R2, a first-order quadratic-regularization method, with Newton-like curvature information.

We also introduce R2DH, a variant with a diagonal curvature model: when \(h\) is separable (ℓ0 or ℓ1, for example), each step is then computed in closed form, without any inner solver. R2DH can be used on its own or to solve the subproblems of R2N. Both methods also come in non-monotone versions.

What we proved

  • Global convergence of a first-order stationarity measure to zero, without assuming that the gradient of \(f\) is Lipschitz continuous, and while allowing the quasi-Newton models to grow unbounded, which is what happens in practice.
  • Tight worst-case complexity. When \(\nabla f\) is Lipschitz and the model Hessians grow at most like a power \(p\) of the number of successful iterations, R2N needs \(O(\varepsilon^{-2/(1-p)})\) iterations to bring the stationarity measure below ε when \(0 \le p < 1\), and \(\exp(c\,\varepsilon^{-2})\) when \(p = 1\). These bounds cannot be improved. They extend to nonsmooth problems what we established for trust-region methods with unbounded Hessian approximations.

In practice

R2N and R2DH are implemented in Julia in RegularizedOptimization.jl. We tested them on basis pursuit (sparse recovery), image denoising, minimum-rank matrix completion and a nonlinear support vector machine. The minimum-rank problem, in particular, cannot currently be solved directly by a trust-region approach, because the proximal operators it would require are not known analytically.

This project started during my 2023 research internship at GERAD, and I presented it at ISMP 2024 in Montréal.

For specialists

Problem: \(\min_x f(x) + h(x)\), with \(f \in C^1\) and \(h\) proper, lower semicontinuous and prox-bounded; both may be nonconvex. At iteration \(k\), R2N computes a step by approximately minimizing

\[ \nabla f(x_k)^T s + \tfrac12 s^T B_k s + \psi(s; x_k) + \tfrac12 \sigma_k \|s\|^2, \]

where \(B_k\) is a (possibly quasi-Newton) model Hessian, \(\psi(\cdot\,; x_k)\) models \(h(x_k + \cdot)\) and \(\sigma_k > 0\) is updated from the ratio of actual to predicted decrease. Unlike trust-region methods, the subproblem involves no indicator of a trust region, so a proximal-gradient inner solver only needs \(\mathrm{prox}_{\nu\psi}\).

Complexity: if \(\nabla f\) is Lipschitz and \(\|B_k\| = O(|S_k|^p)\), where \(S_k\) is the set of successful iterations up to \(k\), the stationarity measure falls below ε after \(O(\varepsilon^{-2/(1-p)})\) iterations for \(0 \le p < 1\), and after \(\exp(c\,\varepsilon^{-2})\) iterations for \(p = 1\); both bounds are tight.