Hessian Riemannian Langevin Monte Carlo (HRLMC)

• HRLMC is a geometric sampling algorithm that extends Langevin Monte Carlo using Hessian manifolds and mirror descent for high-dimensional, non-Euclidean problems.
• The algorithm provides non-asymptotic convergence bounds in Wasserstein and KL divergences by accommodating local metric properties and constrained geometries.
• By unifying Riemannian SDEs, interior-point methods, and optimal transport, HRLMC offers robust performance and explicit metric-dependent contraction rates.

The Hessian Riemannian Langevin Monte Carlo (HRLMC) algorithm is a geometric extension of Langevin Monte Carlo for efficient sampling from high-dimensional probability densities, particularly those lacking global regularity in the standard Euclidean setting. HRLMC leverages a Hessian manifold structure and mirror-descent techniques to model non-flat or constrained geometries, providing non-asymptotic convergence bounds in Wasserstein and Kullback–Leibler distances. This framework unifies Riemannian stochastic differential equations (SDEs), interior-point methods, and optimal transport theory, resolving difficulties in the analysis and computation of Riemannian discretizations by introducing relative smoothness, self-concordance, and explicit metric-dependent contraction rates (Zhang et al., 2020, Gatmiry et al., 2022).

1. Hessian Manifold Structure and Riemannian Metrics

Let $X \subset \mathbb{R}^p$ be open and contractible, and let $\varphi : X \to \mathbb{R}$ be a $C^2$ Legendre-type convex function (often called the "mirror potential" or barrier function). The Hessian manifold is endowed with the metric

$g_x(u, v) = u^\top [D^2 \varphi(x)]\, v$

where $D^2 \varphi(x)$ is positive definite for all $x \in X$. This construction allows dual coordinates $y = \nabla \varphi(x)$, in which the metric becomes Euclidean and Riemannian quantities push forward via the mirror map. Important classes of $\varphi$ include the Euclidean case $\varphi(x) = \frac{1}{2} \|x\|^2$ and the log-barrier $\varphi(x) = -\sum_i \log(a_i^\top x - b_i)$ for polytope constraints, which yield $g(x) = \nabla^2 \varphi(x)$ (Zhang et al., 2020, Gatmiry et al., 2022).

2. Riemannian Langevin Diffusion: SDE Formulation

Let $\pi$ be a target density $d\pi(x) = e^{-f(x)} dx$, $f \in C^3(X)$. The Riemannian Langevin diffusion (RLD) on $(X, g)$ is the Itô SDE:

$$dX_t = \Big[\theta(X_t) - D^2\varphi(X_t)^{-1} \nabla f(X_t)\Big]dt + \sqrt{2}\; D^2\varphi(X_t)^{-1/2} dB_t$$

where

$$\theta(x) = -D^2\varphi(x)^{-1} \, \mathrm{Tr}[ D^3 \varphi(x) D^2\varphi(x)^{-1} ]$$

The invariant measure of this diffusion is $\pi$ by Fokker–Planck analysis. In coordinates $y = \nabla\varphi(x)$ the SDE simplifies, and the diffusion is natural Brownian motion on the manifold, with drift $g^{-1}\nabla (\log \nu)$ and diffusion $g^{-1/2}$ (Zhang et al., 2020, Gatmiry et al., 2022).

3. Discrete HRLMC Scheme and Mirror Descent Analogy

The discretization leading to HRLMC is performed in the dual space and mapped back via the Legendre conjugate. With step-size $h_{k+1}$ and i.i.d. $\xi_{k+1} \sim \mathcal{N}(0, I_p)$: $$\begin{align*} \text{Initialize } x_0 \in X \ y_k = \nabla \varphi(x_k) \ y_{k+1} = y_k - h_{k+1} \nabla f(x_k) + \sqrt{2 h_{k+1} D^2\varphi(x_k)}\, \xi_{k+1} \ x_{k+1} = \nabla \varphi^*(y_{k+1}) \end{align*}$$ HRLMC is the sampling analogue of mirror descent: setting $\xi = 0$ yields the mirror descent step for $f$ with mirror map $\varphi$. The stochastic perturbation ensures the correct stationary distribution $\pi$. In the Euclidean case, HRLMC recovers Unadjusted Langevin Algorithm (ULA) (Zhang et al., 2020). For constrained domains (e.g., polytopes), ULA in local charts uses the metric and steps: $$x_{k+1} = x_k - \eta\; g(x_k)^{-1} \nabla f(x_k) + \sqrt{2\eta}\; g(x_k)^{-1/2} \xi_k$$ (Gatmiry et al., 2022).

4. Assumptions, Metric Properties, and Self-Concordance

Non-asymptotic convergence proofs require the following conditions:

• A1 (Self-concordance-like control): $\exists \kappa \geq 0$ such that $$\|D^2\varphi(x)^{1/2} - D^2\varphi(x')^{1/2}\|_F \leq \kappa \|\nabla\varphi(x) - \nabla\varphi(x')\|_2$$.
• A2 (Moment bound): $$R = \mathbb{E}_{X \sim \pi}[\|D^2\varphi(X)\|_2] < \infty$$.
• A3 (Relative strong convexity of $f$): $\exists m \geq 0$ with $$m\|\nabla\varphi(x) - \nabla\varphi(x')\|_2^2 \leq \langle\nabla f(x)-\nabla f(x'), \nabla\varphi(x)-\nabla\varphi(x')\rangle$$.
• A4 (Relative gradient Lipschitzness): $\exists M > 0$ such that $$\|\nabla f(x) - \nabla f(x')\|_2 \leq M \|\nabla\varphi(x) - \nabla\varphi(x')\|_2$$.
• A5 (Commutator bound): $$\|\left[D^{2}\varphi(x)^{-1}, D^2f(x)\right]\|_2 \leq \delta$$.

Self-concordance is extended to second-order conditions: bounds on derivatives of the metric and its inverse, controlling the Jacobian growth and ensuring discretization errors remain manageable (Gatmiry et al., 2022).

5. Non-Asymptotic Convergence Bounds and Wasserstein Geometry

HRLMC delivers explicit, non-asymptotic contraction rates in the Wasserstein ground metric: $$W_{2,\varphi}^2(\mu, \nu) = \inf_{X \sim \mu, X' \sim \nu} \mathbb{E} \left[\|\nabla\varphi(X) - \nabla\varphi(X')\|_2^2\right]$$ With step-size $h_{k+1}$ and contraction factor $\rho_{k+1}$, for all $k$,

$$W_{2,\varphi}(\mu_{k+1}, \pi) \leq \rho_{k+1}\, W_{2,\varphi}(\mu_k, \pi) + h_{k+1}\, p^{1/2}\, \beta_1(R,\kappa) + h_{k+1}^{3/2} p^{1/2} \beta_2(R,M,\kappa)$$

where $$\tilde{\kappa} = \sqrt{\kappa^2 + \delta (4M+\delta)/2(m+M)}$$ and explicit constants $\beta_1, \beta_2$ are supplied as functions of dimension, moment bounds, and regularity parameters. For constant step-size, iterates contract to an explicit Wasserstein ball around $\pi$ and remain inside, with radius proportional to $\kappa$ and local geometry (Zhang et al., 2020).

For second-order self-concordant metrics, the convergence in Kullback–Leibler divergence is: $$H_\nu(\rho_k) \leq e^{-c\,\alpha\,\epsilon\, k} H_\nu(\rho_0) + O(n^{5/2} L_2 + \sqrt{n} L_2^2 + n L_3) \frac{\epsilon}{\alpha}$$ where $H_\nu$ denotes KL divergence relative to $\nu$ and $L_2$, $L_3$ are Lipschitz constants of $F$ and its Hessian (Gatmiry et al., 2022).

6. Specializations and Applications

Euclidean LMC: Setting $\varphi(x) = \frac{1}{2}\|x\|^2$ yields $g = I$, $\kappa = \delta = 0$, recovers classical LMC with optimal Wasserstein bounds.

Constrained Sampling: For polytopes $\mathcal{M} = \{x: Ax \geq b\}$, using the log-barrier $\varphi(x) = -\sum_i \log(a_i^\top x - b_i)$ yields $g(x) = \nabla^2 \varphi(x)$ with self-concordant constants $\gamma_1 = 2, \gamma_2 = 4, \gamma_3 = 6$. HRLMC enables polynomial-time sampling of isoperimetric densities given gradient oracles and avoids global smoothness requirements.

Other Metrics: Adapting $\varphi$ to the geometry of $f$ allows sampling from densities in non-Euclidean geometries, e.g., Dirichlet or gamma posteriors on the simplex (Zhang et al., 2020).

7. Key Analytical Tools and Proofs

HRLMC analysis relies on:

• Extended Baillon–Haddad Inequality: For relative smoothness, establishes cocoercivity bounds between gradients of $f$ and the mirror map.
• Itô’s Isometry and Minkowski Integral Inequality: Manage stochastic error terms and bias in diffusion discretization.
• Chung’s Lemma: Demonstrates contraction to a Wasserstein ball under vanishing step-sizes.
• Bismut–Elworthy–Li Formula: Controls higher-order terms in stochastic flows.
• Self-Concordance Theory: Manages metric derivatives and stability (Nesterov–Nemirovskii).

These tools together establish dimension-dependent rates, validate non-asymptotic guarantees, and ensure robustness to variance in local geometry (Zhang et al., 2020, Gatmiry et al., 2022).

---

HRLMC thus constitutes an advanced geometric Monte Carlo method tailored for high-dimensional, non-Euclidean, and constrained sampling problems, bridging stochastic analysis, convex optimization, and differential geometry. The rigorous contraction rates and metric-dependent bounds enable its application to non-smooth log-concave target distributions and provide theoretical foundations for sampling in complex domains.