🔢 mathematics

Hadamard Langevin dynamics for sampling the l1-prior

This paper proposes Hadamard Langevin dynamics (HLD), a novel sampling method that utilizes a smooth but nonconvex overparameterization of the l1-norm to exactly recover non-smooth posteriors, and establishes a rigorous theoretical foundation proving the existence, uniqueness, geometric ergodicity, and convergence of this approach.

Original authors: Ivan Cheltsov, Federico Cornalba, Clarice Poon, Tony Shardlow

Published 2026-05-05

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Ivan Cheltsov, Federico Cornalba, Clarice Poon, Tony Shardlow

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

The Big Problem: The "Bumpy" Mountain

Imagine you are trying to find the best spot to set up a campsite in a vast, foggy wilderness. You have a map (your data) and a set of rules (your prior) that tell you where you should be. In the world of statistics, this is called sampling from a posterior distribution.

Usually, the "rules" of the game are smooth, like rolling down a gentle hill. But in many modern problems (like finding the most important features in a huge dataset), the rules include a sparsity requirement. This means you want your solution to have as many zeros as possible (ignoring the noise).

Mathematically, this sparsity rule creates a sharp, jagged peak (like a pyramid or a tent) rather than a smooth hill. This is the $\ell_1$ -prior.

The problem is that the standard tools for navigating this terrain (called Langevin dynamics) are like hikers who need smooth paths. If they try to walk up a jagged pyramid, they get stuck or slip off. To fix this, previous methods tried to "sand down" the sharp corners to make them smooth. But this is like smoothing a pyramid into a cone: you change the shape of the mountain, so you end up at the wrong campsite. You get a biased answer.

The New Solution: The "Hadamard" Trick

The authors of this paper propose a clever new way to navigate the jagged mountain without smoothing it out. They use a mathematical trick called Hadamard parameterization.

The Analogy: The Shadow Puppet
Imagine the jagged mountain (the $\ell_1$ -norm) is a shadow puppet.

Old way: You try to walk on the shadow itself. It's sharp and hard to traverse.
The Paper's way: Instead of walking on the shadow, you look at the hands making the shadow.

The authors realize that the jagged shape of the "sparsity mountain" can be created by multiplying two smooth, round variables together (let's call them $u$ and $v$ ).

If you multiply $u$ and $v$ , you get the jagged shape you need.
But $u$ and $v$ themselves are smooth and round (like a gentle hill).

By moving the hikers from the jagged mountain to the smooth "hand" variables ( $u$ and $v$ ), they can use standard, easy-to-use hiking tools (Langevin dynamics) to explore the terrain. Once they find the best spot in the smooth world, they multiply $u$ and $v$ back together to get the correct answer on the jagged mountain.

The Magic: Because this transformation is exact, the campsite they find in the smooth world corresponds perfectly to the campsite on the jagged mountain. There is no "smoothing bias." They get the true answer.

What the Paper Proves (The "Safety Check")

Just because you found a new path doesn't mean it's safe or reliable. The authors spent most of the paper proving that this new path works mathematically.

It Won't Break (Well-Posedness): They proved that the hikers won't get stuck, fall off a cliff, or wander off to infinity. The path is stable.
It Will Get You There (Ergodicity): They proved that if you walk long enough, you will eventually visit every part of the mountain in the correct proportion. You won't get stuck in a loop or miss a section.
The Digital Map Works (Discretization): Since computers can't take infinite tiny steps, the authors also proved that their "digital map" (the computer algorithm) gets closer and closer to the true path the smaller the steps are. They showed the computer simulation converges to the right answer.

Why This Matters (In Simple Terms)

No More "Sandpaper": Previous methods had to sand down the sharp corners of the problem, which changed the answer. This method keeps the sharp corners but changes the perspective so the sharpness doesn't matter.
Exactness: It recovers the exact desired distribution, not an approximation.
New Perspective: It treats a difficult, non-smooth problem as a smooth, over-parameterized one (using more variables than strictly necessary, but in a way that helps).

What They Didn't Do

The paper is a theoretical foundation.

They did not apply this to specific medical diagnoses or stock market predictions.
They did not claim it is faster than all other methods in every single scenario (though their experiments suggest it mixes well).
They focused strictly on the $\ell_1$ -prior (the standard sparsity rule). While they mention it could work for other similar rules, the paper's main proof is specifically for this one case.

Summary

Think of the paper as inventing a new pair of glasses. When you look at a jagged, impossible-to-climb mountain (the sparsity problem) through normal glasses, you can't climb it. Through these new "Hadamard" glasses, the mountain looks like a smooth, climbable hill. The authors proved that if you climb the hill and translate your position back, you are standing exactly where you were supposed to be on the jagged mountain, and they provided the safety manual to ensure the climb is mathematically sound.

Technical Summary: Hadamard–Langevin Dynamics for Sampling the $\ell_1$ -Prior

Problem Statement

The paper addresses the challenge of sampling from high-dimensional Bayesian posterior distributions involving non-smooth log-densities, specifically the $\ell_1$ -prior (Lasso). The target distribution is given by:
$\rho(x) = \frac{1}{Z} \exp\left(-\beta (\lambda \|x\|_1 + G(x))\right),$
where $G(x)$ represents the data likelihood and $\|x\|_1$ is the non-smooth $\ell_1$ -norm promoting sparsity.

Existing Langevin-based sampling methods (e.g., Unadjusted Langevin Algorithm, ULA) struggle with the non-differentiability of the $\ell_1$ term. Standard approaches typically rely on:

Smoothing techniques: Using the Moreau–Yosida envelope to approximate the non-smooth term, which introduces a bias in the target distribution that vanishes only as the smoothing parameter $\gamma \to 0$ .
Proximal mappings: Incorporating proximal operators (e.g., Proximal Gradient Langevin Dynamics), which can result in stationary distributions with point masses (e.g., at zero) mixed with continuous densities, or require inverting high-dimensional matrices in hierarchical formulations.

The authors seek a method that avoids smoothing bias and proximal operators while maintaining a tractable diffusion process that exactly recovers the target posterior $\rho(x)$ .

Methodology: Hadamard Overparameterization

The core innovation is the Hadamard–Langevin Dynamics (HLD), which utilizes a smooth, overparameterized representation of the $\ell_1$ -norm.

1. The Hadamard Representation

The authors exploit the identity:
$\|x\|_1 = \sum_i |x_i| = \min_{u \odot v = x} \left\{ \frac{1}{2}\|u\|_2^2 + \frac{1}{2}\|v\|_2^2 \right\},$
where $\odot$ denotes the element-wise (Hadamard) product. By lifting the variable $x \in \mathbb{R}^d$ to a pair of variables $(u, v) \in \mathbb{R}^d_+ \times \mathbb{R}^d$ such that $x = u \odot v$ , the non-smooth $\ell_1$ term is replaced by a smooth (though non-convex) quadratic term.

2. The Joint Distribution

Starting from the hierarchical Bayesian formulation of the Laplace prior (Gaussian mixture representation), the authors derive a joint distribution $\pi(u, v)$ over the lifted space:
$\pi(u, v) \propto \left( \prod_{i=1}^d u_i \right) \exp\left( -\beta \left( \frac{\lambda}{2}\|(u, v)\|_2^2 + G(u \odot v) \right) \right).$
Crucially, the marginal distribution of $x = u \odot v$ under $\pi$ is exactly the target $\rho$ , without any approximation bias.

3. The Dynamics

The authors define the Langevin dynamics on the lifted space $(u, v)$ to sample from $\pi$ . The resulting Stochastic Differential Equation (SDE) is:
$\begin{aligned} du_t &= -\lambda u_t dt - v_t \odot \nabla G(u_t \odot v_t) dt + \frac{1}{\beta u_t} dt + \sqrt{\frac{2}{\beta}} dW^1_t, \\ dv_t &= -\lambda v_t dt - u_t \odot \nabla G(u_t \odot v_t) dt + \sqrt{\frac{2}{\beta}} dW^2_t, \end{aligned}$
where $W^1_t$ and $W^2_t$ are independent Brownian motions. The term $1/(\beta u_t)$ arises from the log-density of the prior on $u$ . The dynamics are well-defined on the domain $X = \mathbb{R}^d_+ \times \mathbb{R}^d$ .

Key Contributions and Theoretical Results

1. Well-Posedness of Continuous Dynamics

The paper establishes a rigorous well-posedness theory for the continuous HLD system (Eq. 7):

Existence and Uniqueness: The authors prove the existence of a unique strong solution $(u_t, v_t)$ such that $u_t > 0$ almost surely for all $t > 0$ .
Proof Techniques: Two distinct approaches are provided:
1. Cylindrical Coordinates: Transforming the system into Cartesian coordinates $(y, z, v)$ where $u = \sqrt{y^2 + z^2}$ to remove singularities in the drift, allowing the application of standard SDE theory (Stroock and Varadhan).
2. Girsanov Transformation: Linking the system to the Cox–Ingersoll–Ross (CIR) and Ornstein–Uhlenbeck (OU) processes. This approach confirms that the process inherits the positivity property of the CIR process and establishes the Strong Feller property.
Ergodicity: The continuous dynamics are shown to be geometrically ergodic, converging to the invariant measure $\pi$ at an exponential rate, even when the data term $G$ is unbounded (e.g., quadratic loss), provided $x^\top \nabla G(x)$ is bounded below.

2. Numerical Scheme and Convergence

The authors propose an implicit-explicit time-stepping scheme (Eq. 16) to discretize the HLD. The scheme is implicit in the regularization terms (ensuring positivity of $u$ ) and explicit in the data term $G$ .

Strong Convergence: The paper proves that the numerical scheme converges strongly to the continuous solution as the step size $\Delta t \to 0$ . Specifically, for a fixed time interval $T$ , the error is bounded by $O(\Delta t^{1/2 - \epsilon})$ .
Overcoming Technical Challenges: The analysis addresses two main difficulties:
1. Inverse Moments: The drift contains $1/u$ terms. The authors derive inverse-moment bounds of order up to (but excluding) 2, which is necessary because the standard CIR condition $2\kappa\mu > \theta^2$ holds with equality in this setting.
2. Non-Global Lipschitzness: The drift is neither convex nor globally Lipschitz. The authors employ a coupling technique within large balls and a rapidly decaying bound for the escape probability to prove geometric ergodicity of the discrete scheme.
Invariant Measure Convergence: It is proven that the stationary distribution $\pi_{\Delta t}$ of the discretized scheme converges to the true invariant measure $\pi$ in the Wasserstein-1 distance as $\Delta t \to 0$ .

3. Relation to Existing Methods

Riemannian Connection: The HLD is identified as a specific instance of Riemannian Langevin dynamics on a manifold with a metric derived from the Hessian of a convex function, linking the method to Mirror-Langevin dynamics.
Comparison to Proximal Methods: Unlike Proximal Gradient Langevin Dynamics (PGLD), which can introduce point masses at zero, HLD produces a continuous density. Unlike Moreau–Yosida regularized methods (MYULA), HLD introduces no smoothing bias.

Results and Significance

The paper provides the first theoretical foundation for sampling from nonconvex, nonsmooth posteriors using overparameterized Langevin dynamics.

Theoretical Rigor: The work moves beyond heuristic applications of overparameterization in optimization to provide a complete probabilistic analysis of the resulting diffusion process, including existence, uniqueness, and convergence rates.
Exact Sampling: The method recovers the exact target posterior $\rho(x)$ without the bias inherent in smoothing approximations or the structural artifacts (point masses) of proximal methods.
Numerical Evidence: While the theoretical convergence rates for the bias of the discretized invariant measure are not quantified (due to singularities preventing standard Taylor expansion arguments), numerical experiments on a 1D problem suggest the bias is $O(\Delta t)$ . Experiments show HLD achieves faster mixing (higher Effective Sample Size) compared to PGLD and MYULA, with no observed bias in the histograms.

Limitations and Future Work

The authors modestly note that their current analysis does not provide quantitative convergence rates for the bias of the discretized invariant measure $\pi_{\Delta t}$ , as the singular $1/u$ terms prevent the direct application of standard semigroup expansion techniques. Future work aims to quantify this bias and extend the Gaussian mixture representation approach to other sparsity-inducing priors, such as group- $\ell_1$ (which is shown to be a trivial extension) and potentially $\ell_q$ priors for $q \in (0, 1)$ , though the latter currently lacks a closed-form Gaussian mixture representation.