# From stability of Langevin diffusion to convergence of proximal MCMC for non-log-concave sampling

Marien Renaud, Valentin de Bortoli, Arthur Leclaire, Nicolas Papadakis

► To cite this version:

Marien Renaud, Valentin de Bortoli, Arthur Leclaire, Nicolas Papadakis. From stability of Langevin diffusion to convergence of proximal MCMC for non-log-concave sampling. NeurIPS 2025 - 39th Annual Conference on Neural Information Processing Systems, Dec 2025, San Diego (California), United States. hal-05125391

HAL Id: hal-05125391

https://hal.science/hal-05125391v1

Submitted on 23 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# From stability of Langevin diffusion to convergence of proximal MCMC for non-log-concave sampling

Marien Renaud*

Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251
F-33400 Talence, France

Valentin de Bortoli

ENS, CNRS, PSL University
Paris, 75005, FRANCE

Arthur Leclaire

LTCI, Télécom Paris, IP Paris, France

Nicolas Papadakis

Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251
F-33400 Talence, France

## Abstract

We consider the problem of sampling distributions stemming from non-convex potentials with Unadjusted Langevin Algorithm (ULA). We prove the stability of the discrete-time ULA to drift approximations under the assumption that the potential is strongly convex at infinity. In many context, e.g. imaging inverse problems, potentials are non-convex and non-smooth. Proximal Stochastic Gradient Langevin Algorithm (PSGLA) is a popular algorithm to handle such potentials. It combines the forward-backward optimization algorithm with a ULA step. Our main stability result combined with properties of the Moreau envelope allows us to derive the first proof of convergence of the PSGLA for non-convex potentials. We empirically validate our methodology on synthetic data and in the context of imaging inverse problems. In particular, we observe that PSGLA exhibits faster convergence rates than Stochastic Gradient Langevin Algorithm for posterior sampling while preserving its restoration properties.

## 1 Introduction

Sampling a probability distribution $\pi \propto e^{-V}$ on $\mathbb{R}^d$ is a fundamental task in machine learning [2, 10, 13, 27, 41, 48, 87]. It becomes challenging when the potential $V$ is non-convex [38, 42]. A well-known approach is to solve the Langevin stochastic differential equation [56, 75], defined by

$$dx_t = -\nabla V(x_t)dt + \sqrt{2}dw_t, \quad (1)$$

where $(w_t)_{t \ge 0}$ is a Wiener process [50] on $\mathbb{R}^d$. This process can be seen as a continuous gradient descent on $V$ randomly perturbed by $\sqrt{2}dw_t$. Under mild assumptions, equation (1) admits $\pi$ as its invariant distribution [75]. Moreover, the solution converges geometrically in law to $\pi$ when $t$ goes to infinity if $V$ is strongly convex at infinity [29, 22]. In Section 4 we discuss alternative methods that have been developed to sample $\pi$, such as Gibbs Sampling [17] or Diffusion Models [83].

As the function $\nabla V = -\nabla \log \pi$ is often unknown [70], it is approximated by a drift function $b: \mathbb{R}^d \to \mathbb{R}^d$ such that $b \approx \nabla V$. In practice, the Euler-Maruyama discretization of equation (1) leads

<footnote>*corresponding author: marien.renaud@math.u-bordeaux.fr</footnote>
<table><caption>Table 1: Existing convergence results in the literature for PSGLA (3).</caption><thead><tr><th></th><th>f<br>non-convex</th><th>g<br>non-convex</th><th>Prox<sub>g</sub><br>inexact</th><th>W<sub>2</sub><br>metric</th><th>W<sub>p</sub><br>metric</th><th>TV<br>metric</th></tr></thead><tbody><tr><td>[28, 80, 81]</td><td>✗</td><td>✗</td><td>✗</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>[32]</td><td>✗</td><td>✗</td><td>✓</td><td>✓</td><td>✗</td><td>✗</td></tr><tr><td>This paper</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td><td>✓</td></tr></tbody></table>

to the inexact Unadjusted Langevin Algorithm (iULA)

$$X_{k+1} = X_k - \gamma b(X_k) + \sqrt{2\gamma} Z_{k+1}, \quad (2)$$

with a step-size $\gamma > 0$ and $Z_{k+1} \sim \mathcal{N}(0, I_d)$ i.i.d. It is called inexact because $b$ is an inexact approximation of $\nabla V$ and unadjusted because the invariant law of this Markov Chain is $\pi$ up to a discretization error (see Theorem 2 in Section 2). This Langevin Markov Chain could be adjusted by adding a Metropolis-Hastings accept-reject step [74, 66]. With ergodic MCMC, the uniform measure supported on the iterates $(X_k)_{0 \le k \le N}$ will be use as an estimator for $\pi$. To evaluate the quality of iULA, there is a need to quantify the shift between the invariant law of the MCMC $X_k$ and the target distribution $\pi$. This leads to investigate the stability of iULA invariant law to shift in the drift $b$. Existing stability results [94, 71] are limited by a discretization error. With a new proof strategy, we prove the first drift stability bound without any discretization error (Theorem 1 in Section 2).

In many applications the potential has a composite form $V = f + g$. For instance in Bayesian inverse problems [69, 15, 35], the goal is to sample the posterior distribution with $f$ the data-fidelity, which is usually smooth and possibly non-convex, and $g$ the regularization, which might be non-smooth and non-convex. Due to the non-smoothness of $V$, iULA that relies on evaluation of $\nabla V$ can not be implemented. A popular solution is to design a forward-backward Langevin scheme [9, 7, 81, 1], with the so-called Proximal Stochastic Gradient Langevin Algorithm (PSGLA) [28, 80] defined by

$$X_{k+1} = \text{Prox}_{\gamma g} \left( X_k - \gamma \nabla f(X_k) + \sqrt{2\gamma} Z_{k+1} \right), \quad (3)$$

where the proximal operator defined by $\text{Prox}_{\gamma g}(x) = \arg \min_{y \in \mathbb{R}^d} \frac{1}{2\gamma} \|x - y\|^2 + g(y)$ can be seen as an implicit gradient step. PSGLA has recently been studied [80, 32] in the convex setting. However, to our knowledge, there is no convergence analysis in the non-convex setting (see Table 1). We fill this gap by providing a convergence proof based on our stability result for Langevin Markov Chain in the case of $f$ smooth and $g$ non-smooth and weakly convex (Theorem 3 in Section 3).

## 1.1 Contributions

(1) We prove the first stability result for inexact Unadjusted Langevin Algorithm (iULA defined in equation (2)) with respect to the drift without discretization errors (Theorem 1 in Section 2). (2) We give the discretization error of iULA in various metrics including $p$-Wasserstein and Total Variation (Theorem 2 in Section 2). (3) Leveraging properties of the Moreau envelope, we demonstrate the first convergence proof of Proximal Stochastic Gradient Langevin Algorithm (PSGLA defined in equation (3)) (Theorem 3 in Section 3) and Inexact PSGLA (Theorem 4 in Section 3) for non-convex potentials. (4) We apply our theory to posterior sampling and show empirically that the mixing time of PSGLA is significantly shorter than the iULA (Section 5). In particular, we propose a Plug-and-Play version of PSGLA for image restoration that achieves state-of-the art restoration.

## 1.2 Notations and definitions

For a distribution $\mu$ on $\mathbb{R}^d$ and a function $\phi: \mathbb{R}^d \to \mathbb{R}$, we denote $\mu(\phi) = \int_{\mathbb{R}^d} \phi d\mu$. For two distributions $\mu, \nu$ on $\mathbb{R}^d$, the Total Variation (TV) writes

$$\|\mu - \nu\|_{TV} := \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ \|\phi\|_{\infty} \le 1}} \mu(\phi) - \nu(\phi). \quad (4)$$

For a measurable function $\phi: \mathbb{R}^d \to \mathbb{R}^d$, the push-forward of a probability distribution $\mu$ by $\phi$ denoted by $\phi\#\mu$ is the distribution defined for a measurable set $A \subset \mathbb{R}^d$ by $\phi\#\mu(A) = \mu(\phi^{-1}(A))$.
The p-Wasserstein distance between $\mu$ and $\nu$ is defined, for $p \ge 1$, by

$$W_p(\mu, \nu) = \left( \min_{\beta \in \Pi} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x - y\|^p d\beta(x, y) \right)^{\frac{1}{p}}, \quad (5)$$

with $\Pi$ the set of probability law $\beta$ on $\mathbb{R}^d \times \mathbb{R}^d$ such that $\phi_1\#\beta = \mu$ and $\phi_2\#\beta = \nu$, with $\phi_1 : (x, y) \in \mathbb{R}^d \times \mathbb{R}^d \mapsto x$ and $\phi_2 : (x, y) \in \mathbb{R}^d \times \mathbb{R}^d \mapsto y$.

A function $g$ is $\rho$-weakly convex, with $\rho > 0$, if and only if $g + \frac{\rho}{2}\|\cdot\|^2$ is convex. A function $g$ twice differentiable is $m$-strongly convex at infinity if and only if there exists $R \ge 0$ such that $\forall x \in \mathbb{R}^d$, if $\|x\| \ge R$, then $\nabla^2 g(x) \ge \mu I_d$. Note that we define the strong convexity at infinity only for twice differentiable functions, which is sufficient for this paper. This non-trivial notion of strong convexity on non-convex sets is discussed in Appendix E.

## 2 Sampling stability and discretization error

In this section, we study the sampling stability to drift errors of the inexact Unadjusted Langevin Algorithm (iULA) defined in equation (2). We introduce two iULA with two arbitrary drifts as

$$X_{k+1}^i = X_k^i - \gamma b^i(X_k^i) + \sqrt{2\gamma Z_{k+1}^i}, \quad (6)$$

for $i \in \{1, 2\}$, with the step-size $\gamma > 0$, $Z_{k+1}^i \sim \mathcal{N}(0, I_d)$ i.i.d. and the drifts $b^i \in \mathcal{C}^0(\mathbb{R}^d, \mathbb{R}^d)$. A Markov Chain $(X_k)_{k \ge 0}$ has an invariant law if there exists $\mu$ a probability distribution on $\mathbb{R}^d$ such that if $X_0 \sim \mu$, then $\forall k \ge 0$, $X_k \sim \mu$. A Markov Chain $(X_k)_{k \ge 0}$ is said to be geometrically ergodic, if the law of $X_k$ converges geometrically to the invariant law, i.e. there exist $A \ge 0$, $\rho \in (0, 1)$, an invariant law $\mu$, such that $W_1(p_{X_k}, \mu) \le A\rho^k$ with $p_{X_k}$ the law of $X_k$.

**Assumption 1.** There exist $L, R \ge 0$ and $m > 0$ such that the drift $b$ verifies

(i) $b$ is L-Lipschitz

(ii) $\forall x, y \in \mathbb{R}^d$ such that $\|x - y\| \ge R$, we have $\langle b(x) - b(y), x - y \rangle \ge m\|x - y\|^2$.

Assumption 1(i) imposes that the drift is sufficiently regular, without uncontrolled variations. In particular, the L-lipschitzness involves that the drift is L-weakly convex. Assumption 1(ii) is a relaxation of the strong convexity assumption, as it only holds for couple of points sufficiently far away. Assumption 1(ii) is called weak dissipativity [23, 59] or distant dissipativity [64]. Note that if the drift is $m$-strongly convex at infinity and L-weakly convex then it is weakly dissipative.

We define the $\ell_2(\mu)$ distance between two functions $f_1, f_2 : \mathbb{R}^d \to \mathbb{R}$ by

$$\|f_1 - f_2\|_{\ell_2(\mu)} = (\mathbb{E}_{Y \sim \mu} (\|f_1(Y) - f_2(Y)\|^2))^{\frac{1}{2}}. \quad (7)$$

**Theorem 1.** Let $b^1$ and $b^2$ satisfy Assumption 1. $X_k^1$ and $X_k^2$ that satisfy equation (6) are two geometrically ergodic Markov Chains with invariant laws $\pi_\gamma^1, \pi_\gamma^2$. Then for $\gamma_0 = \frac{m}{L^2}$ and $p \in \mathbb{N}^*$ there exist $C_p, C \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, we have

$$W_p(\pi_\gamma^1, \pi_\gamma^2) \le C_p \|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}^{\frac{1}{p}}, \quad (8)$$

$$\|\pi_\gamma^1 - \pi_\gamma^2\|_{TV} \le C \|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}. \quad (9)$$

The proof of Theorem 1 can be found in Appendix F.4. Theorem 1 shows that errors in the drifts do not accumulate, moreover if two drifts are close (in $\ell_2$ distance) then the corresponding sampled distributions are close. We also demonstrate a similar result in the V-norm, a generalization of the TV-distance (see Theorem 5 in Appendix F.3). More importantly, contrary to [71], our bound does not contain any discretization error. Theorem 1 does not only refine existing works, it also generalizes to $W_p$ the bounds of [71], which were $W_1(\pi_\gamma^1, \pi_\gamma^2) \le C_1\|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}^{\frac{1}{2}} + D_1\gamma^{\frac{1}{8}}$ for equation (8) in the case $p = 1$ and $\|\pi_\gamma^1 - \pi_\gamma^2\|_{TV} \le C\|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)} + D\gamma^{\frac{1}{4}}$ for equation (9).

Sketch of the Proof of Theorem 1 First, we prove that $X_k^1, X_k^2$ are geometrically ergodic. Then by using the triangle inequality and the geometric ergodicity we control the shift between the asymptotic laws by a shift between finite time laws. Finally this finite time shift is controlled by the shift between drifts by applying the Girsanov theory with the core Lemma 20 in Appendix F.2.
**Theorem 2.** Let $b = \nabla V$ verifies Assumption 1 and $\pi \propto e^{-V}$. Then, for $\gamma_0 = \frac{m}{L^2}$ and $\gamma \in (0, \gamma_0]$, the Markov Chain defined in (2) with the drift $b = \nabla V$ has an invariant law $\pi_\gamma$ and is geometrically ergodic. Moreover for $p \in \mathbb{N}^*$, there exist $D_p, D \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, we have

$$W_p(\pi_\gamma, \pi) \le D_p \gamma^{\frac{1}{2p}}, \quad \|\pi_\gamma - \pi\|_{TV} \le D \gamma^{\frac{1}{2}}.$$

Theorem 2 is proved in Appendix F.6 by following the same strategy than Theorem 1 with continuous-time Markov process instead of discrete-time Markov Chain. It provides a result on the discretization error of the Unadjusted Langevin Algorithm, which has been extensively studied in the literature for finite time [64], in the strongly convex setting [19, 8], in TV-distance [26] or KL-divergence [12]. To our knowledge, formulating this discretization error in $W_p$ for all $p \in \mathbb{N}$ and in TV under weak dissipativity (Assumption 1(ii)) is a novelty with respect to the literature.

We quantified the sampling shifts due to drift approximations in Theorem 1 and to the discretization error in Theorem 2 for iULA. We can now put together these results.

**Corollary 1.** Let $b$ and $\nabla V$ verify Assumption 1 and $\pi \propto e^{-V}$. Then for $\gamma_0 = \frac{m}{L^2}$, and $\gamma \in (0, \gamma_0]$, the Markov Chain defined in equation (2) with drift $b$ is geometrically ergodic with an invariant law $\hat{\pi}_\gamma$. Moreover for $p \in \mathbb{N}^*$, there exist $C_p, D_p, C, D \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, we have

$$\begin{align*}
W_p(\hat{\pi}_\gamma, \pi) &\le C_p \|b - \nabla V\|_{\ell_2(\hat{\pi}_\gamma)}^{\frac{1}{p}} + D_p \gamma^{\frac{1}{2p}} \\
\|\hat{\pi}_\gamma - \pi\|_{TV} &\le C \|b - \nabla V\|_{\ell_2(\hat{\pi}_\gamma)} + D \gamma^{\frac{1}{2}}.
\end{align*}$$

Corollary 1 is a direct consequence of Theorem 1 combined with Theorem 2. It quantifies the precision of sampling $\pi \propto e^{-V}$ with the practical algorithm (2) instead of the continuous process defined in equation (1). Corollary 1 extends previous known results established in $W_2$ in the strongly convex setting [8] to strongly convex at infinity setting with various metrics. The pseudo-metric $\| \cdot \|_{\ell_2(\hat{\pi}_\gamma)}$ appears to be the appropriate metric to quantify the difference between drifts.

## 3 Convergence of Proximal Stochastic Gradient Langevin Algorithm for non-convex potentials

We now apply the results of Section 2 to derive in Theorem 3 the first convergence result for PSGLA (3), summarized in Algorithm 1, in the non-convex setting specified by Assumption 2. Then, in Theorem 4 we deduce the convergence of Inexact PSGLA (3) when approximating the proximal operator of the non-convex potential.

**Assumption 2.** The potential $V$ is composite, i.e. $V = f + g$ where

(i) $f$ is $L_f$-smooth, i.e. $\nabla f$ is $L_f$-Lipschitz.

(ii) $g$ is $\rho$-weakly convex with $\rho > 0$, i.e. $g + \frac{\rho}{2} \| \cdot \|^2$ is convex.

In the case of posterior sampling, which will be detailed in Section 5, $f$ is typically the data-fidelity and $g$ the regularization. Under Assumption 2(i), $f$ can be non-convex, as it is the case in image despeckling [24] or black-hole interferometric imaging [85]. The weak convexity of $g$ ensures that the proximal operator $\text{Prox}_{\gamma g}$ is well defined as soon as $\gamma < \frac{1}{\rho}$. Assumption 2(ii) is verified for various convex regularizations as total-variation [79] or for some non-convex deep regularizers [45, 82, 39].

<caption>**Algorithm 1** Proximal Stochastic Gradient Langevin Algorithm (PSGLA)</caption>

1: **input:** $N > 0$, $y \in \mathbb{R}^m$, $X_0 \in \mathbb{R}^d$
2: **for** $k = 0, \dots, N-1$ **do**
3: &nbsp;&nbsp;&nbsp;&nbsp;$Z_{k+1} \sim \mathcal{N}(0, I_d)$
4: &nbsp;&nbsp;&nbsp;&nbsp;$X_{k+1} \leftarrow \text{Prox}_{\gamma g}(X_k - \gamma \nabla f(X_k) + \sqrt{2\gamma} Z_{k+1})$
5: **end for**
6: **output:** $(X_k)_{k \in [0, N]}$
**On the Moreau envelope** To analyse PSGLA under Assumption 2, we will need properties on the proximal operator and the Moreau envelope $g^\gamma$ of $g$ defined for $\gamma > 0$ and $x \in \mathbb{R}^d$ by $g^\gamma(x) = \inf_{y \in \mathbb{R}^d} \frac{1}{2\gamma} \|x - y\|^2 + g(y)$. Notice that if $g$ is $\rho$-weakly convex and $\rho\gamma < 1$, then the Moreau envelope is well defined. In the following Lemma 1 we provide important properties relative to the proximal operator and the Moreau envelope of weakly convex functions.

**Lemma 1.** For a $\rho$-weakly convex function $g$ and $\gamma\rho < 1$, we have the following properties:

(i) If $g$ is differentiable at $x$, then we have $\text{Prox}_{\gamma g}(x + \gamma \nabla g(x)) = x$.

(ii) The proximal operator is a gradient-descent step on the Moreau envelope, i.e. $\text{Prox}_{\gamma g}(x) = x - \gamma \nabla g^\gamma(x)$.

(iii) The proximal operator $\text{Prox}_{\gamma g}$ is $\frac{1}{1-\gamma\rho}$-Lipschitz.

(iv) The image of the proximal operator is almost convex, $\text{Leb}(\text{Conv}(\text{Prox}_{\gamma g}(\mathbb{R}^d)) \setminus \text{Prox}_{\gamma g}(\mathbb{R}^d)) = 0$, with $\text{Leb}$ the Lebesgue measure and $\text{Conv}(\text{Prox}_{\gamma g}(\mathbb{R}^d))$ is the convex envelop of $\text{Prox}_{\gamma g}(\mathbb{R}^d)$.

In this lemma, we adapt to the weakly convex case results known in the convex setting [62, 76, 44, 3, 40, 47]. In particular point (iv) requires a subtle adaptation with respect to the convex case. All proofs are postponed in Appendix D, where we provide a self-contained analysis of the Moreau envelope for weakly convex functions.

**ULA and PSGLA: two discretizations of the same continuous process** in the case of $g$ differentiable. First, the Euler-Maruyama discretization of equation (1) with a stepsize $\gamma > 0$ leads to the Unadjusted Langevin Algorithm (ULA) [26, 56], defined by

$$X_{k+1} = X_k - \gamma \nabla f(X_k) - \gamma \nabla g(X_k) + \sqrt{2\gamma} Z_{k+1}, \quad (10)$$

with $Z_{k+1} \sim \mathcal{N}(0, I_d)$ i.i.d. This Markov Chain has been shown to sample the target law $\pi \propto e^{-V}$ [56] for a non-convex potential $V$, up to a discretization error (see Theorem 2 in Section 2). ULA adds noise at each iteration to the iterate $X_k$ leading to noisy samples. Taking advantage of the composite nature of $V$, PSGLA (3) appears to be a good candidate for producing high quality samples. Indeed, for PSGLA the iterate $X_k$ is obtained after applying $\text{Prox}_{\gamma g}$, which, in practice, often corresponds to a denoising operator [88]. However PSGLA has not yet been shown to sample the law $\pi$ for a non-convex potential.

Secondly, a semi-implicit discretization of equation (1) with $\gamma \in (0, \frac{1}{\rho})$ and Lemma 1(i) gives

$$\begin{align*}
X_{k+1} &= X_k - \gamma \nabla f(X_k) - \gamma \nabla g(X_{k+1}) + \sqrt{2\gamma} Z_{k+1} \\
(I_d + \gamma \nabla g)(X_{k+1}) &= X_k - \gamma \nabla f(X_k) + \sqrt{2\gamma} Z_{k+1} \\
X_{k+1} &= \text{Prox}_{\gamma g} (X_k - \gamma \nabla f(X_k) + \sqrt{2\gamma} Z_{k+1}).
\end{align*}$$

Our postulate is as follows. As the ULA and PSGLA are two discretizations of the same continuous process, we expect that PSGLA also approximately samples the target distribution $\pi \propto e^{-V}$.

**Reformulation of PSGLA as a two-point algorithm** PSGLA (3) can be reformulated as

$$\begin{align*}
Y_{k+1} &= X_k - \gamma \nabla f(X_k) + \sqrt{2\gamma} Z_{k+1} \\
X_{k+1} &= \text{Prox}_{\gamma g}(Y_{k+1}).
\end{align*}$$

By Lemma 1(ii), we get that the iterates $Y_k$ verify the equation

$$\begin{align*}
Y_{k+1} &= \text{Prox}_{\gamma g}(Y_k) - \gamma \nabla f(\text{Prox}_{\gamma g}(Y_k)) + \sqrt{2\gamma} Z_{k+1} \\
&= Y_k - \gamma b^\gamma(Y_k) + \sqrt{2\gamma} Z_{k+1},
\end{align*}$$

where the drift $b^\gamma$ is defined for $y \in \mathbb{R}^d$ as

$$b^\gamma(y) = \nabla f(y - \gamma \nabla g^\gamma(y)) + \nabla g^\gamma(y). \quad (11)$$

Therefore, the shadow sequence $Y_k$ verifies a Markov Chain of the form of equation (2) with a drift that depends on the step-size $\gamma > 0$. In order to apply the result of Section 2, we need to make additional technical assumptions on $g$.
**Assumption 3.**

(i) $\forall \gamma \in (0, \frac{1}{\rho})$, $g$ is $L_g$-smooth on $\text{Prox}_{\gamma g}(\mathbb{R}^d)$.

(ii) $g^\gamma$ is $\mu$-strongly convex at infinity with $\mu \ge 8L_f + 4L_g$, i.e. there exists $\gamma_1 > 0$ and $R_0 \ge 0$ such that $\forall \gamma \in (0, \gamma_1]$, $\nabla^2 g^\gamma \ge \mu I_d$, on $\mathbb{R}^d \setminus B(0, R_0)$.

Assumption 3(i) is verified by non-smooth functions $g$ such as the characteristic function of a convex set. This assumption is not usual for $g$ that is generally not expected to be smooth on a subset of $\mathbb{R}^d$ [80, 28, 32]. Note that $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ can be a non-convex set if $g$ is not finite everywhere. However, Lemma 1(iv) shows that $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ only differs from being convex by a negligible set. Assumption 3(ii) is key to ensure that $Y_k$ is geometrically ergodic. Assumption 3(ii) could be removed by slightly modifying PSGLA as detailed in Appendix B. Lemma 16 in Appendix E shows that if $g$ is strongly convex at infinity and globally smooth, then Assumption 3(ii) is verified.

**Theorem 3.** Under Assumptions 2-3, there exist $r \in (0, 1)$, $C_1, C_2 \in \mathbb{R}_+$ such that $\forall \gamma \in (0, \bar{\gamma}]$, with $\bar{\gamma} = \min(\frac{1}{2L_g}, \frac{\mu}{32(L_f+L_g)^2}, \frac{1}{2\rho}, \gamma_1)$, where $L_g, L_f, \rho, \gamma_1$ are defined in Assumptions 2-3, and $\forall k \in \mathbb{N}$, we have

$$ \mathbf{W}_p(p_{Y_k}, \mu_\gamma) \le C_1 r^{k\gamma} + C_2 \gamma^{\frac{1}{2p}}, \quad (12) $$

with $p_{Y_k}$ the distribution of $Y_k$ and $\mu_\gamma \propto e^{-f-g^\gamma}$.

Moreover there exist $C_3, C_4 \in \mathbb{R}_+$ such that $\forall \gamma \in (0, \bar{\gamma}]$ and $\forall k \in \mathbb{N}$, we have

$$ \mathbf{W}_p(p_{X_k}, \nu_\gamma) \le C_3 r^{k\gamma} + C_4 \gamma^{\frac{1}{2p}}, \quad (13) $$

with $p_{X_k}$ the distribution of $X_k$ and $\nu_\gamma \propto \text{Prox}_{\gamma g} \# e^{-f-g^\gamma}$.

Theorem 3 is proved in Appendix G.2. It shows the convergence of PSGLA (3) in a non-convex setting. This convergence is exponential with an asymptotic bias due to discretization errors. Theorem 3 can also be proved in TV-distance following the same reasoning. Note that the distribution $\mu_\gamma$ and $\nu_\gamma$ are not exactly the target distribution $\pi$, the following Proposition 1 fills this gap by showing that $\mu_\gamma$ and $\nu_\gamma$ both converge to $\pi$ when $\gamma$ goes to zero.

Sketch of the proof of Theorem 3 First, we prove that under Assumptions 2-3, the drift $b^\gamma$ verifies Assumption 1 with constants independent of $\gamma$. So we can apply the results of Section 2. In particular, $(Y_k)_{k \ge 0}$ is geometrically ergodic, with invariant law $p_\infty$. Then we apply Theorem 1 to quantify the distance between $p_\infty$ and the invariant law of the MCMC with the drift $\nabla f + \nabla g^\gamma$, named $p_\gamma$. Finally, we apply Theorem 2 to quantify the distance between $p_\gamma$ and $\mu_\gamma \propto e^{-f-g^\gamma}$. The triangle inequality and Lemma 1(iii) conclude the proof.

**Proposition 1.** With $\pi \propto e^{-f-g}$, $\mu_\gamma \propto e^{-f-g^\gamma}$ and $\nu_\gamma = \text{Prox}_{\gamma g} \# \mu_\gamma$, for $p \ge 1$, we have

$$ \lim_{\gamma \to 0} \mathbf{W}_p(\mu_\gamma, \pi) = 0, \quad \lim_{\gamma \to 0} \mathbf{W}_p(\nu_\gamma, \pi) = 0. $$

Moreover, if $g$ is L-Lipschitz, there exists $E_p \in \mathbb{R}_+$ such that $\forall \gamma \in [0, \frac{2}{L^2}]$

$$ \mathbf{W}_p(\mu_\gamma, \pi) \le E_p(L_\gamma)^{\frac{1}{p}}, \quad \mathbf{W}_p(\nu_\gamma, \pi) \le E_p(L_\gamma)^{\frac{1}{p}} + L_\gamma. $$

Proposition 1 is proved in Appendix G.3. We choose to separate Proposition 1 from Theorem 3 as $g$ is not L-Lipschitz is many applications. Moreover, Proposition 1 is only based on the Moreau envelope approximation study, independently from the Markov Chain context of this paper.

**Inexact PSGLA** In practice, many functions $g$, such as the Total-Variation regularisation [79], have no closed-form $\text{Prox}_{\gamma g}$. Therefore, the proximal operator is approximated by an operator $S_\gamma \in C^0(\mathbb{R}^d, \mathbb{R}^d)$ such that $S_\gamma \approx \text{Prox}_{\gamma g}$. The authors of [32] studied the inexact PSGLA defined by

$$ \hat{X}_{k+1} = S_\gamma \left( \hat{X}_k - \gamma \nabla f(\hat{X}_k) + \sqrt{2\gamma} Z_{k+1} \right), \quad (14) $$

with $Z_{k+1} \sim N(0, I_d)$ i.i.d. We now extend the convergence result of [32, Theorem 3.14] obtained in a convex setting to the strong convexity at infinity assumption. To that end, we observe that the shadow sequence $\hat{Y}_k = \hat{X}_k - \gamma \nabla f(\hat{X}_k) + \sqrt{2\gamma} Z_{k+1}$ verifies equation (2) with the drift $\hat{b}^\gamma(y) = \nabla f(y - G_\gamma(y)) + G_\gamma(y)$, where $G_\gamma = \gamma^{-1}(I_d - S_\gamma)$.
<figure><img src="image_2.png" /><figcaption>Figure 1: Posterior Sampling with PnP-ULA (top row) and PSGLA (bottom row) with three different Gaussian Mixture prior distributions (in green) and three observations (in red) leading to three posterior distributions (in black, mode probabilities are represented by transparency). Algorithms are run with 10,000 steps. Quantitative metrics, discrete Wasserstein distance and $L_2$ distance between the estimated density and the true density, are computed for each algorithm output. The Wasserstein distance during the iterations (right column) is averaged between the three experiments. Note that SnoPnP-ULA succeeds to explore the posterior even if there are various modes.</figcaption></figure>

**Theorem 4.** Under Assumption 2-3, if $\hat{b}^\gamma$ verifies Assumption 1, there exist $\hat{\gamma} > 0$ and $C_5, C_6, C_7 \in \mathbb{R}_+$ such that $\forall \gamma \in (0, \hat{\gamma}]$, the shadow sequence $\hat{Y}_k$ has an invariant law $\hat{\mu}_\gamma$ and $\forall k \in \mathbb{N}$, the distribution $p_{\hat{X}_k}$ of $\hat{X}_k$ defined in equation (14) verifies

$$ \mathbf{W}_p(p_{\hat{X}_k}, \mu_\gamma) \le C_5 r^{k\gamma} + C_6 \gamma^{\frac{1}{2p}} + \frac{C_7}{\gamma^{\frac{1}{p}}} \|S_\gamma - \text{Prox}_{\gamma g}\|_{\ell_2(\hat{\mu}_\gamma)}^{\frac{1}{2p}}. $$

Theorem 4 is proved in Appendix G.4 by reformulating inexact PSGLA as a two-point algorithm and applying Theorem 3. Theorem 4 states that approximations in the proximal operator evaluation do not accumulate and imply a quantified shift in the PSGLA law $p_{\hat{X}_k}$.

# 4 Related Works

An alternative method to sample composite potential is Gibbs sampling that have been developed in [37, 84]. Gibbs sampler splits the sampling problem into two easier sampling sub-problems. More recently, the authors of [90, 17, 57] propose a Gibbs sampler that, combined with deep learning, provides accurate sampling. However their theoretical guarantees are only developed in the convex setting. A non-convex proof of a coarse Gibbs algorithm, without annealing, is proposed in [85].

For image applications, diffusion models [83] are another recent strategy to sample a law $\pi$ with convergence guarantees [11, 21]. Many works have focused on adapting diffusion models for image posterior sampling, see for instance [51, 15, 14, 78, 97]. However, to our knowledge, there is no drift stability results have been shown in this setting, therefore preventing the establishment of convergence results for posterior sampling. More details on related works are provided in Appendix A.

# 5 Application to posterior sampling

Recovering a signal $x \in \mathbb{R}^d$ from a degraded observation of this signal $y \in \mathbb{R}^m$, with typically $m < d$, is called inverse problem. In many settings, a degradation model is given, i.e. $y = A(x) + n$, with $n$ some noise. Even though our approach is not limited to this setting, for the sake of simplicity, in this section we will focus on the linear inverse problem $y = Ax + n$ with $A \in \mathbb{R}^{m \times d}$ and additive Gaussian noise $n \sim N(0, \sigma^2 I_d)$.

In a Bayesian paradigm, we suppose that $x$ is a realization of a random variable having the density $p(x)$, called prior distribution. Then, we aim to sample the distribution $p(x|y)$, called posterior
distribution. The Bayes formula gives that $p(x|y) \propto p(y|x)p(x)$ and the physical model implies that $p(y|x) \propto e^{-\frac{1}{\sigma^2}\|Ax-y\|^2}$. By denoting $f(x) = \frac{1}{\sigma^2}\|Ax-y\|^2$ and $g(x) = -\log p(x)$, we get that

$$p(x|y) \propto e^{-V(x)},$$

with $V = f + g$ a composite potential. Thus, we can use inexact ULA (2) (with $b = \nabla f + \nabla g$) or PSGLA (3) (with $\nabla f$ and $\text{Prox}_{\gamma g}$), to sample this posterior distribution.

**Plug-and-Play** Sampling the posterior distribution is a challenging task as it can be a high dimension problem (e.g. image inverse problems) as the prior distribution $p(x)$ is often unknown. Therefore, it is not possible to compute $\nabla g$ or $\text{Prox}_{\gamma g}$. The authors of [88] proposed to replace $\text{Prox}_{\gamma g}$ by a pretrained denoiser $D_{\sqrt{\gamma}}$, learnt to remove additive Gaussian noise of standard-deviation $\sqrt{\gamma} > 0$. In fact, $\text{Prox}_{\gamma g}(x) = \arg\min_{y \in \mathbb{R}^d} \frac{1}{2\gamma} \|x - y\|^2 + g(y)$ is the maximum-a-posteriori denoising with noise level $\sqrt{\gamma}$ and prior distribution $e^{-g}$. Replacing $\text{Prox}_{\gamma g}$ by a pretrained denoiser $D_{\sqrt{\gamma}}$ is called Plug-and-Play (PnP) and $D_{\sqrt{\gamma}}$ implicitly encodes the prior knowledge.

We thus introduce the PnP-PSGLA algorithm defined by

$$X_{k+1} = D_{\sqrt{\gamma}} \left( X_k - \frac{\gamma}{\lambda} \nabla f(X_k) - \sqrt{2\gamma} Z_{k+1} \right), \quad (15)$$

with $\lambda > 0$ a regularization parameter, $D_{\sqrt{\gamma}}$ a pre-trained denoiser of noise-level $\sqrt{\gamma}$ and $Z_{k+1} \sim N(0, I_d)$. This parameter $\lambda$ can be seen as a temperature parameters by replacing the prior $p$ by $p^\lambda$ which leads to sample a distribution proportional to $e^{-\frac{1}{\lambda}f-g}$ instead of $e^{-f-g}$. This parameter allows more flexibility in practice to improve PnP-PSGLA performance. Note that at each iteration, we inject $\sqrt{2}$ times more noise than we denoise; therefore, we can expect the algorithm to explore the space. If we injected as much noise as we denoise, we would expect the algorithm to behave like a stochastic optimization method rather than a sampling algorithm as detailed in [55, 73, 72].

# 6 Experiments

In this section, we present numerical validation of PnP-PSGLA (15) for posterior sampling. More details and additional experiments are provided in Appendix C.

**Gaussian Mixture in 2D** In order to test experimentally the quality of PnP-PSGLA, compared to PnP-ULA [56], we first look at a simplified case with a Gaussian Mixture prior in 2D. PnP-ULA is a Plug-and-Play version of the ULA (10). In this case, the posterior distribution, $\nabla g$ and the MMSE denoiser can be written in closed form. Therefore, we can run and compare the sampling algorithms PnP-PSGLA and PnP-ULA with the exact prior information.

As illustrated in Figure 1, PnP-PSGLA outperforms PnP-ULA both quantitatively and qualitatively on this Gaussian mixture problem. More precisely, PnP-ULA fails to discover minor modes of the posterior distribution in $N = 10^4$ steps, whereas PnP-PSGLA successfully identifies them. Our experiments suggest that the mixing time of PnP-PSGLA is shorter than for PnP-ULA. If the proximal operator is interpreted as a projection operator, PnP-PSGLA rigidly project the iterates on the "clean data" manifold, whereas PnP-ULA softly guides the iterates. This might help PnP-PSGLA to discover modes faster that PnP-ULA.

**Image restoration with Plug-and-Play approach** We use the Denoising Convolutional Neural Network (DnCNN) denoiser proposed by the authors of [67] and trained to be firmly non-expansive. This denoiser has been trained for one level of noise $\sqrt{\gamma} = 2/255$, which fixes the step-size $\gamma$ of PSGLA. We also considered the proximal deep neural network denoiser of [46], based on the DRUNet architecture [95]. Interestingly, we found that the choice of denoiser greatly influences the quality of the result. As DnCNN provided much better results than DRUNet we consistently use this architecture and report DRUNet restoration results in Appendix C.2.1.

In Figure 2, we present restorations with different methods and denoisers. DiffPIR [97] is a diffusion-inspired algorithm to restore images. It uses the denoiser with various noise levels, excluding the use of the DnCNN network [67] trained with only one level of noise. Plug-and-Play Unadjusted Langevin Algorithm (PnP-ULA) [56] is a Plug-and-Play adaptation of the Unadjusted Langevin Algorithm (10)
<figure><img src="image_3.png" /><figcaption>Figure 2: Qualitative result for image inpainting with 50% masked pixels and a noise level of $\sigma = 1/255$. PnP-ULA is run with 1,000,000 iterations and PSGLA with 10,000 iterations.</figcaption></figure>

<table><caption>Table 2: Quantitative results for image inpainting with 50% masked pixels on CBSD68 dataset. The number of inference of the denoiser *N* and the time of restoration in second on one GPU NVIDIA A100 Tensor Core are given for one image. Best and second best results are respectively displayed in bold and underlined. Note that PnP-PSGLA (15) has comparable performances with state-of-the art methods as DiffPIR and PnP and is significantly faster than PnP-ULA.</caption><thead><tr><th>Algorithm</th><th>Denoiser</th><th>PSNR ↑</th><th>SSIM ↑</th><th>LPIPS ↓</th><th>N ↓</th><th>time (s) ↓</th><th>convergent</th></tr></thead><tbody><tr><td>DiffPIR [97]</td><td>GSDRUNet</td><td>29.99</td><td>0.88</td><td><u>0.06</u></td><td><b>20</b></td><td><b>1</b></td><td><b>X</b></td></tr><tr><td>RED [77]</td><td>DnCNN</td><td>30.49</td><td>0.89</td><td><u>0.06</u></td><td><u>500</u></td><td><u>6</u></td><td>✓</td></tr><tr><td>RED [77]</td><td>GSDRUNet</td><td>29.26</td><td>0.88</td><td>0.12</td><td><u>500</u></td><td>20</td><td>✓</td></tr><tr><td>PnP [88]</td><td>DnCNN</td><td>30.50</td><td><u>0.91</u></td><td><u>0.06</u></td><td><u>500</u></td><td><u>6</u></td><td>✓</td></tr><tr><td>PnP [88]</td><td>GSDRUNet</td><td><u>30.52</u></td><td><b>0.92</b></td><td>0.07</td><td><u>500</u></td><td>20</td><td>✓</td></tr><tr><td>PnP-ULA [56]</td><td>DnCNN</td><td>27.89</td><td>0.82</td><td>0.12</td><td>100,000</td><td>1,200</td><td>✓</td></tr><tr><td>PnP-PSGLA [32]</td><td>TV</td><td>29.24</td><td>0.89</td><td>0.08</td><td>1,000</td><td>25</td><td>✓</td></tr><tr><td>PnP-PSGLA</td><td>DnCNN</td><td><b>30.81</b></td><td><b>0.92</b></td><td><b>0.05</b></td><td>10,000</td><td>120</td><td>✓</td></tr></tbody></table>

that is proved to sample the posterior distribution. PnP-PSGLA is run with two denoisers. The Total Variation (TV) denoiser [79] is a classical denoising method that leads to staircasing effects and over-smoothes the denoised image. We use TV-denoiser as a baseline with limited computational time. PnP-PSGLA is run with the DnCNN denoiser proposed by [67] to improve the restoration quality. Note that PnP-PSGLA, when run with DnCNN, achieves good perceptual quality using 100 times fewer steps than PnP-ULA. We again interpret this acceleration by the ability of PnP-PSGLA to rigidly project the iterates whereas PnP-ULA softly guides them.

On Table 2, we present quantitative comparisons between various restoration methods on the CBSD68 dataset [61]. Results are presented with two distortion metrics, Peak Signal to Noise Ratio (PSNR) and Structural Similarity (SSIM) [91], and one deep perceptual metric, Learned Perceptual Image Patch Similarity (LPIPS) [96]. Regularization by Denoising (RED) [77] and Plug-and-Play (PnP) [88] are restoration algorithms that solve an optimization problem to find the most probable image knowing the degraded observation. The Gradient-Step DRUNet (GSDRUNet) [45] denoiser, based on the DRUNet architecture [95], is a deep neural network denoiser that ensures convergence for RED and PnP algorithms. For PnP-ULA and PnP-PSGLA, the mean of the iterates is taken as an estimator of the expectation of the posterior distribution. Note that PnP-PSGLA with DnCNN achieves performances that are competitive with state-of-the-art methods. Moreover, the PnP-PSGLA achieves better performances that PnP-ULA in less iterations. PnP-PSGLA with DnCNN requires more computational resources than DiffPIR, PnP or RED. Using the TV-denoiser is an alternative to accelerate PnP-PSGLA at the expense of visual performance.

## 7 Conclusion

In this paper, we aim to sample a distribution $\pi \propto e^{-V}$ with non-convex potential $V$. We improve existing results for the stability (Theorem 1 in Section 2) and the discretization error (Theorem 2 in Section 2) of iULA (2). We use these results to derive the first convergence proof of PSGLA (3) for non-log-concave sampling (Theorem 3 in Section 3). We apply our theoretical findings to posterior sampling and show the experimental benefit, especially in computational time, of PSGLA with
respect to the original iULA. We also propose a Plug-and-Play version of PSGLA (PnP-PSGLA) that achieves state-of-the art image restoration.

PnP-PSGLA is limited by its computational time, as is the case with most Markov Chain Monte Carlo (MCMC) methods. Annealing methods [85] or more efficient denoisers could offer promising directions for future works to reduce the computational time of PnP-PSGLA. Convergent posterior sampling algorithms could have a positive societal impact by providing image restoration with quantified uncertainty. Also, such a deep image restoration algorithm could also have negative societal impacts as it may hallucinate details and might be used to generate false information.

## Acknowledgements

This study has been carried out with financial support from the French Direction Générale de l'Armement and the French National Research Agency through Projects PEPR PDE-AI, ANR-19-CE40-005 MISTIC and ANR-23-CE40-0017 SOCOT. Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, supported by Inria, CNRS (LABRI and IMB), Université de Bordeaux, Bordeaux INP and Conseil Régional d'Aquitaine (see https://www.plafrim.fr). We thank Andrés Almansa, Pascal Bianchi, Kaplan Desbouis, Gersende Fort, Erell Gachon, Julien Hermant, Samuel Hurault, Rémi Laumont, Éloan Rapion, Adrien Richou for their time and discussions.

## References

[1] Patrice Abry, Gersende Fort, Barbara Pascal, and Nelly Pustelnik. Credibility intervals for the reproduction number of the Covid-19 pandemic using proximal Langevin samplers. In *31th European Signal Processing Conference (EUSIPCO)*, 2023.

[2] Guillaume Alain and Yoshua Bengio. What regularized auto-encoders learn from the data-generating distribution. *The Journal of Machine Learning Research*, 15(1):3563–3593, 2014.

[3] Heinz H Bauschke, Patrick L Combettes, Heinz H Bauschke, and Patrick L Combettes. *Correc- tion to: convex analysis and monotone operator theory in Hilbert spaces*. Springer, 2017.

[4] Matej Benko, Iwona Chlebicka, Jorgen Endal, and Błażej Miasojedow. Langevin Monte Carlo beyond Lipschitz gradient continuity. In *Proceedings of the AAAI Conference on Artificial Intelligence*, volume 39, pages 15541–15549, 2025.

[5] Espen Bernton. Langevin Monte Carlo and JKO splitting. In *Conference on Learning Theory*, pages 1777–1798. PMLR, 2018.

[6] Radu Ioan Boț and Axel Böhm. Alternating proximal-gradient steps for (stochastic) nonconvex-concave minimax problems. *SIAM Journal on Optimization*, 33(3):1884–1913, 2023.

[7] Nicolas Brosse, Alain Durmus, Éric Moulines, and Marcelo Pereyra. Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo. In *Conference on Learning Theory*, pages 319–342. PMLR, 2017.

[8] Nicolas Brosse, Alain Durmus, Éric Moulines, and Sotirios Sabanis. The tamed unadjusted Langevin algorithm. *Stochastic Processes and their Applications*, 129(10):3638–3663, 2019.

[9] Sebastien Bubeck, Ronen Eldan, and Joseph Lehec. Finite-time analysis of projected Langevin Monte Carlo. *Advances in Neural Information Processing Systems*, 28, 2015.

[10] Olivier Cappé, Simon J Godsill, and Eric Moulines. An overview of existing methods and recent advances in sequential monte carlo. *Proceedings of the IEEE*, 95(5):899–924, 2007.

[11] Sitan Chen, Sinho Chewi, Jerry Li, Yuanzhi Li, Adil Salim, and Anru R Zhang. Sampling is as easy as learning the score: theory for diffusion models with minimal data assumptions. *arXiv preprint arXiv:2209.11215*, 2022.

[12] Xiang Cheng and Peter Bartlett. Convergence of Langevin MCMC in KL-divergence. In *Algorithmic learning theory*, pages 186–211. PMLR, 2018.
[13] Sinho Chewi, Murat A Erdogdu, Mufan Li, Ruoqi Shen, and Shunshi Zhang. Analysis of Langevin Monte Carlo from Poincare to Log-Sobolev. In Po-Ling Loh and Maxim Raginsky, editors, *Proceedings of Thirty Fifth Conference on Learning Theory*, volume 178 of *Proceedings of Machine Learning Research*, pages 1–2. PMLR, 02–05 Jul 2022.

[14] Hyungjin Chung and Jong Chul Ye. Score-based diffusion models for accelerated mri. *Medical image analysis*, 80:102479, 2022.

[15] Hyungjin Chung, Jeongsol Kim, Michael Thompson McCann, Marc Louis Klasky, and Jong Chul Ye. Diffusion posterior sampling for general noisy inverse problems. In *The Eleventh International Conference on Learning Representations*, 2022.

[16] Frank H Clarke. Nonsmooth analysis and optimization. In *Proceedings of the international congress of mathematicians*, volume 5, pages 847–853. Citeseer, 1983.

[17] Florentin Coeurdoux, Nicolas Dobigeon, and Pierre Chainais. Plug-and-play split gibbs sampler: embedding deep generative priors in bayesian inference. *IEEE Transactions on Image Processing*, 2024.

[18] Francesca R Crucinio, Alain Durmus, Pablo Jiménez, and Gareth O Roberts. Optimal scaling results for Moreau-Yosida Metropolis-adjusted Langevin algorithms. *Bernoulli*, 31(3):1889–1907, 2025.

[19] Arnak Dalalyan and Avetik Karagulyan. User-friendly guarantees for the langevin monte carlo with inaccurate gradient. *Stochastic Processes and their Applications*, 129(12):5278–5311, 2019.

[20] Damek Davis and Dmitriy Drusvyatskiy. Stochastic model-based minimization of weakly convex functions. *SIAM Journal on Optimization*, 29(1):207–239, 2019.

[21] Valentin De Bortoli. Convergence of denoising diffusion models under the manifold hypothesis. *Transactions on Machine Learning Research*, 2022.

[22] Valentin De Bortoli and Alain Durmus. Convergence of diffusions and their discretizations: from continuous to discrete processes and back, 2019. URL https://arxiv.org/abs/1904.09808.

[23] Arnaud Debussche, Ying Hu, and Gianmario Tessitore. Ergodic bsdes under weak dissipative assumptions. *Stochastic Processes and their applications*, 121(3):407–426, 2011.

[24] Charles-Alban Deledalle, Loic Denis, Sonia Tabti, and Florence Tupin. Mulog, or how to apply gaussian denoisers to multi-channel sar speckle reduction? *IEEE Transactions on Image Processing*, 26(9):4389–4403, 2017.

[25] Randal Douc, Eric Moulines, Pierre Priouret, Philippe Soulier, Randal Douc, Eric Moulines, Pierre Priouret, and Philippe Soulier. *Markov chains: Basic definitions*. Springer, 2018.

[26] Alain Durmus and Eric Moulines. Nonasymptotic convergence analysis for the unadjusted Langevin algorithm. *The Annals of Applied Probability*, 2017.

[27] Alain Durmus, Eric Moulines, and Marcelo Pereyra. Efficient bayesian computation by proximal markov chain monte carlo: when langevin meets moreau. *SIAM Journal on Imaging Sciences*, 11(1):473–506, 2018.

[28] Alain Durmus, Szymon Majewski, and Błażej Miasojedow. Analysis of Langevin Monte Carlo via convex optimization. *Journal of Machine Learning Research*, 20(73):1–46, 2019.

[29] Andreas Eberle. Reflection couplings and contraction rates for diffusions. *Probability theory and related fields*, 166:851–886, 2016.

[30] David A Edwards. On the kantorovich–rubinstein theorem. *Expositiones Mathematicae*, 29(4): 387–398, 2011.

[31] Bradley Efron. Tweedie's formula and selection bias. *Journal of the American Statistical Association*, 106(496):1602–1614, 2011.
[32] Matthias J Ehrhardt, Lorenz Kuger, and Carola-Bibiane Schönlieb. Proximal Langevin sampling with inexact proximal mapping. *SIAM Journal on Imaging Sciences*, 17(3):1729–1760, 2024.

[33] Paula Cordero Encinar, Francesca R Crucinio, and O Deniz Akyildiz. Proximal interacting particle langevin algorithms, 2024. URL https://arxiv.org/abs/2406.14292.

[34] Rémi Flamary, Nicolas Courty, Alexandre Gramfort, Mokhtar Z Alaya, Aurélie Boisbunon, Stanislas Chambon, Laetitia Chapel, Adrien Corenflos, Kilian Fatras, Nemo Fournier, et al. Pot: Python optimal transport. *Journal of Machine Learning Research*, 22(78):1–8, 2021.

[35] Gersende Fort, Barbara Pascal, Patrice Abry, and Nelly Pustelnik. Covid19 reproduction number: Credibility intervals by blockwise proximal monte carlo samplers. *IEEE Transactions on Signal Processing*, 71:888–900, 2023.

[36] Anne Gagneux, Mathurin Massias, Emmanuel Soubies, and Rémi Gribonval. Convexity in ReLU Neural Networks: beyond ICNNs?, 2025. URL https://arxiv.org/abs/2501.03017.

[37] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. *IEEE Transactions on pattern analysis and machine intelligence*, 1984.

[38] Wally R Gilks, Nicky G Best, and Keith KC Tan. Adaptive rejection metropolis sampling within gibbs sampling. *Journal of the Royal Statistical Society Series C: Applied Statistics*, 44(4): 455–472, 1995.

[39] Alexis Goujon, Sebastian Neumayer, and Michael Unser. Learning weakly convex regularizers for convergent image-reconstruction algorithms. *SIAM Journal on Imaging Sciences*, 17(1): 91–115, 2024.

[40] Rémi Gribonval and Mila Nikolova. A characterization of proximity operators. *Journal of Mathematical Imaging and Vision*, 62(6):773–789, 2020.

[41] Bichuan Guo, Yuxing Han, and Jiangtao Wen. Agem: Solving linear inverse problems via deep priors and sampling. *Advances in Neural Information Processing Systems*, 32, 2019.

[42] Wei Guo, Molei Tao, and Yongxin Chen. Provable benefit of annealed Langevin Monte Carlo for non-log-concave sampling. *International Conference of Learning Representation*, 2025.

[43] Jean-Baptiste Hiriart-Urruty and Claude Lemaréchal. Convex analysis and minimization algorithms I: Fundamentals, volume 305. Springer science & business media, 1996.

[44] Tim Hoheisel, Maxime Laborde, and Adam Oberman. On proximal point-type algorithms for weakly convex functions and their connection to the backward euler method. *Optimization Online*, 2010.

[45] Samuel Hurault, Arthur Leclaire, and Nicolas Papadakis. Gradient step denoiser for convergent plug-and-play. In *International Conference on Learning Representations (ICLR'22)*, 2022.

[46] Samuel Hurault, Antonin Chambolle, Arthur Leclaire, and Nicolas Papadakis. Convergent plug-and-play with proximal denoiser and unconstrained regularization parameter. *Journal of Mathematical Imaging and Vision*, pages 1–23, 2024.

[47] Valdines Leite de Sousa Junior, Lucas Vidal de Meireles, Samara Costa Lima, and Gilson do Nascimento Silva. Local convergence of an inexact proximal algorithm for weakly convex functions, 2023. URL https://arxiv.org/abs/2308.02014.

[48] Zahra Kadkhodaie and Eero Simoncelli. Stochastic solutions for linear inverse problems using the prior implicit in a denoiser. *Advances in Neural Information Processing Systems*, 34: 13242–13254, 2021.

[49] Leonid Vasilevich Kantorovich and SG Rubinshtein. On a space of totally additive functions. *Vestnik of the St. Petersburg University: Mathematics*, 13(7):52–59, 1958.

[50] Ioannis Karatzas and Steven Shreve. *Brownian motion and stochastic calculus*, volume 113. Springer Science & Business Media, 1991.
[51] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song. Denoising diffusion restoration models. *Advances in Neural Information Processing Systems*, 35:23593–23606, 2022.

[52] Teresa Klatzer, Paul Dobson, Yoann Altmann, Marcelo Pereyra, Jesús María Sanz-Serna, and Konstantinos C Zygalakis. Accelerated Bayesian imaging by relaxed proximal-point Langevin sampling. *SIAM Journal on Imaging Sciences*, 17(2):1078–1117, 2024.

[53] Andrew Lamperski. Projected stochastic gradient Langevin algorithms for constrained sampling and non-convex learning. In *Conference on Learning Theory*, pages 2891–2937. PMLR, 2021.

[54] Robert Lang. A note on the measurability of convex sets. *Archiv der Mathematik*, 47:90–92, 1986.

[55] Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra. On maximum a posteriori estimation with plug & play priors and stochastic gradient descent. *Journal of Mathematical Imaging and Vision*, 65(1):140–163, 2023.

[56] Rémi Laumont, Valentin De Bortoli, Andrés Almansa, Julie Delon, Alain Durmus, and Marcelo Pereyra. Bayesian imaging using Plug & Play Priors: When Langevin meets Tweedie. *SIAM Journal on Imaging Sciences*, 15(2):701–737, 2022.

[57] Jiaming Liang and Yongxin Chen. A proximal algorithm for sampling from non-smooth potentials. In *2022 Winter Simulation Conference (WSC)*, pages 3229–3240, 2022.

[58] Robert Liptser, Albert Nikolaevich Shiriaev, and Albert N Shiryaev. *Statistics of random processes: I. General theory*, volume 1. Springer Science & Business Media, 2001.

[59] Pierre-Yves Madec. Ergodicbsdes and related pdes with neumann boundary conditions under weak dissipative assumptions. *Stochastic Processes and their Applications*, 125(5):1821–1860, 2015.

[60] Xuerong Mao. *Stochastic differential equations and applications*. Elsevier, 2007.

[61] David Martin, Charless Fowlkes, Doron Tal, and Jitendra Malik. A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics. In *Proceedings eighth IEEE international conference on computer vision. ICCV 2001*, volume 2, pages 416–423. IEEE, 2001.

[62] Jean-Jacques Moreau. Proximité et dualité dans un espace hilbertien. *Bulletin de la Société mathématique de France*, 93:273–299, 1965.

[63] Jean Jacques Moreau. Inf-convolution, sous-additivité, convexité des fonctions numériques. *Journal de Mathématiques Pures et Appliquées*, 1970.

[64] Wenlong Mou, Nicolas Flammarion, Martin J Wainwright, and Peter L Bartlett. Improved bounds for discretization of Langevin diffusions: Near-optimal rates without convexity. *Bernoulli*, 28(3):1577–1601, 2022.

[65] Nina Ovcharova. Second-order analysis of the moreau–yosida and the lasry–lions regularizations. *Optimization Methods & Software*, 25(1):109–116, 2010.

[66] Marcelo Pereyra. Proximal markov chain monte carlo algorithms. *Statistics and Computing*, 26: 745–760, 2016.

[67] Jean-Christophe Pesquet, Audrey Repetti, Matthieu Terris, and Yves Wiaux. Learning maximally monotone operators for image recovery. *SIAM Journal on Imaging Sciences*, 14(3): 1206–1237, 2021.

[68] Natesh S Pillai. Optimal scaling for the proximal Langevin algorithm in high dimensions. *Journal of Machine Learning Research*, 25(404):1–32, 2024.

[69] Roland Potthast. A survey on sampling and probe methods for inverse problems. *Inverse Problems*, 22(2):R1, 2006.
[70] Miklós Rásonyi and Kinga Tikosi. On the stability of the stochastic gradient Langevin algorithm with dependent data stream. *Statistics & Probability Letters*, 182:109321, 2022.

[71] Marien Renaud, Jiaming Liu, Valentin De Bortoli, Andrés Almansa, and Ulugbek S Kamilov. Plug-and-play posterior sampling under mismatched measurement and prior models. *The Twelfth International Conference on Learning Representations*, 2023.

[72] Marien Renaud, Julien Hermant, and Nicolas Papadakis. Convergence analysis of a proximal stochastic denoising regularization algorithm, 2024. URL https://arxiv.org/abs/2412.08262.

[73] Marien Renaud, Jean Prost, Arthur Leclaire, and Nicolas Papadakis. Plug-and-play image restoration with stochastic denoising regularization, 2024. URL https://arxiv.org/abs/2402.01779.

[74] Gareth O Roberts and Osnat Stramer. Langevin diffusions and Metropolis-Hastings algorithms. *Methodology and computing in applied probability*, 4:337–357, 2002.

[75] Gareth O Roberts and Richard L Tweedie. Exponential convergence of Langevin distributions and their discrete approximations. *Bernoulli*, pages 341–363, 1996.

[76] R Tyrrell Rockafellar and Roger J-B Wets. *Variational analysis*, volume 317. Springer Science & Business Media, 2009.

[77] Yaniv Romano, Michael Elad, and Peyman Milanfar. The little engine that could: Regularization by denoising (red). *SIAM Journal on Imaging Sciences*, 10(4):1804–1844, 2017.

[78] Litu Rout, Negin Raoof, Giannis Daras, Constantine Caramanis, Alex Dimakis, and Sanjay Shakkottai. Solving linear inverse problems provably via posterior sampling with latent diffusion models. *Advances in Neural Information Processing Systems*, 36:49960–49990, 2023.

[79] Leonid I Rudin, Stanley Osher, and Emad Fatemi. Nonlinear total variation based noise removal algorithms. *Physica D: nonlinear phenomena*, 60(1-4):259–268, 1992.

[80] Adil Salim and Peter Richtarik. Primal dual interpretation of the proximal stochastic gradient Langevin algorithm. *Advances in Neural Information Processing Systems*, 33:786–796, 2020.

[81] Adil Salim, Dmitry Kovalev, and Peter Richtárik. Stochastic proximal Langevin algorithm: Potential splitting and nonasymptotic rates. *Advances in Neural Information Processing Systems*, 32, 2019.

[82] Zakhar Shumaylov, Jeremy Budd, Subhadip Mukherjee, and Carola-Bibiane Schönlieb. Weakly convex regularisers for inverse problems: Convergence of critical points and primal-dual optimisation. *Forty-first International Conference on Machine Learning*, 2024.

[83] Yang Song and Stefano Ermon. Generative modeling by estimating gradients of the data distribution. *Advances in neural information processing systems*, 32, 2019.

[84] D Andersen Sorensen, S Andersen, D Gianola, and I Korsgaard. Bayesian inference in threshold models using gibbs sampling. *Genetics Selection Evolution*, 27(3):229–249, 1995.

[85] Yu Sun, Zihui Wu, Yifan Chen, Berthy T Feng, and Katherine L Bouman. Provable probabilistic imaging using score-based generative priors. *IEEE Transactions on Computational Imaging*, 2024.

[86] J Tachella, D Chen, S Hurault, M Terris, and A Wang. Deepinverse: A deep learning framework for inverse problems in imaging. url: https://deepinv.github.io/deepinv, 2023.

[87] Aki Vehtari, Daniel Simpson, Andrew Gelman, Yuling Yao, and Jonah Gabry. Pareto smoothed importance sampling. *Journal of Machine Learning Research*, 25(72):1–58, 2024.

[88] Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg. Plug-and-play priors for model based reconstruction. In *2013 IEEE global conference on signal and information processing*, pages 945–948. IEEE, 2013.
[89] Cédric Villani et al. *Optimal transport: old and new*, volume 338. Springer, 2008.

[90] Maxime Vono, Nicolas Dobigeon, and Pierre Chainais. Split-and-augmented gibbs sampler—application to large-scale inference problems. *IEEE Transactions on Signal Processing*, 67(6):1648–1661, 2019.

[91] Zhou Wang, Alan C Bovik, Hamid R Sheikh, and Eero P Simoncelli. Image quality assessment: from error visibility to structural similarity. *IEEE transactions on image processing*, 13(4): 600–612, 2004.

[92] Deliang Wei, Peng Chen, and Fang Li. Learning pseudo-contractive denoisers for inverse problems. In *Forty-first International Conference on Machine Learning*, 2024.

[93] Andre Wibisono. Sampling as optimization in the space of measures: The Langevin dynamics as a composite optimization problem. In *Conference on Learning Theory*, pages 2093–3027. PMLR, 2018.

[94] Kaylee Yingxi Yang and Andre Wibisono. Convergence of the inexact Langevin algorithm and score-based generative models in KL divergence. *NeurIPS 2022 Workshop on Score-Based Methods*, 2022.

[95] Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte. Plug-and-play image restoration with deep denoiser prior. *IEEE Transactions on Pattern Analysis and Machine Intelligence*, 44(10):6360–6376, 2021.

[96] Richard Zhang, Phillip Isola, Alexei A Efros, Eli Shechtman, and Oliver Wang. The unreasonable effectiveness of deep features as a perceptual metric. In *Proceedings of the IEEE conference on computer vision and pattern recognition*, pages 586–595, 2018.

[97] Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool. Denoising diffusion models for plug-and-play image restoration. In *Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition*, pages 1219–1229, 2023.

# Appendix

In the Appendix, we provide all the proof for the theoretical results of the paper, additional details on experiments and more intuition on the theoretical objects that appears in the paper. In Appendix A, we detail related works on Langevin MCMC with a focus on Proximal Langevin methods. In Appendix B, we discuss each assumptions in details. In Appendix C, we provide additional experiments and implementation details. In Appendix D, we give a self-contained analysis of the Moreau envelope in the weakly convex setting. In Appendix E, we discuss our notion of strong convexity on non-convex sets. In Appendix F, we prove all the results of Section 2 and give intuition on Markov Chain and Markov Process analysis. In Appendix G, we prove the convergence results of Section 3 for PSGLA.

## A Related works

In this part, we detail related works on Langevin Diffusion with a focus on methods involving proximal steps.

**On the stability of inexact Unadjusted Langevin Algorithm** The authors of [94] study the stability of inexact Unadjusted Langevin Algorithm (iULA defined in equation (2)) with respect to the KL divergence. This study can be seen as complementary to our analysis developed in p-Wasserstein and Total Variation distance. Note that our strategy does not generalize easily to KL divergence as it relies massively on the triangle inequality that is not verified by the KL divergence. The proof strategy of [94] relies on the log-Sobolev inequality for the potential V and bounds between the inexact drift b and $\nabla V$ making this stability analysis less general than our analysis. The stability of iULA has been analysed in the same setting than ours by the authors of [71]. However, their stability analysis is only provided in the 1-Wasserstein and Total Variation distances, and yields a coarser bound than Theorem 1, which in particular implies a discretization error.
**PnP-ULA** The authors of [56] propose a Plug-and-Play strategy based on inexact ULA and the Tweedie formula [31], named PnP-ULA, and defined by

$$X_{k+1} = X_k - \gamma \nabla f(X_k) - \frac{\gamma \lambda}{\epsilon} (X_k - D_{\epsilon}(X_k)) - \frac{\gamma}{\alpha} (X_k - \Pi_S(X_k)) - \sqrt{2\gamma Z_{k+1}}, \quad (16)$$

with $Z_{k+1} \sim N(0, I_d)$, $D_{\epsilon}$ a pretrained denoiser of noise-level $\epsilon > 0$, $\Pi_S$ the projection on the convex set $S \subset \mathbb{R}^d$ and $\lambda > 0$ the regularization weight. The term $\frac{\lambda}{\epsilon}(X_k - D_{\epsilon}(X_k))$ is used to encode the prior information. The projection term $\frac{1}{\alpha}(X_k - \Pi_S(X_k))$, controlled by the weight parameter $\alpha > 0$, is added to ensure that the sampled potential is strongly convex at infinity for convergence properties. This Markov Chain has a slow mixing time as shown on Figure 5 where it needs more than $5.10^5$ iterations to converge.

**MYULA** The authors of [66, 27] propose the Moreau-Yoshida Unadjusted Langevin Algorithm (MYULA) defined by

$$X_{k+1} = X_k - \gamma \nabla f(X_k) - \frac{\gamma}{\lambda} (X_k - \text{Prox}_{\lambda g}(X_k)) + \sqrt{2\gamma Z_{k+1}}$$

to sample $e^{-f-g^{\lambda}}$, with $g^{\lambda}$ the Moreau transform of $g$. This algorithm is based on the Moreau envelope approximation (see more details in Appendix D) to approximate $\nabla g$ using the proximal operator $\text{Prox}_{\lambda g}$. It is particularly relevant if the proximal operator on $g$ can be computed easily. Contrary to PSGLA, in MYULA, the noise is added outside the proximal operator. Recently the authors of [52] consider algorithms based on MYULA to accelerate Langevin samplers. Note that if a pre-trained denoiser is plugged in place of the proximal operator, then we recover PnP-ULA.

**Gibbs sampling** The authors of [90, 17] propose another strategy to sample a composite function, named Gibbs sampling. A Gibbs sampler [37, 84] splits the sampling problem into two easier sampling sub-problems. The first sub-problem encodes the prior information and the second the degradation model, following the Plug-and-Play idea. However no convergence guarantee is given for this algorithm. The authors of [85] propose a similar algorithm with some theoretical results in the non-convex setting. However, they employ annealing techniques, and it is unclear how these affect convergence. Finally, the authors of [57] propose another Gibbs sampler and derive convergence guarantees only in the case of convex potentials.

**Proximal Langevin Monte Carlo** The authors of [5] study the Proximal Langevin Monte Carlo (PLMC) sequence defined by

$$X_{k+1} = \text{Prox}_{\gamma g}(X_k) + \sqrt{2\gamma Z_{k+1}}$$

to sample $e^{-g}$ with $g$ convex everywhere. The authors of [4] study the same sequence and focus on developing convergence results beyond the Lipschitz Gradient continuity. PLMC is a particular case of MYULA in the case of $\lambda = \gamma$. A Metropolis-Hasting mechanism has also been proposed in [68, 18] to accelerate PLMC or MYULA in the case of convex potentials.

**Projected Langevin** A particular case of PSGLA, where the function $g$ is the characteristic function of a convex set $K$, leads to the projected Langevin algorithm defined by

$$X_{k+1} = P_K(X_k - \gamma \nabla f(X_k) + \sqrt{2\gamma Z_{k+1}}),$$

which was introduced over a decade ago [9]. This paper studies the convergence in the log-concave setting, i.e. if the potential $f$ is convex. The author of [53] studies a generalization of this algorithm when $\nabla f$ is estimated randomly. In this paper, convergence guarantees are provided if $K$ is convex and $\nabla f$ is Lipschitz and perturbed by Gaussian noise.

**Symmetrized Langevin Algorithm** The Symmetrized Langevin Algorithm (SLA) defined by $X_{k+1} = \text{Prox}_{\gamma V}(X_k - \gamma \nabla V(X_k) + \sqrt{4\gamma Z_{k+1}})$ has been proposed by the authors of [93]. This algorithm does not sample from a composite potential, but instead targets the distribution proportional to $e^{-V}$. This algorithm is studied for the 2-Wasserstein metric in the case of a strongly convex potential $V$.
**Other articles on PSGLA** The PSGLA (3) was first introduced in [28], with an analysis for strongly convex potentials $V = f + g$. The authors of [81] study a generalization of PSGLA with multiple proximal operators applied in sequence. Their analysis also assume that $f$ is strongly convex and $g$ convex. Another convergence proof with $f$ strongly convex and $g$ convex was derived by making an analogy between PSGLA and the forward-backward scheme [80]. If the proximal mapping can only been approximated, the authors of [32] show that the algorithm is stable in the case of convex potentials $V$. Recently, a particle interaction algorithm based on PSGLA has been proposed [33] to improve the sampling quality. The convergence of this particle interaction algorithm have been derived in the strongly convex setting.

## B Detailed discussion on assumptions

In this section, we discuss in details the different assumptions made in the paper.

* **Assumption 1(i).** Assumption 2(i) is standard to obtain the existence of the strong solution of the stochastic differentiable equation (1), see [50, Theorem 2.5] for instance. It implies in particular that the drift $b$ is sub-linear, i.e. there exists $\|b(x)\| \le \|B(0)\| + L\|x\|$. Moreover the drift $b$ is $L$-weakly convex. If the drift is $C^2$, Assumption 1(i) is equivalent to $-LI_d \le \nabla^2 b \le LI_d$.

* **Assumption 1(ii).** It is called weak dissipativity [23, 59] or distant dissipativity [64] in the stochastic differentiable equation community. In order to have a geometric ergodicity of the Markov Chain, it is necessary to control the drift in the tails of the sampled distribution. Assumption 1(ii) is a relaxation of the strong convexity of $b$ that allows to control the behavior of the MCMC at infinity. In particular, having a form of strong convexity at infinity of the drift, with Assumption 1(ii), ensures that the MCMC is bounded and has all its moments bounded [22, 56]. If the drift is $m$-strongly convex at infinity and $L$-weakly convex then it is weakly dissipative.

* **Assumption 2(i).** It is verified for convex functions $f$ such as data-fidelity terms associated to linear inverse problems with additive Gaussian noise $f(x) = \|Ax - y\|^2$ or non-convex $f$, as it is the case in image despeckling [24] or black-hole interferometric imaging [85]. This condition is standard in posterior sampling [56, 71, 52] or in imaging inverse problems [73, 92].

* **Assumption 2(ii).** Many works have focused on convex regularizations [79, 67, 36]. However, it is known that the non-convexity allows to make better restoration [45]. Therefore, we only assume in this paper that the regularization is $\rho$-weakly convex. It ensures that the proximal operator $\text{Prox}_{\gamma g}$ is defined for $\gamma\rho < 1$. Assumption 2(ii) is then verified for all convex regularization and some state-of-the art deep neural networks regularizers [67, 45, 82, 39].

* **Assumption 3(i).** Assumption 3(i) can be verified by non-smooth functions $g$ such as the characteristic function of a convex set. It is not usual to assume that the function which is a proximal operator is smooth on a subset of $\mathbb{R}^d$ [80, 28, 32]. However, the Moreau envelope study leads naturally to look at the behavior of $g$ on the set $\text{Prox}(\mathbb{R}^d)$ as detailed in Appendix D.7. Moreover, the deep neural network denoiser proposed in [46] has been demonstrated to verify Assumption 3(i). Assumption 3(i) ensures that the drift $b^\gamma$ is Lipschitz to verify Assumption 1(i) (based on Lemma 15 in Appendix D.7). Note that $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ can be a non-convex set if $g$ is not finite everywhere. However, Proposition 2 in Appendix D.7 shows that $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ only differs from being convex by a negligible set.

* **Assumption 3(ii).** Assumption 3(ii) implies that the potential $V$ is strongly convex at infinity, which ensures that the drift $b^\gamma$ verifies Assumption 1(ii). Lemma 16 in Appendix E shows that if $g$ is strongly convex at infinity and globally smooth, then Assumption 3(ii) is verified. Assumption 3(ii) is technical and hard to verify in practice but it is important to study PSGLA. An alternative to this assumption is to modify PSGLA by adding a projection term as in [56] leading to the Projected PSGLA

$$X_{k+1} = \text{Prox}_{\gamma g} \left( X_k - \gamma \nabla f - \frac{\gamma}{\alpha} (X_k - \Pi_K(X_k)) + \sqrt{2\gamma Z_{k+1}} \right),$$

with $\Pi_K$ the projection on the convex set $K \subset \mathbb{R}^d$. As detailed in [56, Appendix F.2], by choosing a small enough $\alpha > 0$, this additional projection term ensures that
the drift of the shadow sequence of Projected PSGLA $\nabla f(y - \gamma\nabla g^{\gamma}(y)) + \nabla g^{\gamma} + \frac{1}{\alpha}(y - \gamma\nabla g^{\gamma}(y) - \Pi_K(y - \gamma\nabla g^{\gamma}(y)))$ verifies Assumption 1(ii). Therefore, it is possible to demonstrate the convergence of Projected PSGLA with a small enough parameter $\alpha > 0$ without Assumption 3(ii). Note that if $X_k$ is an image, typically in $[0, 1]^d$, then choosing $K = [-1, 2]^d$ as in [56] involves the projection term being rarely activated. However, in this paper, we chose not to modify the PSGLA, so Assumption 3(ii) is necessary.

# C More details on experiments

## C.1 Experiments with Gaussian Mixture prior in 2D

We suppose that the prior is a Gaussian Mixture, that can be expressed as

$$p(x) = \sum_{i=1}^{p} p_i \mathcal{N}(x; m_i, \Sigma_i),$$

with $\mathcal{N}(x; m_i, \Sigma_i) = \frac{1}{(2\pi|\Sigma_i|)^{\frac{d}{2}}} e^{-\frac{1}{2}(x-m_i)^T\Sigma_i^{-1}(x-m_i)}$ and $\Sigma_i \in \mathbb{R}^{d \times d}$ a symmetric positive-definite covariance matrixcovariance matrix, $m_i \in \mathbb{R}^d$ the mean of the distribution and $p_i \in [0, 1]$ are the mode weights, such that $\sum_{i=1}^p p_i = 1$.

Then the gradient of the log prior is

$$\nabla \log p(x) = - \frac{\sum_{i=1}^{p} p_i \Sigma_i^{-1} (x - m_i) \mathcal{N}(x; m_i, \Sigma_i)}{\sum_{i=1}^{p} p_i \mathcal{N}(x; m_i, \Sigma_i)}.$$

Thanks to the Tweedie formula [31], the Minimum Mean Square Error (MMSE) denoiser $D_{\epsilon}$ can be computed by

$$D_{\epsilon}(x) = x - \epsilon^2 \nabla \log p_{\epsilon}(x),$$

with $\epsilon$ the noise level and $p_{\epsilon} = p \star \mathcal{N}(0, \epsilon^2 I_d)$ the convolution between $p$ and the centered Gaussian kernel $\mathcal{N}(0, \epsilon^2 I_d)$.

Due to the specific form of $p$, we know that

$$\nabla \log p_{\epsilon}(x) = - \frac{\sum_{i=1}^{p} p_i (\Sigma_i + \epsilon^2 I_d)^{-1} (x - m_i) \mathcal{N}(x; m_i, \Sigma_i + \epsilon^2 I_d)}{\sum_{i=1}^{p} p_i \mathcal{N}(x; m_i, \Sigma_i + \epsilon^2 I_d)}.$$

Thus the closed-form expression of the MMSE denoiser is

$$D_{\epsilon}(x) = x + \epsilon^2 \frac{\sum_{i=1}^{p} p_i (\Sigma_i + \epsilon^2 I_d)^{-1} (x - m_i) \mathcal{N}(x; m_i, \Sigma_i + \epsilon^2 I_d)}{\sum_{i=1}^{p} p_i \mathcal{N}(x; m_i, \Sigma_i + \epsilon^2 I_d)}.$$

For Gaussian Mixture prior, the gradient of the regularization $-\log p$ and the denoiser $D_{\epsilon}$ are known in closed-form. Therefore, we can compute PnP-ULA and PnP-PSGLA exactly.

**Parameters setting** We run PnP-ULA and PnP-PSGLA for a 2D denoising problem, i.e. $d = m = 2$, $A = I_2$ and $\sigma = 1$. PnP-ULA defined in equation (16) is run with a step-size $\gamma = 0.1$, a regularization parameter $\lambda = 1.5$ and a denoiser parameter $\epsilon = 0.5$. In this case, the sampled distribution is log-concave at infinity. Therefore, we choose to set $\alpha = +\infty$, i.e. there is no projection. PSGLA is run with a step-size $\gamma = 0.3$ and a regularization parameter $\lambda = 0.67$. These parameters have been chosen to maximize the performance of each method.

**Metrics** We compute two quantitative metric between the discrete uniform measure supported on the iterates and the exact posterior. The first one is the discrete Wasserstein distance computed using the Python Optimal Transport (POT) library [34]. The second one is the $\ell_2$ distance between the empirical density estimated using Gaussian kernels and the exact posterior density. This $\ell_2$ distance is computed on a square including most of the density ($[-8, 8]^2$ in our experiments).
<table><caption>Table 3: Parameters setting for image inpainting for the different implemented methods. $\epsilon$ is the noise level of the denoiser, $\lambda$ the regularization parameter, $\gamma$ the step-size. The parameter settings for PnP-ULA are those suggested in the original paper [56].</caption><thead><tr><th>Parameters</th><th>RED GSDRUNet</th><th>RED DnCNN</th><th>PnP GSDRUNet</th><th>PnP DnCNN</th><th>PnP-ULA DnCNN</th><th>PSGLA TV</th><th>PSGLA DnCNN</th></tr></thead><tbody><tr><td>255 &times; &epsilon;</td><td>7</td><td>2</td><td>5</td><td>2</td><td>2</td><td>10</td><td>2</td></tr><tr><td>&lambda;</td><td>70,000</td><td>150,000</td><td>0.5</td><td>1</td><td>1/2 + 1/2</td><td>10</td><td>5</td></tr><tr><td>&gamma;</td><td>1.10<sup>-5</sup></td><td>1.10<sup>-5</sup></td><td>1.10<sup>-5</sup></td><td>1.10<sup>-5</sup></td><td>1/3 &times; (1/2 + 1/3 + 1/3)</td><td>&epsilon;<sup>2</sup></td><td>&epsilon;<sup>2</sup></td></tr></tbody></table>

## C.2 More details on image restoration experiments

In this section, we present more details on the experiments for image restoration. First we recall the definition of each compared method. Then we give the parameters setting and additional experiments.

**PnP** We compare PSGLA into method that are mixing information from learning and the physics of degradation, named Plug-and-Play (PnP). This algorithm has been proposed for image restoration by [88] and is defined by

$$X_{k+1} = D_{\epsilon} \left( X_k - \frac{\gamma}{\lambda} \nabla f(X_k) \right).$$

The authors of [67] shows that with a maximal monotone DnCNN denoiser PnP converges. The convergence of this algorithm have been extensively studied, see a detailed review in [45, Part 2]. This algorithm have been designed to optimize and compute an approximation of the most probable image knowing the degraded observation. Thus, PnP is not sampling the posterior distribution and there is no information about uncertainty that can be extracted from this algorithm.

**RED** Another optimization algorithm based mixing prior information and physical information is Regularization by Denoising (RED). This algorithm has been proposed for image restoration by [77] and is defined by

$$X_{k+1} = X_k - \gamma \nabla f(X_k) - \gamma \lambda (X_k - D_{\epsilon}(X_k)).$$

The authors of [45] shows that with a Gradient-Step denoiser RED converges.

**DiffPIR** DiffPIR definition is detailed in Algorithm 2. This algorithm is based on a diffusion model method applied to image restoration. To our knowledge, no theoretical analysis of this algorithm has been provided and it is not clear if this algorithm is solving an optimization problem or sampling a target law. This algorithm requires a denoiser trained with various level of noise. We use the Gradient-Step DRUNet in our experiment with the pre-trained weights provided by the authors of [45]. The precomputed parameters $(\bar{\alpha}_t)_{0<t<T}$, $(\sigma_t)_{0<t<T}$ and $(\rho_t)_{0<t<T}$ are set as proposed in [97]. In our experiments for image inpainting, we set $T = 1000$, $t_{start} = 200$, $\lambda = 0.13$, $\zeta = 0.999$. These parameters have been chosen after running a grid-search to maximize the PSNR on the 3 natural images (butterfly, leaves and starfish) from the set3c dataset [45].

**Algorithm 2** DiffPIR [97]

1: **input:** denoiser $D$, $T > 0$, $y \in \mathbb{R}^m$, $0 < t_{start} < T$, $\zeta > 0$, $(\beta_t)_{0<t<T}$, $\lambda > 0$

2: Initialize $\epsilon_{t_{start}} \sim N(0, I_d)$, pre-calculate $(\bar{\alpha}_t)_{0<t<T}$, $(\sigma_t)_{0<t<T}$ and $(\rho_t)_{0<t<T}$

3: $x_{t_{start}} = \sqrt{\bar{\alpha}_{t_{start}}}y + \sqrt{1 - \bar{\alpha}_{t_{start}}}\epsilon_{t_{start}}$

4: **for** $t = t_{start}$, $t_{start} - 1, \dots, 1$ **do**

5: $x_0^t \leftarrow D_{\sigma_t}(x_t)$

6: $\hat{x}_0^t \leftarrow \text{Prox}_{2f(\cdot)/\rho_t}(x_0^t)$

7: $\hat{\epsilon} \leftarrow (x_t - \sqrt{\bar{\alpha}_t}\hat{x}_0^t) / \sqrt{1 - \bar{\alpha}_t}$

8: $\epsilon_t \leftarrow N(0, I_d)$

9: $x_{t_1} \leftarrow \sqrt{\bar{\alpha}_t}\hat{x}_0^t + \sqrt{1 - \bar{\alpha}_t}(\sqrt{1 - \bar{\zeta}}\epsilon + \sqrt{\bar{\zeta}}\epsilon_t)$

10: **end for**

**Parameters setting** In Table 3, we provide the parameter choices used in our experiments for the implemented methods.
<figure><img src="image_4.png" /><figcaption>Figure 3: Qualitative result for image inpainting with 50% masked pixels and a noise level of $\sigma = 1/255$ with PnP, DiffPIR and PnP-PSGLA methods.</figcaption></figure>

**Computational time** In Table 2, we give the mean computation time to restore one image on a GPU NVIDIA A100 Tensor Core. The total computational time to generate this table is thus around 27 hours of computational time on a GPU NVIDIA A100 Tensor Core. Note that the PnP-ULA requires most of the computational resources for a poor restoration quality. The total computational time for PnP-PSGLA with DnCNN to restore the CBSD68 dataset is 136 minutes. The computational cost of the entire project, including the parameter grid search, has not been rigorously computed. We estimate it to be around 40 hours of computational time on an NVIDIA A100 Tensor Core GPU.

**On the denoiser choice** In this paper, we only use pre-trained denoisers with open-source weights. We focus on a few denoisers that are known to be stable and provide good restoration. First the DnCNN with the weights trained by [67] is a maximal monotone operator, thus it is a proximal operator of some convex regularization $g$. Therefore our theory (Theorem 3) applies with this denoiser. We also run the methods with the Gradient Step denoiser [45], a denoiser based on the DRUNet architecture [95], that is known to provide state-of-the-art restoration for Plug-and-Play methods. As this denoiser is trained on natural images with various noise level $\epsilon \in [0, 50]/255$, we can use it for the DiffPIR algorithm. The open access weights of the DnCNN and the GSDRUNet used in our experiment for natural color images can be found in the library DeepInv [86] in this link. In order to reduce the computational time of PnP-PSGLA, we also use the TV-denoiser [79]. The Total-Variation regularization has not a closed form proximal operator so we use an iterative algorithm (implemented
<figure><img src="image_5.png" /><figcaption>Figure 4: Restoration of PSGLA with Prox DRUNet for 50% missing pixels and $\sigma = 1/255$. We chose the aggressive parameters $\epsilon = 50/255$, $\lambda = 5000$ to over-smooth the solution for more stability. PSGLA is run twice on the same observation with two seed in the randomness of the algorithm. Note that the Markov Chain explores outside the image space $[0, 1]^d$ and the instability depends on the realization of the noise. This algorithm is not stable.</figcaption></figure>

in the library DeepInv [86]) with 10 inner iterations, which is enough in our experiments to obtain the best restoration with TV-denoiser. By using such an approximation of PnP-PSGLA with the Total Variation denoiser, we recover the restoration algorithm proposed in [32].

### C.2.1 Why we do not use Prox-DRUNet with PSGLA

Based on the architecture of Gradient-Step DRUNet, the authors of [46] proposed a denoiser that is a proximal operator. This proximal operator satisfies Assumption 2(ii) and Assumption 3(i). The DRUNet architecture [95] is also known to provide a state-of-the-art denoiser. Thus Prox-DRUNet seems to be an ideal denoiser candidate for applying our theory and having empirical performances. However, we observe in our experiments some instabilities with the weight provided by the authors [46].

On Figure 4, we observe the restored image and the standard deviation of the Markov Chain run twice on the same observation with two seeds in the randomness of the PSGLA. A different area of the image is not restored properly and the Markov Chain explores outside the image space $[0, 1]^d$ with a large standard deviation in these areas. The mosaic motif that appears may be due to periodic colorization outside $[0, 1]^d$. The fact that the problematic areas depend of the realization of the noise indicates that the issue lies with the instability of the algorithm, not with the conditioning of the target distribution. Then, adding a projection in $[0, 1]^d$ at each iteration of $X_k$ does not prevent from these instabilities. Moreover, these instabilities appear for all the parameters setting that we test in our large grid-search. Finally, compare with DnCNN, trained with a fix level of noise, PSGLA with Prox-DRUNet has more parameters to fine-tune. For all these arguments, we chose to not use PSGLA with the Prox-DRUNet denoiser.

### C.2.2 PSGLA seems to effectively sample the posterior law

For image inverse problem, as the posterior distribution is unknown, it is not possible to evaluate if an algorithm effectively samples the posterior distribution. A standard strategy [56, 52] to have an indication of the sampling performance of a particular algorithm is to look at the standard-deviation of the Markov Chain.

In Figure 5, we present the mean and the standard deviation of the PSGLA Markov Chain for different numbers of iterations. We observe that when the Markov Chain is too short ($N = 10^4$) the mean is relevant but the standard deviation is not informative. On the other hand, when the number of iterations is large ($N = 10^6$), the standard deviation is both informative and relevant. One can thus observe that the image edges are more uncertain than uniform areas, such as the sky. The posterior distribution in image inverse problems is known to concentrate uncertainty in the edges [56, 52]. This is a good indication that the PSGLA might sample the posterior distribution in practice. In order to compute relevant uncertainty information, $N = 10^5$ seems to be a good compromise between computational efficiency and accuracy.
<figure><img src="image_6.png" /><figcaption>Figure 5: PnP-ULA and PSGLA with DnCNN for 50% missing pixels and $\sigma = 1/255$ with various number of iterations $N \in \{10^4, 10^5, 10^6\}$. The standard deviation of the Markov Chains are shown for each number of iterations and the evolution of the PSNR for $N \in [0, 10^6]$. Note that the PnP-ULA Markov Chain is slow to converge. The effective convergence seems to occur around $5 \cdot 10^5$ iterations. Note that the PSGLA Markov Chain seems to well explore the posterior distribution as the uncertainty relies on the edges and the background, as the blue uniform sky, is certain.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: PSGLA with DnCNN for 50% missing pixels and $\sigma = 1/255$ with two number of iterations $N \in \{10^4, 10^5\}$. We run PSGLA with 100 random seeds in the algorithm randomness for $N = 10^4$ and 10 random seeds in the algorithm randomness for $N = 10^5$ on the same observation. We can compute an estimator of the standard deviation of the MMSE estimator. This MMSE restored image is expected to be deterministic, so we expect this standard deviation to be significatly smaller than the standard deviation of the Markov Chain itself (see in Figure 5). Note that the standard deviation of the estimator is 10 times smaller than the standard deviation of PSGLA Markov Chain.</figcaption></figure>

In Figure 6, we show the standard deviation of the MMSE estimator. We observe that the MMSE estimator has a small standard deviation compared with the standard deviation of the PSGLA Markov Chain. Moreover the standard deviation of the MMSE estimator decreases with the number of iterations. This observation support the fact the MMSE estimator is deterministic and that the PSGLA effectively solves a sampling problem.

### C.2.3 PnP-ULA is slow to converge

In Table 2, we observe that the performance of PnP-ULA is significantly lower than that of the other methods. However, the authors of [56, 71] show that PnP-ULA can provide good restoration. This apparent paradox is explained by the fact that PnP-ULA does not succeed to converge in $N = 10^5$ iterations. On Figure 5, we observe that the Markov Chain indeed converges in approximately $5 \cdot 10^5$ iterations. In order to save computational resources, we choose to not run with such a large number of iterations PnP-ULA on the whole CBSD68 dataset to generate Table 2.
# D On the Moreau envelope for weakly convex function

In this section, we recall and detail results on the Moreau envelope, also called Moreau-Yosida envelop, in the case of weak convexity. For clarity of this paper, we detail all the proofs. More details on the Moreau envelope of convex function can be found in [76, Part 1.G] or in [62]; and on the inf-convolution in [63].

All the following results are known by experts in the field. However, to our knowledge, there is no reference that states and proves all the results given in this appendix, especially for **weakly convex functions**. Different proofs can be found in [62, 76, 44, 3, 40, 47]. For simplicity, we choose not to introduce the Clarke sub-differential [16] which is necessary to analyse the Moreau envelope with non-differentiable functions. Therefore, some of the proofs (Lemma 11 and Lemma 13) are only given for differentiable functions, references being given for the general case.

## D.1 Definitions

We first define the notion of inf-convolution as well as the Moreau envelope.

**Definition 1.** The inf-convolution between two functions $f, g: \mathbb{R}^d \to \mathbb{R}$ is defined for $x \in \mathbb{R}^d$ by

$$ (f \Box g)(x) = \inf_{y \in \mathbb{R}^d} f(y) + g(x - y), \qquad (17) $$

with the convention $(f \Box g)(x) = -\infty$ if the right-hand-side of equation (17) is not bounded from below.

**Definition 2.** For $\gamma > 0$, the Moreau envelope of $g: \mathbb{R}^d \to \mathbb{R}$, noted $g^\gamma$ is defined for $x \in \mathbb{R}^d$ by

$$ g^{\gamma}(x) = \left( g \Box \frac{1}{2\gamma} \| \cdot \|^2 \right) (x) = \inf_{y \in \mathbb{R}^d} \frac{1}{2\gamma} \|x - y\|^2 + g(y). \qquad (18) $$

The Moreau envelope is finite on $\mathbb{R}^d$ if $g$ is $\rho$-weakly convex, i.e. $\rho > 0$ and $g + \frac{\rho}{2} \| \cdot \|^2$ is convex, with $\gamma\rho < 1$. In fact, the function $y \mapsto \frac{1}{2\gamma} \|x - y\|^2 + g(y)$ is $(\frac{1}{\gamma} - \rho)$ strongly convex so it admits a unique minimum for $\gamma\rho < 1$. The minimal point is called the proximal point $\text{Prox}_{\gamma g}(x) = \arg\min_{y \in \mathbb{R}^d} \frac{1}{2\gamma} \|x - y\|^2 + g(y)$. The Moreau envelope is an approximation of $g$ in the sense that $\forall x \in \mathbb{R}^d$, $\lim_{\gamma \to 0} g^\gamma(x) = g(x)$ (see Lemma 2) and it is a regularization as $g^\gamma$ is differentiable (see Lemma 12) even if $g$ is not.

## D.2 On the dependence of $g^\gamma$ in $\gamma > 0$

The following result shows the monotonicity of $g^\gamma$ in $\gamma$ and the fact that the Moreau envelope is an approximation of $g$ when $\gamma \to 0$.

**Lemma 2.** If $g$ is $\rho$-weakly convex and continuous, then for $x \in \mathbb{R}^d$, $\gamma \in (0, \frac{1}{\rho}) \to g^\gamma(x)$ is decreasing and

$$ \lim_{\gamma \to 0} g^\gamma(x) = g(x). \qquad (19) $$

Moreover, $\gamma \in (0, \frac{1}{\rho}) \to g(\text{Prox}_{\gamma g}(x))$ is decreasing and

$$ \lim_{\gamma \to 0} g(\text{Prox}_{\gamma g}(x)) = g(x). \qquad (20) $$

Lemma 2 can be generalized in the case where $g$ is lower semicontinuous instead of continuous.

*Proof.* We need $\gamma \le \frac{1}{\rho}$ to ensure that the Moreau envelope is well defined. We obtain from the definition of the Moreau envelope (see relation (18)) that $\gamma \to g^\gamma(x)$ is decreasing and $g^\gamma(x) \le g(x)$. The rest of the proof is a generalization of [3, Proposition 12.33] in the weakly convex case.

For $x, y \in \mathbb{R}^d$, we denote by $p = \text{Prox}_{\gamma g}(x)$ and, with $\alpha \in (0, 1)$, $p_\alpha = \alpha y + (1 - \alpha)p$. By definition of the proximal operator, we have

$$ g(p) + \frac{1}{2\gamma} \|x - p\|^2 \le g(p_\alpha) + \frac{1}{2\gamma} \|x - p_\alpha\|^2. $$
By the previous inequality and the $\rho$-weak convexity of $g$, we get

$$ 
\begin{aligned}
g(p) &\leq g(p_{\alpha}) + \frac{1}{2\gamma} \|x - p_{\alpha}\|^2 - \frac{1}{2\gamma} \|x - p\|^2 \\
&\leq \alpha g(y) + (1 - \alpha)g(p) + \frac{\rho}{2}\alpha(1 - \alpha)\|y - p\|^2 - \frac{\alpha}{\gamma}\langle x - p, y - p \rangle + \frac{\alpha^2}{2\gamma}\|y - p\|^2.
\end{aligned}
$$

By rearranging the terms, we get

$$ \alpha \langle x - p, y - p \rangle \leq \alpha \gamma g(y) - \alpha \gamma g(p) + \frac{\alpha}{2} (\alpha + \rho \gamma (1 - \alpha)) \|y - p\|^2. $$

By dividing by $\alpha$ and taking $\alpha \to 0$, we get

$$ \langle x - p, y - p \rangle \leq \gamma g(y) - \gamma g(p) + \frac{\gamma \rho}{2} \|y - p\|^2. $$

For $\frac{1}{\rho} > \mu \geq \gamma$, we define $q = \text{Prox}_{\mu g}(x)$. By the same reasoning, we have that $\forall y \in \mathbb{R}^d$:

$$ \langle x - q, y - q \rangle \leq \gamma g(y) - \gamma g(q) + \frac{\mu \rho}{2} \|y - q\|^2. $$

Taking respectively $y = q$ and $y = p$ in the two previous inequalities, we obtain

$$ \langle q - p, x - p \rangle \leq \gamma g(q) - \gamma g(p) + \frac{\gamma \rho}{2} \|q - p\|^2 $$

$$ \langle p - q, x - q \rangle \leq \mu g(p) - \mu g(q) + \frac{\mu \rho}{2} \|p - q\|^2. $$

By summing these two relations and rearranging the terms, we get

$$ \left(1 - \frac{\rho}{2}(\mu - \gamma)\right) \|p - q\|^2 \leq (\mu - \gamma)(g(p) - g(q)). $$

Therefore, if $\gamma, \mu \in (0, \frac{1}{\rho})$, for $\mu \geq \gamma$, we get $g(\text{Prox}_{\mu g}(x)) \leq g(\text{Prox}_{\gamma g}(x))$. It proves that $\gamma \in (0, \frac{1}{\rho}) \to g(\text{Prox}_{\gamma g}(x))$ is decreasing.

Now, we prove that $g^{\gamma}(x) \to g(x)$ when $\gamma \to 0$. We denote $M = \sup_{\gamma \in (0, \frac{1}{\rho})} g^{\gamma}(x)$, we have $\forall \gamma \in (0, \frac{1}{\rho})$,

$$ M \geq g^{\gamma}(x) = g(\text{Prox}_{\gamma g}(x)) + \frac{1}{2\gamma} \|x - \text{Prox}_{\gamma g}(x)\|^2. $$

Since, $y \to g(y) + \rho \|x - y\|^2$ is $\rho$-strongly convex thus coercive and given that $\text{Prox}_{\gamma g}(x)$ has a value lower than $M$, we get that $m = \sup_{\gamma \in (0, \frac{1}{\rho})} \text{Prox}_{\gamma g}(x) < +\infty$. Moreover, the function $g + \frac{\rho}{2} \|\cdot\|^2$ is convex and lower bounded by an affine function. So there exist $a_1 \in \mathbb{R}^d$ and $a_2 \in \mathbb{R}$, such that $\forall x \in \mathbb{R}^d$, $g(x) + \frac{\rho}{2} \|x\|^2 \geq \langle a_1, x \rangle + a_2$. Combining the previous facts, we get

$$ 
\begin{aligned}
M &\geq g(\text{Prox}_{\gamma g}(x)) + \frac{1}{2\gamma} \|x - \text{Prox}_{\gamma g}(x)\|^2 \\
&\geq \langle a_1, \text{Prox}_{\gamma g}(x) \rangle + a_2 - \frac{\rho}{2} \|\text{Prox}_{\gamma g}(x)\|^2 + \frac{1}{2\gamma} \|x - \text{Prox}_{\gamma g}(x)\|^2 \\
&\geq -\|a_1\|m + a_2 - \frac{\rho}{2}m^2 + \frac{1}{2\gamma} \|x - \text{Prox}_{\gamma g}(x)\|^2.
\end{aligned}
$$

So we get $\lim_{\gamma \to 0} \|x - \text{Prox}_{\gamma g}(x)\| = 0$. Finally, by the continuity of $g$, we have

$$ g(x) \geq \lim_{\gamma \to 0} g^{\gamma}(x) \geq \lim_{\gamma \to 0} g(\text{Prox}_{\gamma g}(x)) = g(x), $$

which proves the desired result.

□
### D.3 On the convexity and weak-convexity of the Moreau envelope

In this part, we prove that the Moreau envelope preserves the convexity and weak-convexity properties.

We will first need the general result of Lemma 3.

**Lemma 3.** If the function $(x, y) \in A \times B \to f(x, y)$, with A, B convex sets, is convex on the variable $(x, y) \in A \times B$, then the function $x \in A \mapsto \inf_{y \in B} f(x, y)$ is convex.

*Proof.* By the convexity of $f$, we get that $\forall t \in [0, 1]$, $x, x' \in A$, $y, y' \in B$,

$$f(tx + (1-t)x', ty + (1-t)y') \leq tf(x, y) + (1-t)f(x', y').$$

By taking the infimum on $y, y' \in B$, we get

$$\inf_{y,y' \in B} f(tx + (1-t)x', ty + (1-t)y') \leq t \inf_{y \in B} f(x, y) + (1-t) \inf_{y' \in B} f(x', y').$$

Because B is convex, $\{ty+(1-t)y'|y, y' \in B\} \subset B$, so $\inf_{y,y' \in B} f(tx+(1-t)x', ty+(1-t)y') \geq \inf_{y \in B} f(tx+(1-t)x', y)$ and we get the convexity of $x \in A \mapsto \inf_{y \in B} f(x, y)$. $\square$

**Lemma 4.** If $g$ is convex, then for all $\gamma > 0$, $g^\gamma$ is convex.

*Proof.* We denote $v(x, y) = \frac{1}{2\gamma}\|x - y\|^2 + g(y)$. For $x, x', y, y' \in \mathbb{R}^d$, we have

$$v(\lambda x + (1 - \lambda)x', \lambda y + (1 - \lambda)y') \\ = \frac{1}{2\gamma} \|\lambda(x - y) + (1 - \lambda)(x' - y')\|^2 + g(\lambda y + (1 - \lambda)y').$$

From the convexity of $\|\cdot\|^2$ and $g$, we get that $v$ is convex in $(x, y) \in \mathbb{R}^d \times \mathbb{R}^d$. Then by applying Lemma 3, we get that $g^\gamma$ is convex. $\square$

In order to study the weak convexity of the Moreau envelope, we will need the following technical lemma.

**Lemma 5.** If $g$ is $\rho$-weakly convex, and denoting $g_\rho = g + \frac{\rho}{2}\|\cdot\|^2$ which is convex, we have

$$g^{\gamma}(x) = g_{\rho}^{\frac{\gamma}{1-\gamma\rho}} \left( \frac{x}{1-\rho\gamma} \right) - \frac{\rho}{2(1-\rho\gamma)} \|x\|^2,$$

and

$$\mathrm{Prox}_{\gamma g}(x) = \mathrm{Prox}_{\frac{\gamma}{1-\gamma\rho} g_{\rho}} \left( \frac{x}{1-\gamma\rho} \right).$$

*Proof.* Let us denote as $g_\rho$ the convex function $g_\rho = g + \frac{\rho}{2}\|\cdot\|^2$. One has $g^\gamma = (g_\rho - \frac{\rho}{2}\|\cdot\|^2)\square\frac{1}{2\gamma}\|\cdot\|^2$, so

$$g^{\gamma}(x) = \inf_{y} g_{\rho}(y) - \frac{\rho}{2} \|y\|^2 + \frac{1}{2\gamma} \|x - y\|^2 \quad (21)$$

$$= \inf_{y} \left( g_{\rho}(y) + \frac{1 - \rho\gamma}{2\gamma} \left( \|y\|^2 - \frac{2}{1 - \rho\gamma} \langle x, y \rangle + \frac{1}{(1 - \rho\gamma)^2} \|x\|^2 \right) \right) \\ + \frac{1}{2\gamma} \left( 1 - \frac{1}{(1 - \rho\gamma)} \right) \|x\|^2$$

$$= \inf_{y} \left( g_{\rho}(y) + \frac{1 - \rho\gamma}{2\gamma} \left\| y - \frac{1}{1 - \rho\gamma} x \right\|^2 \right) - \frac{\rho}{2(1 - \rho\gamma)} \|x\|^2 \quad (22)$$

$$= g_{\rho}^{\frac{\gamma}{1-\gamma\rho}} \left( \frac{x}{1-\rho\gamma} \right) - \frac{\rho}{2(1-\rho\gamma)} \|x\|^2.$$

Then, by the definition of the proximal operator as the minimum point $y \in \mathbb{R}^d$ in the previous optimizations problems (21) and (22), we get

$$\mathrm{Prox}_{\gamma g}(x) = \mathrm{Prox}_{\frac{\gamma}{1-\gamma\rho} g_{\rho}} \left( \frac{x}{1-\gamma\rho} \right). \quad \square$$
**Lemma 6.** If $g$ is $\rho$-weakly convex, with $\gamma\rho < 1$, then $g^\gamma$ is $\frac{\rho}{1-\gamma\rho}$-weakly convex.

Lemma 6 generalizes Lemma 4, for $\rho > 0$. Note that the Moreau envelope of $g$ is less weakly convex that $g$ because $\frac{\rho}{1-\gamma\rho} > \rho$.

*Proof.* We give two proofs of Lemma 6. The first one is more technical and requires $g$ to be twice differentiable, based on an Hessian computation. This proof will be useful for future analysis. The second proof is more elegant and straightforward.

**Proof in the case where $g$ is twice differentiable** We introduce $v(x, y) = \frac{1}{2\gamma}\|x - y\|^2 + g(y)$. By definition, $g^\gamma(x) = \inf_{y \in \mathbb{R}^d} v(x, y)$.

For $\lambda > 0$, we study the Hessian of the function $x, y \in \mathbb{R}^d \to v(x, y) + \frac{\lambda}{2}\|x\|^2$ that writes

$$M(x, y) = \begin{pmatrix} \left(\frac{1}{\gamma} + \lambda\right) I_d & -\frac{1}{\gamma} I_d \\ -\frac{1}{\gamma} I_d & \frac{1}{\gamma} I_d + \nabla^2 g(y) \end{pmatrix}.$$

For $z, t \in \mathbb{R}^d$, by the $\rho$-weakly convexity of $g$, we have

$$(z^T \quad t^T) M(x, y) \begin{pmatrix} z \\ t \end{pmatrix} = \left(\frac{1}{\gamma} + \lambda\right) \|z\|^2 - \frac{2}{\gamma} z^T t + \frac{1}{\gamma} \|t\|^2 + t^T \nabla^2 g(y)t \\ \geq \left(\frac{1}{\gamma} + \lambda\right) \|z\|^2 - \frac{2}{\gamma} z^T t + \left(\frac{1}{\gamma} - \rho\right) \|t\|^2.$$

Moreover, $\forall \alpha > 0$, $2z^T t \le 2\|z\|\|t\| \le \frac{1}{\alpha}\|z\|^2 + \alpha\|t\|^2$. By applying this inequality with $\alpha = 1 - \rho\gamma > 0$, we get

$$(z^T \quad t^T) M(x, y) \begin{pmatrix} z \\ t \end{pmatrix} \\ \geq \left(\frac{1}{\gamma} + \lambda\right) \|z\|^2 - \frac{1}{(1 - \rho\gamma)\gamma} \|z\|^2 - \frac{1}{\gamma} \|t\|^2 + \left(\frac{1}{\gamma} - \rho\right) \|t\|^2 \\ \geq \left(\lambda - \frac{\rho}{1 - \rho\gamma}\right) \|z\|^2.$$

So, we get that for $\lambda = \frac{\rho}{1-\rho\gamma}$, the matrix $M(x, y)$ is positive. Then the function $(x, y) \in \mathbb{R}^d \times \mathbb{R}^d \mapsto v(x, y) + \frac{\rho}{2(1-\rho\gamma)}\|x\|^2$ is convex in $(x, y) \in \mathbb{R}^d \times \mathbb{R}^d$. From Lemma 3, we get that $g^\gamma$ is $\frac{\rho}{1-\rho\gamma}$-weakly convex.

**General proof** Let us denote as $g_\rho$ the convex function $g_\rho = g + \frac{\rho}{2}\|\cdot\|^2$. By Lemma 5, we get

$$g^{\gamma}(x) = g_{\rho}^{\frac{\gamma}{1-\gamma\rho}} \left( \frac{x}{1-\rho\gamma} \right) - \frac{\rho}{2(1-\rho\gamma)} \|x\|^2$$

Since $g_\rho$ is convex, so does $g_\rho^{\frac{\gamma}{1-\gamma\rho}}(\frac{x}{1-\rho\gamma})$ for $\gamma\rho < 1$. Hence we get that $g^\gamma + \frac{\rho}{2(1-\rho\gamma)}\|\cdot\|^2$ is convex, so that $g^\gamma$ is $\frac{\rho}{1-\rho\gamma}$-weakly convex for $\gamma\rho < 1$. $\square$

### D.4 The intrinsic link between the Moreau envelope and the convex conjugate

In this part, we introduce the convex conjugate and shows its link with the Moreau envelope. In particular, Lemma 9 shows a duality result between the convex conjugate and the Proximal operator.

**Definition 3.** The convex conjugate of a proper function $f: \mathbb{R}^d \to \mathbb{R}$, denoted by $f^*$, is defined for $x \in \mathbb{R}^d$ by

$$f^*(x) = \sup_{y \in \mathbb{R}^d} \langle x, y \rangle - f(y).$$

The notion of convex conjugate of $f$ is related with the inf-convolution and the Moreau envelope, as shown by the following Lemmas 7 to 9.
**Lemma 7.** For $f: \mathbb{R}^d \to \mathbb{R}$, a proper function, we have the following properties

(i) $f^*$ is convex.

(ii) For $f, g: \mathbb{R}^d \to \mathbb{R}$ convex, then $(f \Box g)^* = f^* + g^*$.

(iii) $f: \mathbb{R}^d \to \mathbb{R}^d$ is convex if and only if $f^{**} = f$.

(iv) $(\frac{\alpha}{2} \| \cdot \|^2)^* = \frac{1}{2\alpha} \| \cdot \|^2$.

Note that, by Lemma 7(i), $f^*$ is convex even if $f$ is non-convex. Moreover, Lemma 7(ii) shows the compatibility of the convex conjugate with the inf-convolution operation.

*Proof. (i)* For $x, y \in \mathbb{R}^d$ and $\lambda \in [0, 1]$, we have

$$ 
\begin{aligned}
f^*(\lambda x + (1-\lambda)y) &= \sup_{z \in \mathbb{R}^d} \langle \lambda x + (1-\lambda)y, z \rangle - f(z) \\
&= \sup_{z \in \mathbb{R}^d} \lambda (\langle x, z \rangle - f(z)) + (1-\lambda) (\langle y, z \rangle - f(z)) \\
&\leq \lambda \sup_{z \in \mathbb{R}^d} (\langle x, z \rangle - f(z)) + (1-\lambda) \sup_{z \in \mathbb{R}^d} (\langle y, z \rangle - f(z)) \\
&\leq \lambda f^*(x) + (1-\lambda) f^*(y),
\end{aligned}
$$

which shows the convexity of $f^*$.

(ii) For $x \in \mathbb{R}^d$,

$$ 
\begin{aligned}
(f \Box g)^*(x) &= \sup_{y \in \mathbb{R}^d} \langle x, y \rangle - (f \Box g)(y) \\
&= \sup_{y \in \mathbb{R}^d} \langle x, y \rangle - \inf_{z \in \mathbb{R}^d} f(z) + g(y-z) \\
&= \sup_{y \in \mathbb{R}^d} \sup_{z \in \mathbb{R}^d} \langle x, y \rangle - f(z) - g(y-z) \\
&= \sup_{z \in \mathbb{R}^d} -f(z) + \sup_{y \in \mathbb{R}^d} \langle x, y \rangle - g(y-z) \\
&= \sup_{z \in \mathbb{R}^d} -f(z) + \langle x, z \rangle + \sup_{y \in \mathbb{R}^d} \langle x, y-z \rangle - g(y-z) \\
&= \sup_{z \in \mathbb{R}^d} -f(z) + \langle x, z \rangle + g^*(x) \\
&= f^*(x) + g^*(x).
\end{aligned}
$$

(iii) If $f^{**} = f$, then by point (i), $f$ is convex.

If $f$ is convex, we have

$$ f(x) = \sup_{(a,b) \in \Sigma} \langle a, x \rangle + b, $$

with $\Sigma = \{(a, b) \in \mathbb{R}^d \times \mathbb{R} \mid \forall x \in \mathbb{R}^d, \langle a, x \rangle + b \leq f(x)\}$. $(a, b) \in \Sigma$, if and only if $\forall x \in \mathbb{R}^d$, $\langle a, x \rangle - f(x) \leq -b$, i.e. $-b \geq f^*(a)$. Then, we have

$$ 
\begin{aligned}
f(x) &= \sup_{(a,b) \in \Sigma} \langle a, x \rangle + b = \sup_{a \in \mathbb{R}^d, b \leq -f^*(a)} \langle a, x \rangle + b \\
&= \sup_{a \in \mathbb{R}^d} \langle a, x \rangle - f^*(a) = f^{**}(x).
\end{aligned}
$$

(iv) For $x \in \mathbb{R}^d$, $(\frac{\alpha}{2} \|x\|^2)^* = \sup_{y \in \mathbb{R}^d} \langle x, y \rangle - \frac{\alpha}{2} \|y\|^2$. The sup is reached at $y = \frac{x}{\alpha}$, so that $(\frac{\alpha}{2} \| \cdot \|^2)^* = \frac{1}{2\alpha} \| \cdot \|^2$.

<figure><img src="image_8.png" /></figure>

**Lemma 8.** For $f: \mathbb{R}^d \to \mathbb{R}$ a proper function, we have the following properties

(i) $\forall x, y \in \mathbb{R}^d$, $f(x) + f^*(y) \geq \langle x, y \rangle$.
(ii) For $f$ convex and differentiable and $x \in \mathbb{R}^d$, with $y = \nabla f(x)$, we have $f(x) + f^*(y) = \langle x, y \rangle$.

(iii) For $f$ convex and differentiable, we have $y = \nabla f(x)$ if and only if $x = \nabla f^*(y)$.

(iv) For $f$ convex and differentiable, $f$ is $\mu$-strongly convex if and only if $f^*$ is $1/\mu$-smooth.

Lemma 8(iii) can be reformulated as $\nabla f \circ \nabla f^* = I_d$.

*Proof. (i)* For $x, y \in \mathbb{R}^d$

$$f(x) + f^*(y) = \sup_{z \in \mathbb{R}^d} \langle y, z \rangle + f(x) - f(z) \geq \langle y, x \rangle.$$

(ii) If $y = \nabla f(x)$, then by the convexity of $f$, we get, $\forall z \in \mathbb{R}^d$,

$$\begin{aligned}
f(z) &\geq f(x) + \langle y, z - x \rangle \\
\langle x, y \rangle - f(x) &\geq \langle y, z \rangle - f(z) \\
\langle x, y \rangle - f(x) &\geq f^*(y).
\end{aligned}$$

So, by definition of $f^*$, we get $f^*(y) = \langle x, y \rangle - f(x)$.

(iii) If $y = \nabla f(x)$, by the point (ii), we have $f(x) + f^*(y) = \langle x, y \rangle$. Then for all $t \in \mathbb{R}^d$

$$\begin{aligned}
f^*(t) &= \sup_{s \in \mathbb{R}^d} \langle t, s \rangle - f(s) \\
&\geq \langle t, x \rangle - f(x) \\
&\geq \langle t - y, x \rangle + \langle y, x \rangle - f(x) \\
&\geq \langle t - y, x \rangle + f^*(y).
\end{aligned}$$

Thanks to Lemma 7(i) $f^*$ is convex. Hence we get that $x = \nabla f^*(y)$. The reverse is true because $f$ is convex so $f^{**} = f$ by Lemma 7(iii).

(iv) We denote $u(x, y) = \langle x, y \rangle - f(y)$. The $\mu$-strong convexity of $f$ implies that $\forall x, y \in \mathbb{R}^d$,

$$\langle \nabla f(x) - \nabla f(y), x - y \rangle \geq \mu \|x - y\|^2. \quad (23)$$

Thus $\nabla f: \mathbb{R}^d \to \mathbb{R}^d$ is injective. As $u$ is strongly concave w.r.t. $y$, it admits a unique maximal point denoted by $y_0$. On this point the optimal condition gives $\nabla_y u(x, y_0) = 0$, so $\nabla f(y_0) = x$. Therefore, $\nabla f$ is surjective, so it is a bijective function and

$$f^*(x) = \langle x, (\nabla f)^{-1}(x) \rangle - f((\nabla f)^{-1}(x)). \quad (24)$$

By applying Equation (23) with $(\nabla f)^{-1}$, we get $\forall x, y \in \mathbb{R}^d$,

$$\langle x - y, (\nabla f)^{-1}(x) - (\nabla f)^{-1}(y) \rangle \geq \mu \|(\nabla f)^{-1}(x) - (\nabla f)^{-1}(y)\|^2.$$

The Cauchy-Schwarz inequality then gives

$$\| (\nabla f)^{-1}(x) - (\nabla f)^{-1}(y) \| \leq \frac{1}{\mu} \|x - y\|. \quad (25)$$

Using (24), we have for $x, h \in \mathbb{R}^d$

$$\begin{aligned}
f^*(x+h) - f^*(x) &= \langle x+h, (\nabla f)^{-1}(x+h) \rangle - f((\nabla f)^{-1}(x+h)) - \langle x, (\nabla f)^{-1}(x) \rangle \\
&\quad + f((\nabla f)^{-1}(x)) \\
&= \langle x, (\nabla f)^{-1}(x+h) - (\nabla f)^{-1}(x) \rangle + \langle h, (\nabla f)^{-1}(x+h) \rangle \\
&\quad - f((\nabla f)^{-1}(x+h)) + f((\nabla f)^{-1}(x)).
\end{aligned}$$

By denoting $u = (\nabla f)^{-1}(x+h) - (\nabla f)^{-1}(x)$, we know that $u = O(\|h\|)$ thanks to Equation (25). Thus, we get

$$\begin{aligned}
f^*(x+h) - f^*(x) &= \langle x, u \rangle + \langle h, (\nabla f)^{-1}(x) \rangle + \langle h, u \rangle - f((\nabla f)^{-1}(x) + u) + f((\nabla f)^{-1}(x)) \\
&= \langle x, u \rangle + \langle h, (\nabla f)^{-1}(x) \rangle + \langle h, u \rangle - \langle \nabla f((\nabla f)^{-1}(x)), u \rangle + O(\|h\|^2) \\
&= \langle h, (\nabla f)^{-1}(x) \rangle + O(\|h\|^2).
\end{aligned}$$
Therefore $f^*$ is differentiable and $\nabla f^* = (\nabla f)^{-1}$. Combining this result with equation (25) proves that $f^*$ is $\frac{1}{\mu}$-smooth.

Conversely, if $f^*$ is $\frac{1}{\mu}$-smooth, then $\forall x, y \in \mathbb{R}^d$,

$$ \langle \nabla f^*(x) - \nabla f^*(y), x - y \rangle \geq \mu \| \nabla f^*(x) - \nabla f^*(y) \|^2. \quad (26) $$

By denoting $z = \nabla f^*(x)$ and $t = \nabla f^*(y)$, by Lemma 8(iii), we have $x = \nabla f(z)$ and $y = \nabla f(t)$, so

$$ \langle z - t, \nabla f(z) - \nabla f(t) \rangle \geq \mu \|z - t\|^2, $$

which implies that $f$ is $\mu$-strongly convex.

□

**Lemma 9 ([62]).** For $f$ convex and $x, y, z \in \mathbb{R}^d$, the two following statements are equivalent

$$ (i) \quad z = x + y \text{ and } f(x) + f^*(y) = \langle x, y \rangle. $$

$$ (ii) \quad x = \operatorname{Prox}_f(z) \text{ and } y = \operatorname{Prox}_{f^*}(z). $$

Lemma 9 shows the intrinsic link between the proximal operator and the convex conjugate. In particular, it shows that the inequality of Lemma 8(i) is attained only for couple that are of the form $(\operatorname{Prox}_f(z), \operatorname{Prox}_{f^*}(z))$. Moreover, it shows also the equality $\operatorname{Prox}_f(z) + \operatorname{Prox}_{f^*}(z) = z$.

*Proof. (i) ⇒ (ii)* For $u \in \mathbb{R}^d$, by the convex conjugate definition, we have

$$ f^*(y) \geq \langle u, y \rangle - f(u). $$

By using that $f(x) + f^*(y) = \langle x, y \rangle$, we get

$$ \langle x, y \rangle - f(x) \geq \langle u, y \rangle - f(u). $$

Then, using the previous inequality and $z = x + y$,

$$ \begin{aligned} \frac{1}{2} \|x - z\|^2 + f(x) &= \frac{1}{2} \|x\|^2 + \frac{1}{2} \|z\|^2 - \langle x, z \rangle + f(x) \\ &\leq \frac{1}{2} \|x\|^2 + \frac{1}{2} \|z\|^2 - \langle x, z \rangle + \langle x, y \rangle + f(u) - \langle u, y \rangle \\ &\leq \frac{1}{2} \|x\|^2 + \frac{1}{2} \|z\|^2 - \|x\|^2 + f(u) - \langle u, z \rangle + \langle u, x \rangle \\ &\leq \frac{1}{2} \|u - z\|^2 + f(u) - \frac{1}{2} \|x - u\|^2. \end{aligned} $$

Then, necessarily $x$ is the only minimum of the functional $u \mapsto \frac{1}{2}\|u - z\|^2 + f(u)$, which means that $x = \operatorname{Prox}_f(z)$. A similar computation gives that $y = \operatorname{Prox}_{f^*}(z)$.

(ii) ⇒ (i) We note $y' = z - x$, with $x = \operatorname{Prox}_f(z)$. By the convexity of $f$, for $t \in (0, 1)$ and $u \in \mathbb{R}^d$, we have

$$ f(tu + (1-t)x) \leq tf(u) + (1-t)f(x). $$

Then, because $x = \operatorname{Prox}_f(z)$, we have $\forall u \in \mathbb{R}^d$ and $\forall t \in (0, 1)$,

$$ \frac{1}{2} \|x - z\|^2 + f(x) \leq \frac{1}{2} \|tu + (1-t)x - z\|^2 + f(tu + (1-t)x). $$

By combining the two previous inequalities, we get

$$ \begin{aligned} \frac{1}{2} \|x - z\|^2 + f(x) &\leq \frac{1}{2} \|tu + (1-t)x - z\|^2 + tf(u) + (1-t)f(x) \\ &\quad - \frac{t^2}{2} \|u - x\|^2 + t\langle u - x, y' \rangle \leq tf(u) - tf(x). \end{aligned} $$

By dividing by $t$ and letting $t \to 0$, we obtain that for all $u \in \mathbb{R}^d$

$$ \begin{aligned} \langle u - x, y' \rangle &\leq f(u) - f(x) \\ \langle u, y' \rangle - f(u) &\leq \langle x, y' \rangle - f(x). \end{aligned} $$

So

$$ f^*(y') = \langle x, y' \rangle - f(x), $$

which ends the proof.

□
### D.5 Properties of the proximal operator

Based on the previous link between the convex conjugate and the Moreau envelope, we can now deduce some properties on the proximal operator (Corollary (2) and Lemma 11) of weakly convex functions. Moreover, we will study how the Moreau envelope preserves strong convexity (Lemma 10).

**Corollary 2.** For $g$ $\rho$-weakly convex with $\gamma\rho < 1$ and differentiable at $x \in \mathbb{R}^d$, we have

$$ \operatorname{Prox}_{\gamma g}(x + \gamma \nabla g(x)) = x $$

Corollary 2 leads to the notation $\operatorname{Prox}_{\gamma g} = (I_d + \gamma \nabla g)^{-1}$.

*Proof.* We first demonstrate the relation for convex functions $g$ and then use it for weakly convex ones.

* For $g$ convex and $x \in \mathbb{R}^d$, by Lemma 8(ii), we have for $y = \nabla g(x)$, $g(x) + g^*(y) = \langle x, y \rangle$. Then by Lemma 9, we get $x = \operatorname{Prox}_g (x + \nabla g(x))$. Then, if $g$ is convex, $\gamma g$ is convex and we get $\forall \gamma > 0$,

$$ \operatorname{Prox}_{\gamma g}(x + \gamma \nabla g(x)) = x. \qquad (27) $$

* If $g$ is $\rho$-weakly convex with $\gamma\rho < 1$, then we introduce the convex function $g_\rho = g + \frac{\rho}{2}\|\cdot\|^2$. By Lemma 5, we have

$$ \operatorname{Prox}_{\gamma g}(x) = \operatorname{Prox}_{\frac{\gamma}{1-\gamma\rho}g_{\rho}} \left( \frac{x}{1-\gamma\rho} \right). $$

Then

$$ \begin{aligned} \operatorname{Prox}_{\gamma g}(x + \gamma \nabla g(x)) &= \operatorname{Prox}_{\frac{\gamma}{1-\gamma\rho}g_{\rho}} \left( \frac{x + \gamma \nabla g(x)}{1 - \gamma\rho} \right) \\ &= \operatorname{Prox}_{\frac{\gamma}{1-\gamma\rho}g_{\rho}} \left( \frac{x + \gamma \nabla g_{\rho}(x) - \gamma\rho x}{1 - \gamma\rho} \right) \\ &= \operatorname{Prox}_{\frac{\gamma}{1-\gamma\rho}g_{\rho}} \left( x + \frac{\gamma}{1 - \gamma\rho} \nabla g_{\rho}(x) \right) \\ &= x. \end{aligned} $$

The last equality is obtained by applying (27) on the convex function $g_\rho$.

□

**Lemma 10.** For a function $g$ $\rho$-strongly convex, then $g^\gamma$ is $\frac{\rho}{1+\gamma\rho}$ strongly convex.

Note that $g^\gamma$ is less strongly convex than $g$ as $\frac{\rho}{1+\gamma\rho} < \rho$.

*Proof.* Assume that $g$ is $\rho$-strongly convex, then, thanks to Lemma 7(ii) and Lemma 7(iv), $(g^\gamma)^* = (g \Box \frac{1}{2\gamma} \|x\|^2)^* = g^* + \frac{\gamma}{2} \|x\|^2$ that is convex and $\frac{1}{\rho} + \gamma = \frac{1+\rho\gamma}{\rho}$-smooth thanks to Lemma 8(iv). Combining Lemma 8(iv) and Lemma 8(iii), we deduce that $(g^\gamma)^{**} = g^\gamma$ is $\frac{\rho}{1+\rho\gamma}$ strongly convex. □

**Lemma 11.** For $g$ $\rho$-weakly convex with $\gamma\rho < 1$, $\operatorname{Prox}_{\gamma g}$ is $\frac{1}{1-\gamma\rho}$ Lipschitz.

Note that in particular, if $g$ is convex, then $\rho = 0$ and $\forall \gamma > 0$, $\operatorname{Prox}_{\gamma g}$ is 1-Lipschitz. The fact that the proximal operator is Lipschitz is important for guaranteeing its stability.

*Proof.* The general proof of this Lemma 11 can be found in [40, Proposition 2] where the authors used [62, Proposition 5.b]. We provide here a more direct proof in the case of a differentiable function $g$.
As $g$ is $\rho$-weakly convex and differentiable, then $\forall x, y \in \mathbb{R}^d$, we have by the convexity of $g + \frac{\rho}{2} \| \cdot \|^2$

$$g(x) + \frac{\rho}{2} \|x\|^2 + \langle \nabla g(x) + \rho x, y - x \rangle \leq g(y) + \frac{\rho}{2} \|y\|^2 \\ \langle \nabla g(x), y - x \rangle \leq g(y) - g(x) + \frac{\rho}{2} \|x - y\|^2.$$

For $x, x' \in \mathbb{R}^d$, we denote $z = \text{Prox}_{\gamma g}(x)$ and $z' = \text{Prox}_{\gamma g}(x')$, we get

$$\langle \nabla g(z), z' - z \rangle \leq g(z') - g(z) + \frac{\rho}{2} \|z - z'\|^2 \\ \langle \nabla g(z'), z - z' \rangle \leq g(z) - g(z') + \frac{\rho}{2} \|z - z'\|^2.$$

By the optimal condition of the proximal operator, we also have

$$\frac{1}{\gamma}(z - x) + \nabla g(z) = 0 \\ \frac{1}{\gamma}(z' - x') + \nabla g(z') = 0,$$

Combining the two previous sets of equations, we obtain

$$\frac{1}{\gamma} \langle x - z, z' - z \rangle \leq g(z') - g(z) + \frac{\rho}{2} \|z - z'\|^2 \\ \frac{1}{\gamma} \langle x' - z', z - z' \rangle \leq g(z) - g(z') + \frac{\rho}{2} \|z - z'\|^2.$$

By summing the two previous inequalities and using the Cauchy-Schwarz inequality, we get

$$-\rho \|z - z'\|^2 \leq \frac{1}{\gamma} \langle z - z' - x + x', z' - z \rangle \\ (1 - \gamma \rho) \|z - z'\|^2 \leq \langle x' - x, z' - z \rangle \\ (1 - \gamma \rho) \|z - z'\|^2 \leq \|x - x'\| \|z - z'\|.$$

Therefore, we get that, $\forall x, x' \in \mathbb{R}^d$

$$\| \text{Prox}_{\gamma g}(x) - \text{Prox}_{\gamma g}(x') \| \leq \frac{1}{1 - \gamma \rho} \|x - x'\|,$$

which proves that $\text{Prox}_{\gamma g}$ is $\frac{1}{1-\gamma\rho}$-Lipschitz.

<figure><img src="image_9.png" /></figure>

### D.6 Proximal operator as a gradient step

Now, we can prove the main results of this appendix, that allows to interpret the proximal operator as a gradient descent step on the Moreau envelope. Note that Lemma 12 also proves that the Moreau envelope is differentiable even if $g$ is not.

**Lemma 12.** For a function $g$ $\rho$-weakly convex, for $\gamma\rho < 1$, we have

$$\nabla g^{\gamma}(x) = \frac{1}{\gamma}(x - \text{Prox}_{\gamma g}(x))$$

*Proof.* Our proof is based on the strategy detailed in [76, Theorem 2.6] for $g$ convex, that we generalize for $g$ $\rho$-weakly convex. This Lemma is formulated in the weakly-convex setting in [6, Lemma 2.2] and in [20, Lemma 2.2] without proofs.

For $x \in \mathbb{R}^d$, we introduce $v = \text{Prox}_{\gamma g}(x)$, $w = \frac{1}{\gamma}(x - v)$ and $h(u) = g^{\gamma}(x + u) - g^{\gamma}(x) - \langle w, u \rangle$. We aim to prove that $g^{\gamma}$ is differentiable at $x$ with $\nabla g^{\gamma} = w$, which is equivalent to $h$ differentiable at 0 and $\nabla h(0) = 0$.
By definition of the Moreau envelope, since $v = \text{Prox}_{\gamma g}(x)$, we have $g^{\gamma}(x) = \frac{1}{2\gamma}\|x - v\|^2 + g(v)$ and $g^{\gamma}(x + u) \le \frac{1}{2\gamma}\|x + u - v\|^2 + g(v)$. So, we get

$$
\begin{align*}
h(u) &\le \frac{1}{2\gamma} \|x + u - v\|^2 - \frac{1}{2\gamma} \|x - v\|^2 - \langle w, u \rangle \\
&\le \frac{1}{2\gamma} \|u\|^2 + \frac{1}{\gamma} \langle x - v, u \rangle - \langle w, u \rangle \\
&\le \frac{1}{2\gamma} \|u\|^2.
\end{align*}
$$

Moreover, as $g^{\gamma}(x+u)$ is $\frac{\rho}{1-\rho\gamma}$-weakly convex by Lemma 6, then $h(u) + \frac{\rho}{2(1-\rho\gamma)}\|u\|^2 = g^{\gamma}(x+u) + \frac{\rho}{2(1-\rho\gamma)}\|u\|^2 - g^{\gamma}(x) - \langle w, u \rangle$ is convex in $u$. Therefore

$$
\begin{align*}
\frac{1}{2} \left( h(u) + \frac{\rho}{2(1-\rho\gamma)} \|u\|^2 + h(-u) + \frac{\rho}{2(1-\rho\gamma)} \|u\|^2 \right) &\geq h(0) = 0 \\
h(u) &\geq -\frac{\rho}{(1-\rho\gamma)} \|u\|^2 - h(-u) \geq \left( -\frac{\rho}{1-\rho\gamma} - \frac{1}{2\gamma} \right) \|u\|^2.
\end{align*}
$$

So, we get

$$
\left(-\frac{\rho}{1-\rho\gamma} - \frac{1}{2\gamma}\right) \|u\|^2 \le h(u) \le \frac{1}{2\gamma} \|u\|^2.
$$

Therefore $h$ is differentiable at 0 and $\nabla h(0) = 0$. This concludes the proof. $\square$

<figure><img src="image_10.png" /></figure>

We can deduce that the Moreau envelope is smooth, i.e. $\nabla g^{\gamma}$ is Lipschitz.

**Lemma 13.** If $g$ is $\rho$-weakly convex, with $\gamma\rho < 1$, then $g^{\gamma}$ is $\max\left(\frac{1}{\gamma}, \frac{\rho}{1-\gamma\rho}\right)$ smooth.

*Proof.* By Lemma 6, $g^{\gamma}$ is $rac{\rho}{1-\gamma\rho}$-weakly convex. Moreover, for $f$ a convex function, by [3, Proposition 23.8], $x \mapsto x - \text{Prox}_{\gamma f}(x)$ is 1-Lipschitz. We prove this fact for $f$ convex and differentiable. For $x, y \in \mathbb{R}^d$ and $u = \text{Prox}_{\gamma f}(x)$, $v = \text{Prox}_{\gamma f}(y)$, we get by the optimal condition of the proximal operator

$$
\langle (x - u) - (y - v), u - v \rangle = \frac{1}{\gamma} \langle \nabla f(u) - \nabla f(v), u - v \rangle \geq 0,
$$

by convexity of $f$. Then

$$
\begin{align*}
\|x - y\|^2 &= \|(x - u) - (y - v) + u - v\|^2 \\
&= \|(x - u) - (y - v)\|^2 + \|u - v\|^2 + 2\langle (x - u) - (y - v), u - v \rangle \\
&\geq \|(x - u) - (y - v)\|^2 + \|u - v\|^2.
\end{align*}
$$

So $\|(x-u)-(y-v)\| \le \|x-y\|$, which means that $x \mapsto x - \text{Prox}_{\gamma f}(x)$ is 1-Lipschitz.

Combining the previous property with Lemma 12, we get that $\nabla f^{\gamma}$ is $\frac{1}{\gamma}$-Lipschitz. So, we have the inequality,

$$
\langle \nabla f^{\gamma}(x) - \nabla f^{\gamma}(y), x - y \rangle \le \frac{1}{\gamma} \|x - y\|^2. \quad (28)
$$

However, Lemma 5 gives that

$$
g^{\gamma} = g_{\rho}^{\frac{\gamma}{1-\rho\gamma}} \left( \frac{x}{1-\rho\gamma} \right) - \frac{\rho}{2(1-\rho\gamma)} \|x\|^2.
$$
As $g_\rho$ is convex, we can apply equation (28) and obtain

$$
\begin{align*}
\langle \nabla g^{\gamma}(x) - \nabla g^{\gamma}(y), x - y \rangle \\
&= \frac{1}{1 - \rho\gamma} \langle \nabla g_{\rho}^{\frac{1}{1-\gamma\rho}} \left( \frac{x}{1-\gamma\rho} \right) - \nabla g_{\rho}^{\frac{1}{1-\gamma\rho}} \left( \frac{y}{1-\gamma\rho} \right), x - y \rangle - \frac{\rho}{1-\gamma\rho} \|x - y\|^2 \\
&\le \frac{1}{1-\rho\gamma} \frac{1}{\gamma} \|x - y\|^2 - \frac{\rho}{1-\gamma\rho} \|x - y\|^2 \\
&\le \frac{1}{\gamma} \|x - y\|^2.
\end{align*}
$$

With the weak convexity of $g^\gamma$ (Lemma 6), we get

$$
- \frac{\rho}{1 - \gamma \rho} \|x - y\|^2 \leq \langle \nabla g^{\gamma}(x) - \nabla g^{\gamma}(y), x - y \rangle \leq \frac{1}{\gamma} \|x - y\|^2.
$$

So, we get that $\nabla g^{\gamma}$ is max $\left(\frac{1}{\gamma}, \frac{\rho}{1-\gamma\rho}\right)$ smooth. $\square$

D.7 Second derivative of the Moreau envelope and convexity of the image of the proximal operator

We can now prove that the Moreau envelope is twice differentiable if $g$ is twice differentiable and Lipschitz on the image of the proximal operator. Moreover, we will also show that the border of the the convex envelop of $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ is of measure zero (Proposition 2). In other words, the image of the proximity operator of a weakly convex function is almost convex.

**Lemma 14.** Let $g$ be a $\rho$-weakly convex function. Assume that $g$ that is $C^2$ and $L_g$-smooth on $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ with $\rho\gamma < 1$ and $L_g\gamma < 1$. Then we get, for $x \in \mathbb{R}^d$

$$
\begin{align*}
\nabla \text{Prox}_{\gamma g}(x) &= (I_d + \gamma \nabla^2 g(\text{Prox}_{\gamma g}(x)))^{-1} \\
\nabla^2 g^{\gamma}(x) &= \frac{1}{\gamma} (I_d - (I_d + \gamma \nabla^2 g(\text{Prox}_{\gamma g}(x)))^{-1}).
\end{align*}
$$

*Proof.* This proof is an adaptation of [65, Theorem 3.4] for weakly convex functions.

By the optimal condition of the proximal operator, we get that

$$
\frac{1}{\gamma}(\text{Prox}_{\gamma g}(x) - x) + \nabla g(\text{Prox}_{\gamma g}(x)) = 0.
$$

Combined with Lemma 12, we get that

$$
\nabla g^{\gamma}(x) = \nabla g(\text{Prox}_{\gamma g}(x)), \quad (29)
$$

Then for all $x, h \in \mathbb{R}^d$,

$$
\begin{align*}
\nabla g^{\gamma}(x+h) - \nabla g^{\gamma}(x) &= \nabla g(\text{Prox}_{\gamma g}(x+h)) - \nabla g(\text{Prox}_{\gamma g}(x)) \\
&= \nabla^2 g(\text{Prox}_{\gamma g}(x)) (\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)) + o(\|\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)\|).
\end{align*}
$$

However, by Lemma 11, $\text{Prox}_{\gamma g}$ is $\frac{1}{1-\gamma\rho}$-Lipschitz. Therefore $o(\|\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)\|) = o(\|h\|)$. Then, we get

$$
\nabla g^{\gamma}(x+h) - \nabla g^{\gamma}(x) = \nabla^2 g(\text{Prox}_{\gamma g}(x)) (\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)) + o(\|h\|).
$$

Injecting the characterization from Lemma 12 in the left part of the previous equation, we get

$$
\begin{align*}
& \frac{1}{\gamma}h - \frac{1}{\gamma}(\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)) \\
&= \nabla^2 g(\text{Prox}_{\gamma g}(x)) (\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)) + o(\|h\|) \\
& \frac{1}{\gamma}h = \left(\frac{1}{\gamma}I_d + \nabla^2 g(\text{Prox}_{\gamma g}(x))\right) (\text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x)) + o(\|h\|) \\
& \text{Prox}_{\gamma g}(x+h) - \text{Prox}_{\gamma g}(x) = \frac{1}{\gamma}\left(\frac{1}{\gamma}I_d + \nabla^2 g(\text{Prox}_{\gamma g}(x))\right)^{-1}(h) + o(\|h\|),
\end{align*}
$$
where $\frac{1}{\gamma}I_d + \nabla^2 g(\text{Prox}_{\gamma g}(x))$ is invertible because $g$ is $L_g$-smooth at the point $\text{Prox}_{\gamma g}(x) \in \text{Prox}_{\gamma g}(\mathbb{R}^d)$ with $\gamma L_g < 1$. The last equation shows that $\text{Prox}_{\gamma g}$ is differentiable and its gradient is

$$ \nabla \text{Prox}_{\gamma g} = (I_d + \gamma \nabla^2 g \circ \text{Prox}_{\gamma g})^{-1}. $$

Then, we differentiate Equation (29) to get

$$ \begin{aligned} \nabla^2 g^{\gamma}(x) &= \nabla^2 g(\text{Prox}_{\gamma g}(x)) \nabla \text{Prox}_{\gamma g}(x) \\ &= \nabla^2 g(\text{Prox}_{\gamma g}(x)) (I_d + \gamma \nabla^2 g(\text{Prox}_{\gamma g}(x)))^{-1} \\ &= \frac{1}{\gamma} (I_d - (I_d + \gamma \nabla^2 g(\text{Prox}_{\gamma g}(x)))^{-1}), \end{aligned} $$

which shows the second part of Lemma 14. □

We can now deduce that for $\gamma$ sufficiently small the Moreau envelope is smooth (with a constant independent of $\gamma$) on all $\mathbb{R}^d$ if the function $g$ is smooth on $\text{Prox}_{\gamma g}(\mathbb{R}^d)$.

**Lemma 15.** Let $g$ be a $\rho$-weakly convex function. Assume that $g$ is $C^2$ and $L_g$-smooth on $\text{Prox}_{\gamma g}(\mathbb{R}^d)$ with $L_g\gamma \le \frac{1}{2}$ and $\rho\gamma < 1$. Then $g^\gamma$ is $2L_g$-smooth on $\mathbb{R}^d$.

*Proof.* By Lemma 14, we have

$$ \nabla^2 g^{\gamma}(x) = \frac{1}{\gamma} \left( I_d - \left( I_d + \gamma \nabla^2 g(\text{Prox}_{\gamma g}(x)) \right)^{-1} \right). $$

As we assume that $-L_g I_d \le \nabla^2 g(\text{Prox}_{\gamma g}(x)) \le L_g I_d$, we get

$$ \frac{1}{\gamma} \left( 1 - \frac{1}{1 - L_g \gamma} \right) I_d \le \nabla^2 g^{\gamma}(x) \le \frac{1}{\gamma} \left( 1 - \frac{1}{1 + L_g \gamma} \right) I_d \\ \underbrace{- \frac{L_g}{1 - L_g \gamma}}_{:=u(\gamma)} I_d \le \nabla^2 g^{\gamma}(x) \le \underbrace{\frac{L_g}{1 + L_g \gamma}}_{:=v(\gamma)} I_d $$

Since $u'(\gamma) = -\frac{L_g^2}{(1-L_g\gamma)^2} \le 0$, $u$ is decreasing. As $u(\frac{1}{2L_g}) = -2L_g$, we obtain that $u(\gamma) \ge -2L_g$, for all $\gamma \in [0, \frac{1}{2L_g}]$. In the same way, $v'(\gamma) = -\frac{L_g^2}{(1+L_g\gamma)^2} \le 0$, so $v$ is decreasing. Since $v(0) = L_g$, for $\gamma \ge 0$, we have $v(\gamma) \le L_g$. Finally, for $\gamma \in [0, \frac{1}{2L_g}]$, we get

$$ -2L_g I_d \le \nabla^2 g^{\gamma}(x) \le L_g I_d. $$

So $g^\gamma$ is $2L_g$-smooth on $\mathbb{R}^d$. □

**Proposition 2.** Let $g$ be a $\rho$-weakly convex function with $\rho\gamma < 1$. We have

$$ \text{Leb}(\text{dom}(g) \setminus \text{Prox}_{\gamma g}(\mathbb{R}^d)) = 0, \qquad (30) $$

with $\text{Leb}$ the Lebesgue measure, $\text{dom}(g) = \{x \in \mathbb{R}^d | g(x) < +\infty\}$. Moreover, we have that

$$ \text{Leb}(\text{Conv}(\text{Prox}_{\gamma g}(\mathbb{R}^d)) \setminus \text{Prox}_{\gamma g}(\mathbb{R}^d)) = 0, \qquad (31) $$

where $\text{Conv}(\text{Prox}_{\gamma g}(\mathbb{R}^d))$ is the convex envelop of $\text{Prox}_{\gamma g}(\mathbb{R}^d)$.

*Proof.* We denote by $\text{dom}(g)$ the set where $g < +\infty$. By the weak convexity of $g$, we know that $\text{dom}(g)$ is convex. By [54, Theorem 1], because $\text{dom}(g)$ is convex, we have that $\text{Leb}(\text{dom}(g) \setminus \text{int}(\text{dom}(g))) = 0$.

By [43, Theorem 4.2.3], a convex function is differentiable almost everywhere on the interior of its domain. As $g + \frac{\ell}{2}\|\cdot\|^2$ is convex, it is thus differentiable almost everywhere on the interior of its domain. As a consequence $g$ is differentiable almost everywhere on the interior of its domain,
<figure><img src="image_11.png" /><figcaption>Figure 7: Representation of the function $f$ defined in Equation (32).</figcaption></figure>

$\text{int}(\text{dom}(g))$. Denoting as $\Gamma \subset \text{int}(\text{dom}(g))$ the subset on which $g$ is differentiable, we thus have $\text{Leb}(\text{int}(\text{dom}(g)) \setminus \Gamma) = 0$.

Moreover, for $x \in \Gamma$, $g$ is differentiable at $x$, so by Corollary 2, we have $x = \text{Prox}_{\gamma g}(x + \gamma \nabla g(x)) \in \text{Prox}_{\gamma g}(\mathbb{R}^d)$. So we get $\Gamma \subset \text{Prox}_{\gamma g}(\mathbb{R}^d)$. Therefore, we can deduce that $\text{Leb}(\text{int}(\text{dom}(g)) \setminus \text{Prox}_{\gamma g}(\mathbb{R}^d)) = 0$.

Combining the fact that $\text{Leb}(\text{dom}(g) \setminus \text{int}(\text{dom}(g))) = 0$ and $\text{Leb}(\text{int}(\text{dom}(g)) \setminus \text{Prox}_{\gamma g}(\mathbb{R}^d)) = 0$, we get that $\text{Leb}(\text{dom}(g) \setminus \text{Prox}_{\gamma g}(\mathbb{R}^d)) = 0$. It shows equation (30).

Moreover for $x \in \mathbb{R}^d$ and $y \in \text{dom}(g)$, by definition of the proximal operator, we get $g(\text{Prox}_{\gamma g}(x)) + \frac{1}{2\gamma}\|x - \text{Prox}_{\gamma g}(x)\|^2 \le g(y) + \frac{1}{2\gamma}\|x - y\|^2$. So $g(\text{Prox}_{\gamma g}(x)) < +\infty$, which gives that $\text{Prox}_{\gamma g}(\mathbb{R}^d) \subset \text{dom}(g)$. Because $\text{dom}(g)$ is convex, we get that $\text{Conv}(\text{Prox}_{\gamma g}(\mathbb{R}^d)) \subset \text{dom}(g)$. Combining with equation (30), it gives equation (31). $\square$

## E On the strong convexity on non-convex sets

**Definition 4.** For an open subset $\mathbf{S} \subset \mathbb{R}^d$ and $\mu \ge 0$, a function $f \in C^2(\mathbf{S}, \mathbb{R})$ is said to be $\mu$-strongly convex on $\mathbf{S}$ if and only if $\forall x \in \mathbf{S}, \nabla^2 f \succeq \mu I_d$.

A function that is $\mu$-strongly convex on $\mathbf{S}$ with $\mu = 0$ is called convex on $\mathbf{S}$.

**Proposition 3.** For a convex open set $\mathbf{S}$ and a function $f \in C^2(\mathbf{S}, \mathbb{R})$, $f$ is convex on $\mathbf{S}$ if and only if $\forall t \in [0, 1], \forall x, y \in \mathbf{S}, f(tx + (1-t)y) \le tf(x) + (1-t)f(y)$.

Proposition 3 is standard and we let the proof to the reader.

**Remark 1.** If $\mathbf{S}$ is non-convex, then even if $f$ is convex on $\mathbf{S}$, the inequality $\forall t \in [0, 1], \forall x, y \in \mathbf{S}, f(tx + (1-t)y) \le tf(x) + (1-t)f(y)$ is not necessarily verified.

This is for instance the case with $\mathbf{S} = \{z \in \mathbb{R}^d \| \|z\| < 1\} \cup \{z \in \mathbb{R}^d \| \|z\| > 2\}$ and $f$ defined by

$$f(x) = 1 \quad \text{if} \quad \|z\| < 1 \qquad (32)$$

$$f(x) = 0 \quad \text{if} \quad \|z\| > 2, \qquad (33)$$

see Figure 7 for a representation of $f$. Then, $\forall x \in \mathbf{S}, \nabla^2 f = 0$ so $f$ is convex on $\mathbf{S}$. However, for $u \in \mathbb{R}^d$ such that $\|u\| = 1$, we have

$$1 = f(0) = f\left(\frac{1}{2}3u + \frac{1}{2}(-3u)\right) > \frac{1}{2}f(3u) + \frac{1}{2}f(-3u) = 0.$$

**Lemma 16.** If $g$ is $L_g$-smooth on $\mathbb{R}^d$ with $\gamma L_g < 1$ and there exists $R > 0$, such that $g$ is $\mu$-strongly convex on $\mathbb{R}^d \setminus B(0, R)$, i.e. $\nabla^2 g(x) \succeq \mu I_d$, then there exists $R_0 > 0$, such that $\forall \gamma \in [0, \frac{1}{2L_g}]$, $g^\gamma$ is $\frac{\mu}{1+\mu\gamma}$-strongly convex on $\mathbb{R}^d \setminus B(0, R_0)$.
*Proof.* As $g$ is strongly convex on $\mathbb{R}^d \setminus B(0, R)$, it is coercive, so $g$ is lower bounded on $\mathbb{R}^d$, i.e. there exists $g_{inf} = \min_{x \in \mathbb{R}^d} g(x) > -\infty$ such that $\forall x \in \mathbb{R}^d$, $g(x) \ge g_{inf}$. We recall that the Moreau envelope of $g$ is defined as

$$g^{\gamma}(x) = \inf_{y \in \mathbb{R}^d} \frac{1}{2\gamma} \|x - y\|^2 + g(y). \qquad (34)$$

For $x \in \mathbb{R}^d$, such that $\|x\| \ge R + 1$, for $y \in \mathbb{R}^d$, such that $\|y - x\| \ge \|x\| - R - 1$, we have

$$\frac{1}{2\gamma} \|x - y\|^2 + g(y) \ge \frac{1}{2\gamma} (\|x\| - R - 1)^2 + g_{inf}. \qquad (35)$$

Since $\nabla g$ is $L_g$-Lipschitz, we get that

$$g(x) \le g(0) + \langle \nabla g(0), x \rangle + \frac{L_g}{2} \|x\|^2. \qquad (36)$$

Due to $L_g < \frac{1}{\gamma}$, we get that asymptotically $\frac{1}{2\gamma} (\|x\| - R - 1)^2 + g_{inf}$ is larger than $g(0) + \langle \nabla g(0), x \rangle + \frac{L_g}{2} \|x\|^2$. More formally, we study the quadratic form

$$\begin{aligned} & \frac{1}{2\gamma} (\|x\| - R - 1)^2 + g_{inf} - \left( g(0) + \langle \nabla g(0), x \rangle + \frac{L_g}{2} \|x\|^2 \right) \\ & = \frac{1}{2} \left( \frac{1}{\gamma} - L_g \right) \|x\|^2 - \frac{R+1}{\gamma} \|x\| - \langle \nabla g(0), x \rangle + g_{inf} - g(0) + \frac{(R+1)^2}{2\gamma} \\ & \ge \frac{1}{2} \left( \frac{1}{\gamma} - L_g \right) \|x\|^2 - \left( \frac{R+1}{\gamma} + \|\nabla g(0)\| \right) \|x\| + g_{inf} - g(0) + \frac{(R+1)^2}{2\gamma} \\ & \ge \left( \frac{1}{4} \left( \frac{1}{\gamma} - L_g \right) \|x\| - \left( \frac{R+1}{\gamma} + \|\nabla g(0)\| \right) \right) \|x\| \\ & \quad + \frac{1}{4} \left( \frac{1}{\gamma} - L_g \right) \|x\|^2 + g_{inf} - g(0) + \frac{(R+1)^2}{2\gamma}. \end{aligned}$$

Then, the previous quadratic form is non negative if $x \in \mathbb{R}^d$ verifies

$$\|x\| \ge 4 \frac{\frac{R+1}{\gamma} + \|\nabla g(0)\|}{\frac{1}{\gamma} - L_g} \\ \text{and } \|x\| \ge 2 \sqrt{\frac{g_{inf} - g(0) + \frac{(R+1)^2}{2\gamma}}{\frac{1}{\gamma} - L_g}}.$$

By defining

$$r_{\gamma} = \max \left( 4 \frac{\frac{R+1}{\gamma} + \|\nabla g(0)\|}{\frac{1}{\gamma} - L_g}, 2 \sqrt{\frac{g_{inf} - g(0) + \frac{(R+1)^2}{2\gamma}}{\frac{1}{\gamma} - L_g}} \right),$$

we get that $\forall x \in \mathbb{R}^d \setminus B(0, r_{\gamma})$,

$$g(0) + \langle \nabla g(0), x \rangle + \frac{L_g}{2} \|x\|^2 \le \frac{1}{2\gamma} (\|x\| - R - 1)^2 + g_{inf}.$$

Moreover, for $\gamma \le \frac{1}{2L_g}$, $\frac{1}{\gamma} - L_g \ge L_g$ and $1 - \gamma L_g \ge \frac{1}{2}$. So we get that

$$r_{\gamma} \le R_0 = \max \left( 8(R+1) + \frac{4\|\nabla g(0)\|}{L_g}, 2\sqrt{\frac{g_{inf} - g(0)}{L_g} + (R+1)^2} \right).$$

Therefore, we get from (36) that $\forall \gamma \in [0, \frac{1}{2L_g}]$ and $\forall x \in \mathbb{R}^d \setminus B(0, R_0)$,

$$g(x) \le g(0) + \langle \nabla g(0), x \rangle + \frac{L_g}{2} \|x\|^2 \le \frac{1}{2\gamma} (\|x\| - R - 1)^2 + g_{inf}.$$
By combining the previous inequality and Equation (35), we get that for $x \in \mathbb{R}^d \setminus B(0, R_0)$ and for $y$ such that $\|y - x\| \ge \|x\| - R - 1$, we have

$$g(x) \le \frac{1}{2\gamma} \|x - y\|^2 + g(y).$$

Hence for $x \in \mathbb{R}^d \setminus B(0, R_0)$, the infimum in the computation of $g^\gamma$ (see relation 34) is reached inside the ball $B(x, \|x\| - R - 1)$, i.e.

$$g^\gamma(x) = \inf_{y \in \mathbb{R}^d, \|y-x\| \le \|x\|-R-1} \frac{1}{2\gamma} \|x-y\|^2 + g(y).$$

Moreover, $B(x, \|x\| - R - 1) \subset \mathbb{R}^d \setminus B(0, R)$ and $B(x, \|x\| - R - 1)$ is a convex set. As we assumed that $g$ is $\mu$-strongly convex on $\mathbb{R}^d \setminus B(0, R)$, then $(x, y) \in \mathbb{R}^d \times B(x, \|x\| - R - 1) \mapsto \frac{1}{2\gamma} \|x-y\|^2 + g(y) - \frac{\mu}{2(1+\gamma\mu)} \|x\|^2$ is convex. Then, by Lemma 3, $g^\gamma - \frac{\mu}{2(1+\gamma\mu)} \|x\|^2$ is convex on $\mathbb{R}^d \setminus B(0, R_0)$ which ends the proof. $\square$

# F Proof of Section 2

## F.1 Introduction to Markov Process and Markov Chain: definition and concepts

In this part, we introduce the tools of Markov Chain and Markov processes that will be useful for our analysis. We refer to the book of [25] for a detailed introduction.

We define a continuous-time Markov stochastic differential equation by

$$dX_t = b(X_t)dt + \sqrt{2}dw_t \quad (37)$$

$$X_0 \sim \mu, \quad (38)$$

where $dw_t$ is a Wiener process and $\mu$ a probability law on $\mathbb{R}^d$. The deterministic part of the equation $b(X_t)$ is called the drift. If the drift $b$ is Lipschitz, then the previous equation has a strong solution, i.e. there is a unique solution $(X_t)_{t \ge 0}$ up to negligible fluctuation. In this paper, we will always be in the case $b$ Lipschitz which ensures that the previous equation is well defined and defines a unique $X_t$.

The constant standard deviation of $\sqrt{2}$ ensures that if $b(X_t) = -\nabla f(X_t)$ and $b$ is Lipschitz, then the process converges in law to $\pi \propto e^{-f}$ [75]. The following tools are not restricted to the case of this constant standard deviation case. In the case of $\mu = \delta_x$, we say that the process $X_t$ starts on $x$.

We will study the law of $X_t$ for each time $t \ge 0$. To explore these laws, we introduce the Markov semi group $P_t$ defined for $x \in \mathbb{R}^d$ and $A$ a measurable set by

$$P_t(x, A) = \mathbb{P}[X_t \in A | X_0 = x] \quad (39)$$

$P_t(x, A)$ quantifies the probability that the process starting at the point $x$ and following the stochastic equation (37) is in the set $A$ at time $t$. It is called a semi-group because if we define $(P_t P_s)(x, A) = \int_{\mathbb{R}^d} P_t(y, A) P_s(x, dy)$, then we have $\forall t, s \ge 0$, $P_t P_s = P_{t+s}$. We introduce two operations with the object $P_t$. For a distribution $\mu$ on $\mathbb{R}^d$, we define $\mu P_t(A) := \int_{\mathbb{R}^d} P_t(x, A) d\mu(x)$. Thus $\mu P_t$ is the distribution $X_t$ on $\mathbb{R}^d$ knowing that $X_0$ follow $\mu$. In particular $P_0(x, A) = \delta_x(A)$ and $\mu P_0 = \mu$. For a function $\phi: \mathbb{R}^d \to \mathbb{R}$, we define $P_t \phi(x) := \int_{\mathbb{R}^d} \phi(y) P_t(x, dy) = \mathbb{E}(\phi(X_t)|X_0 = x)$. Thus $P_t \phi$ is a function of $\mathbb{R}^d$ to $\mathbb{R}$. In particular $P_0 \phi = \phi$. The two previous operations can be combined together to get $\mu P_t \phi = \mathbb{E}(\phi(X_t)|X_0 \sim \mu)$.

With a step-size $\gamma > 0$, the Markov process (37) can be discretized with the Euler-Maruyama scheme leading to

$$X_{k+1} = X_k + \gamma b(X_k) + \sqrt{2\gamma} Z_{k+1} \quad (40)$$

$$X_0 \sim \mu, \quad (41)$$

with $Z_{k+1} \sim \mathcal{N}(0, I_d)$. This process is a Markov Chain. The deterministic part $b(X_k)$ is called the drift. Similarly to the Markov semi-group, we introduce the Markov kernel $R_\gamma$ defined for $x \in \mathbb{R}^d$ and $A \subset \mathbb{R}^d$ by $R_\gamma(x, A) = \mathbb{P}[X_1 \in A | X_0 = x]$. We also define

$$R_\gamma^k(x, A) = \mathbb{P}[X_k \in A | X_0 = x]. \quad (42)$$
This Markov kernel $R_{\gamma}^{k}(x, A)$ quantifies the probability of being in a set $A$ at step $k$ for the Markov Chain starting at point $x$. We have a semi-group property in the sense that $\forall n, m \ge 0$, $R_{\gamma}^{n}R_{\gamma}^{m} = R_{\gamma}^{n+m}$. It can be seen as a transition matrix if the state space is not $\mathbb{R}^d$ but a finite space. In fact, for a probability distribution $\mu$, we define $\mu R_{\gamma}^{k}(A) = \int_{\mathbb{R}^d} R_{\gamma}^{k}(x, A)d\mu(x)$ the probability distribution of the discrete-time process $X_k$ at time $k$ knowing that $X_0 \sim \mu$. In particular $R_{\gamma}^{0}(x, A) = \delta_x(A)$ and $\mu R_{\gamma}^{0} = \mu$.

We can also define the application of the Markov kernel on a function $\phi: \mathbb{R}^d \to \mathbb{R}$ by $R_{\gamma}^{k}\phi(x) = \int_{\mathbb{R}^d} \phi(y)R_{\gamma}^{k}(x, dy)$. In particular $R_{\gamma}^{0}\phi = \phi$.

The Markov kernel will be the core object to study the evolution of the law of the iterates $X_k$.

For a function $V: \mathbb{R}^d \to [1, +\infty)$, we introduce the $V$-norm between two distributions $\mu, \nu$ defined on $\mathbb{R}^d$ by

$$ \|\mu - \nu\|_V := \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ |\phi| \le 1}} \int_{\mathbb{R}^d} \phi(d\mu - d\nu), \qquad (43) $$

The supremum is taken on $f$ measurable. Note that the function $V$ that defines the $V$-norm has nothing to do with the potential $V$ of the law $\pi \propto e^{-V}$ that we aim to sample. We keep the notation $V$-norm as it is standard in the literature [25, 56]. In particular for $V = 1$, it defines the total variation norm defined by

$$ \|\mu - \nu\|_{TV} := \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ |\phi| \le 1}} \int_{\mathbb{R}^d} \phi(d\mu - d\nu). \qquad (44) $$

Note that, for any $V \ge 1$, we have $\|\mu - \nu\|_{TV} \le \|\mu - \nu\|_V$.

## F.2 Technical lemmas

In this section we present technical lemmas that link the $V$-norm to the $W_p$ distance or the $TV$-distance and an application of Girsanov's theorem that will be useful for our analysis.

**Lemma 17.** [25, Theorem 19.1.7] For $V: \mathbb{R}^d \to [1, +\infty)$, we have

$$ \|\mu - \nu\|_V = \inf_{(X,Y) \sim \Pi} \mathbb{E}((V(X) + V(Y))1_{X \neq Y}), $$

with the infimum taken on the set $\Pi$ of coupling $\beta$ between $\mu$ and $\nu$. A coupling between $\mu$ and $\nu$ is a probability law $\beta$ on $\mathbb{R}^d \times \mathbb{R}^d$ with marginals $\mu$ and $\nu$.

**Lemma 18.** For $p \ge 1$ and $\mu, \nu$ two probability measures defined on $(\mathbb{R}^d, \mathcal{B}(\mathbb{R}^d))$, we have

$$ \mathbf{W}_p^p(\mu, \nu) \le 2^{p-1} \|\mu - \nu\|_{V_p}, $$

with $V_p = 1 + \| \cdot \|_p$.

*Proof.* By definition of the Wasserstein distance, we have

$$ \mathbf{W}_p^p(\mu, \nu) = \inf_{(X,Y) \sim \beta} \mathbb{E} (\|X - Y\|^p), $$

where the inf is taken on the set of coupling $\beta$ between $\mu$ and $\nu$. Using that the function $x \mapsto x^p$ is convex, we have that $(x + y)^p \le 2^{p-1}(x^p + y^p)$ and Lemma 17. This leads to

$$ \begin{aligned} \mathbf{W}_p^p(\mu, \nu) &= \inf_{(X,Y) \sim \beta} \mathbb{E}((\|X - Y\|^p)1_{X \neq Y}) \\ &\le 2^{p-1} \inf_{(X,Y) \sim \beta} \mathbb{E}((\|X\|^p + \|Y\|^p)1_{X \neq Y}) \\ &\le 2^{p-1} \inf_{(X,Y) \sim \beta} \mathbb{E}((1 + \|X\|^p + 1 + \|Y\|^p)1_{X \neq Y}) \\ &= 2^{p-1} \|\mu - \nu\|_{V_p}, \end{aligned} $$

with $V_p = 1 + \| \cdot \|_p$. $\square$
**Lemma 19.** For $\mu, \nu$ two distributions on $\mathbb{R}^d$, we have

$$ \|\mu - \nu\|_V \le \sqrt{2} (\mu(V^2) + \nu(V^2))^{\frac{1}{2}} \|\mu - \nu\|_{TV}^{\frac{1}{2}}. $$

*Proof.* By Lemma 17 and the Cauchy-Schwarz inequality, we get

$$ \begin{aligned} \|\mu - \nu\|_V &= \inf_{(X,Y) \sim \beta} \mathbb{E}((V(X) + V(Y))1_{X \neq Y}) \\ &\le \inf_{(X,Y) \sim \beta} \mathbb{E}((V(X) + V(Y))^2)^{\frac{1}{2}} \mathbb{E}(1_{X \neq Y})^{\frac{1}{2}} \\ &\le \inf_{(X,Y) \sim \beta} \mathbb{E}(2(V^2(X) + V^2(Y)))^{\frac{1}{2}} \mathbb{E}(1_{X \neq Y})^{\frac{1}{2}} \\ &\le (2(\mu(V^2) + \nu(V^2)))^{\frac{1}{2}} \|\mu - \nu\|_{TV}^{\frac{1}{2}}, \end{aligned} $$

using that $\|\mu - \nu\|_{TV} = \mathbb{E}_{(X,Y) \sim \mu \times \nu}(1_{X \neq Y})$. $\square$

**Lemma 20.** For $T > 0$, let $b^1, b^2: C([0, T], \mathbb{R}^d) \times [0, T] \to \mathbb{R}^d$ be two drifts such that $\forall i \in \{1, 2\}$, $dx_t^i = b^i((x_u^i)_{u \in [0, T]}, t)dt + \sqrt{2}dw_t$ admits a strong solution with $x_0^i \sim \mu$ with the Markov semi-group $(P_t^i)_{t \in \mathbb{R}_+}$ (see definition in equation (39)) and $(w_t)_{t \ge 0}$ a Wiener process. Moreover, assume that $\mathbb{P}[\int_0^T \|b^i((x_u^i)_{u \in [0, T]}, t)\|^2 + \|b^i((w_u)_{u \in [0, T]}, t)\|^2 dt < +\infty] = 1$. Let $V: \mathbb{R}^d \to [1, +\infty)$ measurable, then we have

$$ \|\mu P_T^1 - \mu P_T^2\|_V \le (\mu P_T^1(V^2) + \mu P_T^2(V^2))^{\frac{1}{2}} \left( \int_0^T \|b^1((x_u^1)_{u \in [0,T]}, t) - b^2((x_u^1)_{u \in [0,T]}, t)\|^2 dt \right)^{\frac{1}{2}}. $$

Lemma 20 is a generalization of [56, Lemma 19] for any distribution $\mu$ on $\mathbb{R}^d$ instead of Dirac distributions $\delta_x$.

*Proof.* We follow the sketch of the proof of [56, Lemma 19]. For $i \in \{1, 2\}$, we denote by $\mu^i$ the distribution of $(X_t^i)_{t \in [0, T]}$ on the Wiener space $(C([0, T], \mathbb{R}), \mathcal{B}(C([0, T], \mathbb{R})))$, with $x_0^i \sim \mu$ and $dx_t^i = b^i((x_u^i)_{u \in [0, T]}, t)dt + \sqrt{2}dw_t$. We denote by $\mu_W$ the distribution of $(w_t)_{t \in [0, T]}$, the Wiener process.

By the generalized Pinsker's inequality [26, Lemma 24], we have

$$ \|\mu P_T^1 - \mu P_T^2\|_V \le \sqrt{2} (\mu P_T^1(V^2) + \mu P_T^2(V^2))^{\frac{1}{2}} \mathrm{KL}^{\frac{1}{2}}(\mu^1|\mu^2). \quad (45) $$

We assumed $\mathbb{P}[\int_0^T \|b^i((x_u^i)_{u \in [0, T]}, t)\|^2 + \|b^i((w_u)_{u \in [0, T]}, t)\|^2 dt < +\infty] = 1$ for any $i \in \{1, 2\}$. Hence we can apply Girsanov's Theorem [58, Theorem 7.7] and get $\mu_W$-almost surely that

$$ \begin{aligned} \frac{d\mu^1}{d\mu_W}((w_u)_{u \in [0,T]}, t) &= \exp \left[ \frac{1}{2} \int_0^T \langle b^1((w_u)_{u \in [0,T]}, t), dw_t \rangle - \frac{1}{4} \int_0^T \|b^1((w_u)_{u \in [0,T]}, t)\|^2 dt \right] \\ \frac{d\mu_W}{d\mu^2}((w_u)_{u \in [0,T]}, t) &= \exp \left[ -\frac{1}{2} \int_0^T \langle b^2((w_u)_{u \in [0,T]}, t), dw_t \rangle + \frac{1}{4} \int_0^T \|b^2((w_u)_{u \in [0,T]}, t)\|^2 dt \right], \end{aligned} $$

where $\langle \cdot, \cdot \rangle$ is the canonic scalar product on $\mathbb{R}^d$. Note that $\int_0^T \langle b^i((w_u)_{u \in [0, T]}, t), dw_t \rangle$, for $i \in \{1, 2\}$, are stochastic integrals.

Hence, we get

$$ \begin{aligned} \mathrm{KL}(\mu^1|\mu^2) &= \mathbb{E} \left( \log \frac{d\mu^1}{d\mu^2}((x_u^1)_{u \in [0,T]}, t) \right) \\ &= \frac{1}{4} \int_0^T \mathbb{E} \left[ \|b^1((x_u^1)_{u \in [0,T]}, t) - b^2((x_u^1)_{u \in [0,T]}, t)\|^2 \right]. \end{aligned} \quad (46) $$

By putting together equation (45) and equation (46), we get the desired result.

<figure><img src="image_12.png" /></figure>
### F.3 Sampling stability for Markov Chain in $V_p$-norm

In this part, we demonstrate the sampling stability of discrete Markov processes in the $V_p$-norm. Theorem 5 is a significant refinement of [71, Theorem 1] as it does not involve any discretization error and it generalizes previous known results in $W_1$ and TV-norm to the $V_p$-norms, which will imply the result in TV-norm and $W_p$ distance for any $p \in \mathbb{N}^*$.

We introduce two discrete-time Markov Chains $X_k^i$, $i \in \{1, 2\}$, defined for step-size $\gamma > 0$ by

$$X_{k+1}^i = X_k^i - \gamma b^i(X_k^i) + \sqrt{2\gamma} Z_{k+1}^i \quad (47)$$

with $Z_{k+1}^i \sim \mathcal{N}(0, I_d)$, two independent standard Gaussian distributions, and two drifts $b^i \in \mathcal{C}^0(\mathbb{R}^d, \mathbb{R}^d)$.

We will show a similar result than Theorem 1 in $V_p$-norm for $V_p = 1 + \|x\|^p$ with $p \in \mathbb{N}^*$ and then deduce Theorem 1 in TV-norm and $W_p$ (see Section F.4).

**Theorem 5.** If $b^1$ and $b^2$ satisfy Assumption 1, the two Markov Chains $X_k^1$ and $X_k^2$ are geometrically ergodic, with invariant laws $\pi_\gamma^1, \pi_\gamma^2$. Moreover, for $\gamma_0 = \frac{m}{L^2}$, with $L$, $m$ defined in Assumption 1, $p \in \mathbb{N}^*$ and $V_p(x) = 1 + \|x\|^p$, there exists $G_p \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, we have

$$\|\pi_\gamma^1 - \pi_\gamma^2\|_{V_p} \le G_p \left( \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right) \right)^{\frac{1}{2}}.$$

Note that the pseudo-distance between the drift $\|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)} = \left(\mathbb{E}_{Y \sim \pi_\gamma^1} \left(\|b^1(Y) - b^2(Y)\|^2\right)\right)^{\frac{1}{2}}$ is not symmetric in the drifts $b^1, b^2$ which allows to evaluate this pseudo-distance only by sampling $\pi_\gamma^1$, see [71] for empirical evaluations. It could be symmetries with $\min(\|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}, \|b^1 - b^2\|_{\ell_2(\pi_\gamma^2)})$ instead of $\|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}$.

### Proof. Scheme of the Proof

In a first part, we demonstrate that the two discrete-time Markov Chains $X_k^1, X_k^2$ are geometrically ergodic. Thus, we can define and study their invariant laws $\pi_\gamma^1$ and $\pi_\gamma^2$.

In a second part, we study the distance between these two invariant laws. The idea is to approximate them by the law of the processes at finite time $k\gamma$ starting at a point $x \in \mathbb{R}^d$, i.e. $\delta_x R_{\gamma,1}^k$ and $\delta_x R_{\gamma,2}^k$. Then the ergodicity allows us to quantify the distance between the two laws at finite time $k\gamma$ by the distance between the two laws at a bounded time $(m\gamma \approx 1)$ in the proof. The Girsanov theorem (Lemma 20) gives a quantification of the two laws at a bounded time by a distance on the drift. Finally, we look at the asymptotic $k \to +\infty$ of the obtained control to get the desired result on the invariant laws.

To simplify the notations, without loss of generality, we assume that $b^1$ and $b^2$ verify Assumption 1 with the same constants $L, R \ge 0$ and $m > 0$.

An introduction to Markov continuous-time processes and Markov discrete-time processes is given in Appendix F.1 and technical lemmas proofs can be founded in Appendix F.2.

### Proof that $X_k^1, X_k^2$ are geometrically ergodic

First, by [22, Corollary 2] and under Assumption 1 with $\gamma_0 = \frac{m}{L^2}$ (see [22, Proposition 12] for the expression of $\gamma_0$), there exist $\rho \in (0, 1)$ and $D \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, $\forall i \in \{1, 2\}$, we have

$$\mathbf{W}_1(\delta_x R_{\gamma,i}^k, \delta_y R_{\gamma,i}^k) \le D \rho^{k\gamma} \|x - y\|, \quad (48)$$

with $\delta_x$ the dirac distribution on $x \in \mathbb{R}^d$ and $R_{\gamma,i}$ the Markov kernel of the Markov Chain $i$ defined in equation (47) as introduced in Appendix F.1. Note that $\delta_x R_{\gamma,i}^k$ represents the law of the Markov Chain $X_k^i$ at step $k$ knowing that $X_0^i = x$. We want to generalize the previous inequality for any distributions $\mu, \nu$ on $\mathbb{R}^d$ instead of $\delta_x, \delta_y$. We denote $\beta^*$ the optimal transport plan from $\mu$ to $\nu$, by definition of the $\mathbf{W}_1$ distance, we have

$$\mathbf{W}_1(\mu, \nu) = \int_{x_1, x_2 \in \mathbb{R}^d \times \mathbb{R}^d} \|x_1 - x_2\| \beta^*(dx_1, dx_2). \quad (49)$$
By the Kantorovitch-Rubinstein Theorem [30, Theorem 4.1] [49], the dual formulation of the $W_1$ distance writes

$$W_1(\mu, \nu) = \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ \mathrm{Lip}(\phi) \le 1}} \int_{\mathbb{R}^d} \phi(d\mu - d\nu),$$

with $\mathrm{Lip}(\phi)$ the Lipschitz constant of $\phi$.

Then, by equation (48), equation (49) and the dual formulation of $W_1$ and the fact that $\delta_{x_1} R_{\gamma,i}^k(dx) = R_{\gamma,i}^k(x_1, dx)$ thanks to the Markov kernel definition (42), we get

$$
\begin{align*}
D\rho^{k\gamma} \mathbf{W}_1(\mu, \nu) &= \int_{\mathbb{R}^d \times \mathbb{R}^d} D\rho^{k\gamma} \|x_1 - x_2\| \beta^*(dx_1, dx_2) \\
&\geq \int_{\mathbb{R}^d \times \mathbb{R}^d} \mathbf{W}_1(\delta_{x_1} R_{\gamma,i}^k, \delta_{x_2} R_{\gamma,i}^k) \beta^*(dx_1, dx_2) \\
&= \int_{\mathbb{R}^d \times \mathbb{R}^d} \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ \mathrm{Lip}(\phi) \le 1}} \int_{\mathbb{R}^d} \phi(x) \left( (\delta_{x_1} R_{\gamma,i}^k)(dx) - (\delta_{x_2} R_{\gamma,i}^k)(dx) \right) \beta^*(dx_1, dx_2) \\
&\geq \sup_{f: \mathbb{R}^d \to \mathbb{R} \\ \mathrm{Lip}(\phi) \le 1} \int_{\mathbb{R}^d} \phi(x) \int_{\mathbb{R}^d \times \mathbb{R}^d} \left( (\delta_{x_1} R_{\gamma,i}^k)(dx) - (\delta_{x_2} R_{\gamma,i}^k)(dx) \right) \beta^*(dx_1, dx_2) \\
&= \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ \mathrm{Lip}(\phi) \le 1}} \int_{\mathbb{R}^d} \phi(x) \int_{\mathbb{R}^d \times \mathbb{R}^d} (R_{\gamma,i}^k(x_1, dx) - R_{\gamma,i}^k(x_2, dx)) \beta^*(dx_1, dx_2) \\
&= \sup_{\substack{\phi: \mathbb{R}^d \to \mathbb{R} \\ \mathrm{Lip}(\phi) \le 1}} \int_{\mathbb{R}^d} \phi(x) \left( (\mu R_{\gamma,i}^k)(dx) - (\nu R_{\gamma,i}^k)(dx) \right) \\
&= \mathbf{W}_1(\mu R_{\gamma,i}^k, \nu R_{\gamma,i}^k).
\end{align*}
$$

So we get for any distributions $\mu, \nu$ on $\mathbb{R}^d$,

$$\mathbf{W}_1(\mu R_{\gamma,i}^k, \nu R_{\gamma,i}^k) \le D\rho^{k\gamma} \mathbf{W}_1(\mu, \nu).$$

This proves that for $k$ large enough, the Markov kernel $R_{\delta,i}^k$ is contractive with respect to the 1-Wasserstein distance. Moreover as $(P_1(\mathbb{R}^d), W_1)$ is a complete space [89, Theorem 6.18], the Picard fixed point theorem shows that there exists a unique $\pi_{\gamma}^i$, such that for $\forall k \ge 0$, $R_{\gamma,i}^k \pi_{\gamma}^i = \pi_{\gamma}^i$. This probability law $\pi_{\gamma}^i$ is called the invariant law of the Markov Chain.

## Approximation of the invariant laws by finite time laws

We proved that $\pi_{\gamma}^1$ and $\pi_{\gamma}^2$ are well defined, we can now tackle the problem of quantifying the distance between these two distributions in $V_p$-norm. By the triangle inequality, for $k \in \mathbb{N}$ and

$$m = \left\lfloor \frac{1}{\gamma} \right\rfloor \qquad (50)$$

we have

$$\|\pi_{\gamma}^1 - \pi_{\gamma}^2\|_{V_p} \le \|\pi_{\gamma}^1 - \delta_x R_{\gamma,1}^{km}\|_{V_p} + \|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p} + \|\delta_x R_{\gamma,2}^{km} - \pi_{\gamma}^2\|_{V_p} \quad (51)$$

We need to control two type of term in the previous equation, the contractiveness in $V_p$-norm for the first and the third terms and the shift between the two law at time $k\gamma$ for the second term.

## Contractiveness in $V_p$-norm

We study the contractiveness of the Markov kernel in $V_p$-norm. By Lemma 19 and [22, Corollary 2], there exist $\rho_p \in (0, 1)$ and $E_p \ge 0$ such that $\forall x, y \in \mathbb{R}^d$ and $\forall i \in \{1, 2\}$, we have

$$
\begin{align}
\|\delta_x R_{\gamma,i}^k - \delta_y R_{\gamma,i}^k\|_{V_p} &\le \sqrt{2} (\delta_x R_{\gamma,i}^k(V_p^2) + \delta_y R_{\gamma,i}^k(V_p^2))^{\frac{1}{2}} \|\delta_x R_{\gamma,i}^k - \delta_y R_{\gamma,i}^k\|_{TV}^{\frac{1}{2}} \nonumber \\
&\le \sqrt{2} E_p (\delta_x R_{\gamma,i}^k(V_p^2) + \delta_y R_{\gamma,i}^k(V_p^2))^{\frac{1}{2}} \rho_p^{k\gamma} (1 + \|x - y\|)^{\frac{1}{2}}. \tag{52}
\end{align}
$$
By [56, Lemma 16] and [56, Lemma 17], under Assumption 1, for any $m \in \mathbb{N}$ and $V_{2m}(x) = 1 + \|x\|^{2m}$ there exists $M_m \ge 0$ such that $\forall k \in \mathbb{N}$ and $\forall i \in \{1, 2\}$, we have

$$R_{\gamma,i}^k V_{2m}(x) \le M_{2m} V_{2m}(x). \qquad (53)$$

Thus, given that $V_p^2 \le 2V_{2p}$ and $\delta_x(V_p^2) = V_p^2(x)$, we have $\forall x \in \mathbb{R}^d$

$$\delta_x R_{\gamma,i}^k V_p^2 \le 2M_{2p} V_{2p}(x). \qquad (54)$$

By injecting equation (54) into equation (52), using that $\sqrt{a+b} \le \sqrt{a} + \sqrt{b}$ for $a, b \ge 0$ and that $\sqrt{V_{2p}} \le V_p$, we get that for all $x \in \mathbb{R}^d$

$$\begin{aligned} \|\delta_x R_{\gamma,i}^k - \delta_y R_{\gamma,i}^k\|_{V_p} &\le 2E_p \sqrt{M_{2p}} \sqrt{V_{2p}(x) + V_{2p}(y)} \rho_p^{k\gamma} (1 + \|x-y\|)^{\frac{1}{2}} \\ &\le 2E_p \sqrt{M_{2p}} (V_p(x) + V_p(y)) \rho_p^{k\gamma} (1 + \|x\| + \|y\|)^{\frac{1}{2}}. \end{aligned}$$

Using the inequality $\sqrt{1 + \|x\| + \|y\|} \le 1 + \|x\| + \|y\| \le 1 + \|x\| + 1 + \|y\| \le 2^{1-\frac{1}{p}}(1 + \|x\|^p)^{\frac{1}{p}} + 2^{1-\frac{1}{p}}(1 + \|y\|^p)^{\frac{1}{p}} \le 2(V_p(x) + V_p(y))$, we get

$$\|\delta_x R_{\gamma,i}^k - \delta_y R_{\gamma,i}^k\|_{V_p} \le D_p \rho_p^{k\gamma} (V_p^2(x) + V_p^2(y)), \qquad (55)$$

with $D_p := 8E_p \sqrt{M_{2p}}$. This inequality generalizes [56, Proposition 5] that was demonstrated only for $V_1$.

By integrating equation (55) for $x \sim \mu$ and $y \sim \nu$ and the triangle inequality of the V-norm, we get for any probability measures $\mu, \nu$

$$\|\mu R_{\gamma,i}^k - \nu R_{\gamma,i}^k\|_{V_p} \le D_p \rho_p^{k\gamma} (\mu(V_p^2) + \nu(V_p^2)). \qquad (56)$$

This proves in particular that with $\nu = \pi_\gamma^i$, which satisfies $\pi_\gamma^i R_{\gamma,i}^k = \pi_\gamma^i$, the Markov Chain $X_k^i$ is geometrically ergodic, i.e. the empirical measure of the Markov Chain at time $k \ge 0$ of one realization converges in $V_p$-norm geometrically to the invariant law of the process.

$$\|\pi_\gamma^1 - \delta_x R_{\gamma,1}^{km}\|_{V_p} + \|\delta_x R_{\gamma,2}^{km} - \pi_\gamma^2\|_{V_p} \le D_p(x) \rho_p^{km\gamma} \qquad (57)$$

with $D_p(x) := D_p(2V_p^2(x) + \pi_\gamma^1(V_p^2) + \pi_\gamma^2(V_p^2))$, thanks to $\delta_x(V_p^2) = V_p^2(x)$. This is a control for the first and the third terms of equation (51)

### Control of the shift of the two laws at time $k\gamma$

The distribution $\delta_x R_{\gamma,i}^k$ is the distribution at time $t = k\gamma$ of the discrete-time Markov Chain $X_k^i$ starting at $x$. The term $\|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p}$ quantifies the deformation at time $t = k\gamma$ between the laws of the two processes starting at the same point run with two different drifts.

By the triangle inequality, we get

$$\begin{aligned} \|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p} &= \left\| \sum_{j=0}^{k-1} \delta_x R_{\gamma,1}^{(k-j)m} R_{\gamma,2}^{jm} - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^{(j+1)m} \right\|_{V_p} \\ &\le \sum_{j=0}^{k-1} \|\delta_x R_{\gamma,1}^{(k-j)m} R_{\gamma,2}^{jm} - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^{(j+1)m}\|_{V_p} \\ &= \sum_{j=0}^{k-1} \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m R_{\gamma,2}^{jm} - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m R_{\gamma,2}^{jm}\|_{V_p}. \end{aligned} \qquad (58)$$

For $\phi : \mathbb{R}^d \to \mathbb{R}$ such that $\|\phi\| \le V_p$ and due to $\pi_\gamma^2(\phi) \in \mathbb{R}$, we get

$$\begin{aligned} & (\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m R_{\gamma,2}^{jm} - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m R_{\gamma,2}^{jm})(\phi) \\ &= (\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m) R_{\gamma,2}^{jm}(\phi) \\ &= (\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m) (R_{\gamma,2}^{jm}(\phi) - \pi_\gamma^2(\phi)). \end{aligned} \qquad (59)$$
By equation (56) applied with $\mu = \delta_x$, $\nu = \pi_\gamma^2$, $i = 2$ and $k = jm$, thanks to $V_p(x) \ge 1$, we have

$$ \|\delta_x R_{\gamma,2}^{jm} - \pi_\gamma^2\|_{V_p} \le B \rho_p^{jm\gamma} V_p^2(x), $$

with $B := D_p(1 + \pi_\gamma^2(V_p^2))$.

By the definition of the $V_p$-distance (43), we have

$$ |\delta_x R_{\gamma,2}^{jm}(\phi) - \pi_\gamma^2(\phi)| \le B \rho_p^{jm\gamma} V_p^2(x). $$

Using that $\delta_x R_{\gamma,2}^{jm}(\phi) = R_{\gamma,2}^{jm}(\phi)(x)$ and that $V_p^2 \le 2V_{2p}$, we deduce that

$$ \left| \frac{R_{\gamma,2}^{jm}(f)(x) - \pi_\gamma^2(f)}{2B\rho_p^{jm\gamma}} \right| \le V_{2p}(x). \quad (60) $$

Going back to equation (59), thanks to equation (60), we get

$$ \begin{aligned} & (\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m R_{\gamma,2}^{jm} - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m R_{\gamma,2}^{jm})(\phi) \\ & \le |\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m| |R_{\gamma,2}^{jm}(\phi) - \pi_\gamma^2(\phi)| \\ & \le 2B\rho_p^{jm\gamma} |\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m| \left| \frac{R_{\gamma,2}^{jm}(\phi) - \pi_\gamma^2(\phi)}{2B\rho^{jm\gamma}} \right| \\ & \le 2B\rho_p^{jm\gamma} \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m\|_{V_{2p}}. \end{aligned} $$

Taking the supremum on $f$ in the last inequality, we get

$$ \begin{aligned} & \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m R_{\gamma,2}^{jm} - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m R_{\gamma,2}^{jm}\|_{V_p} \\ & \le 2B\rho_p^{jm\gamma} \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m\|_{V_{2p}}. \end{aligned} \quad (61) $$

By injecting the previous inequality into equation (58), we get

$$ \| \delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km} \|_{V_p} \le 2B \sum_{j=0}^{k-1} \rho_p^{jm\gamma} \| \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m \|_{V_{2p}}. \quad (62) $$

### Application of the Girsanov Theorem to get control at bounded time $m\gamma$

The previous inequality shows that we can bound $\|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p}$ with quantities of the form $\|\mu R_{\gamma,1}^m - \mu R_{\gamma,2}^m\|_{V_p}$. This last quantity has a temporal duration of $m\gamma \in (1 - \gamma, 1]$ due to $m = \lfloor \frac{1}{\gamma} \rfloor$, which allows us to apply Girsanov's theory to bound it by the drift difference. To do so, we need to introduce a continuous version of the Markov Chain:

$$ dx_t^i = -b_\gamma^i((x_u^i)_{u \in [t-\gamma, t]}, t)dt + \sqrt{2}dw_t^i, \quad (63) $$

with $(w_t^i)_{t \ge 0}$ a Wiener process and $b_\gamma^i((x_u)_{u \in [t-\gamma, t]}, t) = b^i(x_{\lfloor \frac{t}{\gamma} \rfloor \gamma})$, with $b^i$ defined in equation (47).

Here the drift $b_\gamma^i$ is non-homogeneous, i.e. it depends on the time $t$, and it is path-dependent. This process is a continuous counterpart of the Markov Chain (47). In fact, by integrating equation (63) for $t \in [k\gamma, (k+1)\gamma)$, we get

$$ \begin{aligned} \int_{k\gamma}^{(k+1)\gamma} dx_t^i &= - \int_{k\gamma}^{(k+1)\gamma} b^i(x_{k\gamma}^i) dt + \sqrt{2} \int_{k\gamma}^{(k+1)\gamma} dw_t^i \\ x_{(k+1)\gamma}^i - x_{k\gamma}^i &= -\gamma b^i(x_{k\gamma}^i) + \sqrt{2\gamma} Z_{k+1}, \end{aligned} $$

with $Z_{k+1} \sim \mathcal{N}(0, I_d)$ because $\int_{k\gamma}^{(k+1)\gamma} dw_t^i \sim \mathcal{N}(0, \gamma I_d)$. Therefore, $x_{k\gamma}^i$ initialized with the probability law $\mu$ has the same law than $X_k^i$ defined in equation (47) initialized with the law $\mu$.
Under Assumption 1, the drift is Lipschitz with respect to the infinite norm, i.e. for $\phi, \psi \in C([0, \gamma], \mathbb{R}^d)$

$$ \begin{aligned} \|b_{\gamma}^{i}((\phi(u))_{u \in [0, \gamma]}, t) - b_{\gamma}^{i}((\psi(u))_{u \in [0, \gamma]}, t)\| &= \left\| b^{i}(\phi(\lfloor \frac{t}{\gamma} \rfloor \gamma)) - b^{i}(\psi(\lfloor \frac{t}{\gamma} \rfloor \gamma)) \right\| \\ &\le L \left\| \phi(\lfloor \frac{t}{\gamma} \rfloor \gamma) - \psi(\lfloor \frac{t}{\gamma} \rfloor \gamma) \right\| \le L \|\phi - \psi\|_{\infty}. \end{aligned} $$

By [60, Theorem 5.2.2], equation (63) has a unique strong solution $x_t^i$ and

$$ \mathbb{E}(\|b_{\gamma}^{i}((x_{u}^{i})_{u \in [t-\gamma, t]}, t)\|) < +\infty. $$

So assumptions of Lemma 20 are verified, and we get that for $T > 0$

$$ \|\mu P_T^1 - \mu P_T^2\|_{V_p} \le (\mu P_T^1(V_p^2) + \mu P_T^2(V_p^2))^{\frac{1}{2}} \left( \int_0^T \|b^1(x_t^1, t) - b^2(x_t^1, t)\|^2 dt \right)^{\frac{1}{2}}, \quad (64) $$

with $V_p = 1 + \| \cdot \|^p, P_t^i$ for $i \in \{1, 2\}$ the Markov semi-group at time $t$ (define in equation (39)) of the process following equation (63) and $(x_t^i)_{t \ge 0}$ the stochastic process following equation (63) with the initial law $x_0^i \sim \mu$. For $T = m\gamma$, we have that $P_{m\gamma}^i = R_{\gamma,i}^m$ so

$$ \begin{aligned} & \|\mu R_{\gamma,1}^m - \mu R_{\gamma,2}^m\|_{V_{2p}} \\ & \le (\mu R_{\gamma,1}^m(V_{2p}^2) + \mu R_{\gamma,2}^m(V_{2p}^2))^{\frac{1}{2}} \left( \int_0^{m\gamma} \|b^1(x_t^1, t) - b^2(x_t^1, t)\|^2 dt \right)^{\frac{1}{2}}. \end{aligned} \quad (65) $$

Now, we can apply inequality (65) with $\mu = \delta_x R_{\gamma,1}^{(k-j-1)m}$ to get

$$ \begin{aligned} & \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m\|_{V_{2p}} \\ & \le \left( \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m(V_{2p}^2) + \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m(V_{2p}^2) \right)^{\frac{1}{2}} \\ & \times \left( \int_0^{m\gamma} \mathbb{E} \left( \|b^1(x_{t,j}^1, t) - b^2(x_{t,j}^1, t)\|^2 \right) dt \right)^{\frac{1}{2}}, \end{aligned} \quad (66) $$

with $x_{t,j}^1$ the process starting with the law $\delta_x R_{\gamma,1}^{(k-j-1)m}$ at time $t = 0$ and verifying equation (63).

Using equation (53) we have $R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m(V_{2p}) \le M_{2p}V_{2p}$. Recalling that $\delta_x(V_p^2) = V_p^2(x)$ and $V_p^2 \le 2V_{2p}$, we obtain a uniform bound

$$ \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m(V_{2p}^2) \le 2M_{4p}V_{4p}(x). $$

We also have that $R_{\gamma,2}^m(V_p^2) \le 2M_{4p}V_{4p}$ so, using again (53), we get

$$ \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m(V_{2p}^2) \le 2M_{4p} \delta_x R_{\gamma,1}^{(k-j-1)m} V_{2p} \le 2M_{4p}^2 V_{4p}(x). $$

Then we obtain

$$ \left( \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m(V_{2p}^2) + \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m(V_{2p}^2) \right)^{\frac{1}{2}} \le \left( 2M_{4p}V_{4p}(x) + 2M_{4p}^2V_{4p}(x) \right)^{\frac{1}{2}}. $$

By injecting the previous moment bounds in equation (66), we get

$$ \begin{aligned} & \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m\|_{V_{2p}} \\ & \le A_1(x) \left( \int_0^{m\gamma} \mathbb{E} \left( \|b^1(x_{t,j}^1, t) - b^2(x_{t,j}^1, t)\|^2 \right) dt \right)^{\frac{1}{2}} \\ & = A_1(x) \left( \gamma \sum_{l=0}^{m-1} \mathbb{E} \left( \|b^1(x_{l,\gamma,j}^1) - b^2(x_{l,\gamma,j}^1)\|^2 \right) \right)^{\frac{1}{2}}, \end{aligned} $$
with $A_1(x) = (2M_{4p}V_{4p}(x) + 2M_{4p}^2V_{4p}(x))^{\frac{1}{2}}$ a constant independent of $j$ and $k$. Moreover, $x_{l\gamma,j}^1$ has the same law that $X_{(k-j-1)m+l}^1$.

Hence we get that

$$ \begin{aligned} & \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m\|_{V_{2p}} \\ & \le A_1(x) \left( \gamma \sum_{l=0}^{m-1} \mathbb{E} \left( \|b^1(X_{(k-j-1)m+l}^1) - b^2(X_{(k-j-1)m+l}^1)\|^2 \right) \right)^{\frac{1}{2}}, \end{aligned} \quad (67) $$

**From bounded time $m\gamma$ to any finite time $k\gamma$ with ergodicity**

We will massively use in this paragraph the ergodicity of the discrete-time Markov Chain to deduce from the control at time $m\gamma$ in equation (67) a control for any time $k$ on $\|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p}$.

We define $Y_k$ the Markov Chain following equation (47) and starting with the invariant law $\pi_\gamma^1$ at time $k=0$. Then the law of the Markov Chain $Y_k$ is $\pi_\gamma^1$ at each time $k \ge 0$. Next, by Assumption 1, $b^1$ and $b^2$ are L-Lipschitz and we have

$$ \begin{aligned} & \sum_{l=0}^{m-1} \mathbb{E} \left( \|b^1(X_{(k-j-1)m+l}^1) - b^2(X_{(k-j-1)m+l}^1)\|^2 \right) \\ & \le 2 \sum_{l=0}^{m-1} \mathbb{E} \left( \|b^1(X_{(k-j-1)m+l}^1) - b^2(X_{(k-j-1)m+l}^1) - b^1(Y_k) + b^2(Y_k)\|^2 \right) \\ & + 2 \sum_{l=0}^{m-1} \mathbb{E}_{Y_k \sim \pi_\gamma^1} \left( \|b^1(Y_k) - b^2(Y_k)\|^2 \right) \\ & \le 8L^2 \sum_{l=0}^{m-1} \mathbb{E} \left( \|X_{(k-j-1)m+l}^1 - Y_k\|^2 \right) dt + 2m \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right). \end{aligned} \quad (68) $$

To bound the first term of the above relation, we use [22, Corollary 2], which states, under Assumption 1, that there exist $A_2(x) \ge 0$ and $\tilde{\rho} \in (0,1)$ such that $\forall k \ge 0$

$$ \|\delta_x R_{\gamma,i}^k - \pi_\gamma^i\|_{TV} \le A_2(x) \tilde{\rho}^{k\gamma}. $$

Then, the previous inequality, the moment bound in equation (53) and the inequality $(x+y)^4 \le 8(x^4 + y^4)$ give that

$$ \begin{aligned} & \mathbb{E} \left( \|X_{(k-j-1)m+l}^1 - Y_k\|^2 \right) = \mathbb{E} \left( \|X_{(k-j-1)m+l}^1 - Y_k\|^2 \mathbb{1}_{X_{(k-j-1)m+l}^1 \neq Y_k} \right) \\ & \le \mathbb{E} \left( \|X_{(k-j-1)m+l}^1 - Y_k\|^4 \right)^{\frac{1}{2}} \mathbb{E} \left( \mathbb{1}_{X_{(k-j-1)m+l}^1 \neq Y_k} \right)^{\frac{1}{2}} \\ & = \mathbb{E} \left( \|X_{(k-j-1)m+l}^1 - Y_k\|^4 \right)^{\frac{1}{2}} \|\delta_x R_{\gamma,1}^{(k-j-1)m+l} - \pi_\gamma^1\|_{TV}^{\frac{1}{2}} \\ & \le 2\sqrt{2} \mathbb{E} \left( \|X_{(k-j-1)m+l}^1\|^4 + \|Y_k\|^4 \right)^{\frac{1}{2}} C_2^{\frac{1}{2}}(x) \tilde{\rho}^{\frac{(k-j-1)m+l}{2}} \\ & \le 4\sqrt{2M_4V_4(x)A_2(x)} \tilde{\rho}^{\frac{(k-j-1)m+l}{2}}, \end{aligned} $$

with $V_4(x) = 1 + \|x\|^4$ and $M_4$ a moment bound defined in equation (53).

By combining the previous inequality with equation (68), we get

$$ \begin{aligned} & \sum_{l=0}^{m-1} \mathbb{E} \left( \|b^1(X_{(k-j-1)m+l}^1) - b^2(X_{(k-j-1)m+l}^1)\|^2 \right) \\ & \le 32L^2 \sqrt{2C_4W_4(x)A_2(x)} \sum_{l=0}^{m-1} \tilde{\rho}^{\frac{(k-j-1)m+l}{2}} + 2m \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right) \\ & \le A_3(x) \tilde{\rho}^{\frac{(k-j-1)m}{2}} + 2m \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right), \end{aligned} \quad (69) $$
with $A_3(x) = 32L^2\sqrt{2C_4W_4(x)A_2(x)}\frac{1}{1-\tilde{\rho}^{\frac{m}{2}}}$ a constant independent of $k$ and $j$. By injecting the previous equation in equation (67) and given that $\gamma m \le 1$ thanks to $m$ definition in equation (50), we get

$$ \begin{aligned} & \|\delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,1}^m - \delta_x R_{\gamma,1}^{(k-j-1)m} R_{\gamma,2}^m\|_{V_{2p}} \\ & \le A_1(x) \left(A_3(x)\gamma \tilde{\rho}^{\frac{(k-j-1)m\gamma}{2}} + 2\gamma m \mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2)\right)^{\frac{1}{2}} \\ & \le A_1(x) A_3^{\frac{1}{2}}(x) \gamma^{\frac{1}{2}} \tilde{\rho}^{\frac{(k-j-1)m\gamma}{4}} + \sqrt{2} A_1(x) \left(\mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2)\right)^{\frac{1}{2}}. \end{aligned} $$

We now combine the previous inequality with equation (62) to obtain

$$ \begin{aligned} & \|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p} \\ & \le 2BA_1(x) \sum_{j=0}^{k-1} \rho_p^{jm\gamma} \left( A_3^{\frac{1}{2}}(x) \gamma^{\frac{1}{2}} \tilde{\rho}^{\frac{(k-j-1)m\gamma}{4}} + \sqrt{2} \left( \mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2) \right)^{\frac{1}{2}} \right) \\ & \le 2BA_1(x) A_3^{\frac{1}{2}}(x) \gamma^{\frac{1}{2}} \frac{1}{\left|1 - \frac{\rho_p^{m\gamma}}{\tilde{\rho}^{\frac{m\gamma}{4}}}\right|} \tilde{\rho}^{\frac{(k-j-1)m\gamma}{4}} \\ & \quad + \sqrt{22}BA_1(x) \frac{1}{1 - \rho_p^{m\gamma}} \left( \mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2) \right)^{\frac{1}{2}}, \end{aligned} $$

where, at the cost of slightly increasing $\tilde{\rho}$, we assume that $\frac{\rho_p^{m\gamma}}{\tilde{\rho}^{\frac{m\gamma}{4}}} \ne 1$. Denoting $A_4(x) = 2BA_1(x)A_3^{\frac{1}{2}}(x)\gamma^{\frac{1}{2}}\frac{\tilde{\rho}^{\frac{m\gamma}{4}}}{\left|1 - \frac{\rho_p^{m\gamma}}{\tilde{\rho}^{\frac{m\gamma}{4}}}\right|}$ and $A_5(x) = \sqrt{22}BA_1(x)\frac{1}{1-\rho_p^{m\gamma}}$, we get

$$ \|\delta_x R_{\gamma,1}^{km} - \delta_x R_{\gamma,2}^{km}\|_{V_p} \le A_4(x) \tilde{\rho}^{\frac{km\gamma}{4}} + A_5(x) \left(\mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2)\right)^{\frac{1}{2}}. \quad (70) $$

## Putting all together and asymptotic computations

We combine equation (51), equation (57) and equation (70) to obtain

$$ \begin{aligned} \|\pi_{\gamma}^1 - \pi_{\gamma}^2\|_{V_p} \le A_4(x) \tilde{\rho}^{\frac{km\gamma}{4}} + A_5(x) \left( \mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2) \right)^{\frac{1}{2}} \\ + D_p(x) \rho_p^{km\gamma}. \end{aligned} $$

By taking $k \to +\infty$, we get

$$ \|\pi_{\gamma}^1 - \pi_{\gamma}^2\|_{V_p} \le A_5(x) \left( \mathbb{E}_{Y\sim\pi_{\gamma}^1} (\|b^1(Y) - b^2(Y)\|^2) \right)^{\frac{1}{2}}, $$

with $A_5(x) = \sqrt{22}BA_1(x)\frac{1}{1-\rho_p^{m\gamma}}$, with $m\gamma \ge 1 - \gamma \ge 1 - \gamma_0$, since we defined $m = \lfloor \frac{1}{\gamma} \rfloor$ in relation (50) and we supposed that $\gamma \in (0, \gamma_0)$ in the theorem's assumptions. Note that $x \in \mathbb{R}^d$ is free. We choose $x = 0$, then $B = D_p(1 + \pi_{\gamma}^2(V_p^2))$ which is defined after equation (59) and $A_1(0) = (2M_{4p} + 2M_{4p}^2)^{\frac{1}{2}}$ which is defined before equation (67). By taking $k \to +\infty$ in equation (54), we get that $\forall x \in \mathbb{R}^d$, $\pi_{\gamma}^2(V_p^2) \le 2M_{2p}V_{2p}(x)$. In particular for $x = 0$, we get $\pi_{\gamma}^2(V_p^2) \le 2M_{2p}$, thus $B \le D_p(1 + 2M_{2p})$.

This proves Theorem 1 in p-Wasserstein distance with the constant

$$ G_p = \frac{\sqrt{22}}{1 - \rho_p^{1-\gamma_0}} D_p(1 + 2M_{2p}) (2M_{4p} + 2M_{4p}^2)^{\frac{1}{2}}, $$

where the constants $D_p$ and $\rho_p$ are defined by the geometric convergence in equation (55), $M_p$ by the moment bound in (53), $\gamma_0 = \frac{m}{L^2}$ with $m, L$ defined in Assumption 1 and $V_p = 1 + \|\cdot\|^p$. $\square$
## F.4 Proof of Theorem 1

In this part, we prove Theorem 1. First, we recall its statement.

**Theorem 6.** Let $b^1$ and $b^2$ satisfy Assumption 1. $X_k^1$ and $X_k^2$ are two geometrically ergodic Markov Chains with invariant laws $\pi_\gamma^1, \pi_\gamma^2$. Then for $\gamma_0 = \frac{m}{L^2}$ and $p \in \mathbb{N}^*$ there exist $C_p, C \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, we have

$$ 
\begin{aligned}
\mathbf{W}_p(\pi_\gamma^1, \pi_\gamma^2) &\le C_p \|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}^{\frac{1}{p}} \\
\|\pi_\gamma^1 - \pi_\gamma^2\|_{TV} &\le \tilde{C} \|b^1 - b^2\|_{\ell_2(\pi_\gamma^1)}.
\end{aligned}
$$

Theorem 1 is demonstrated for the $\mathbf{W}_p$ and TV metrics. In [28], the authors analyse the convergence of the Langevin dynamics using the Kullback-Leibler divergence. However, because our approach depends heavily on the triangle inequality property of the metric, it does not readily extend to the Kullback-Leibler divergence, which does not satisfy this inequality.

*Proof.* By Theorem 5, under Assumption 1, we have that $X_k^1$ and $X_k^2$ are geometrically ergodic, with invariant laws $\pi_\gamma^1, \pi_\gamma^2$. Moreover, with $\gamma_0 = \frac{m}{L^2}$, for $p \in \mathbb{N}^*$, there exists $G_p \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, we have

$$ \|\pi_\gamma^1 - \pi_\gamma^2\|_{V_p} \le G_p \left( \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right) \right)^{\frac{1}{2}}, $$

with $V_p = 1 + \| \cdot \|^p$.

However, by definition of the V-norm in equation (43), we have

$$ \|\pi_\gamma^1 - \pi_\gamma^2\|_{TV} \le \|\pi_\gamma^1 - \pi_\gamma^2\|_{V_1} \le G_1 \left( \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right) \right)^{\frac{1}{2}}, $$

which proves the first part of the Theorem relative to the TV-distance.

Next, by Lemma 18, we have for $p \in \mathbb{N}^*$

$$ 
\begin{aligned}
\mathbf{W}_p (\pi_\gamma^1, \pi_\gamma^2) &\le 2^{1-\frac{1}{p}} \|\pi_\gamma^1 - \pi_\gamma^2\|_{V_p}^{\frac{1}{p}} \\
&\le 2^{1-\frac{1}{p}} G_p^{\frac{1}{p}} \left( \mathbb{E}_{Y \sim \pi_\gamma^1} \left( \|b^1(Y) - b^2(Y)\|^2 \right) \right)^{\frac{1}{2p}},
\end{aligned}
$$

which proves the second part for the $\mathbf{W}_p$ distance. $\square$

## F.5 On the constant $C_p, C$ of Theorem 1

Note that $C = C_1$, so we only need to study $C_p$ for $p \in \mathbb{N}$. The constant $C_p$ of Theorem 1 is defined by

$$ C_p = 2^{1-\frac{1}{p}} \left( \frac{2\sqrt{2}}{1-\rho_p^{1-\gamma_0}} D_p(1+2M_{2p}) (2M_{4p}+2M_{4p}^2)^{\frac{1}{2p}} \right)^{\frac{1}{p}}, \quad (71) $$

where the constants $D_p$ and $\rho_p$ are defined by the geometric convergence in equation (55), $M_p$ by the moment bound in (53), $\gamma_0 = \frac{m}{L^2}$ with $m, L$ defined in Assumption 1 and $V_p = 1 + \| \cdot \|^p$.

$C_p$ is independent of the step-size $\gamma$ and of the drifts $b^1, b^2$.

The authors of [21] proved that the constant $\rho_p$ is independent of the dimension $d$ and that $D_p$ is polynomial in the dimension. Therefore, $C_p$ has a polynomial dependence in the dimension $d$ with a dependence on the moment of the target law.

The dependency on $p$ is unclear as it depends highly on the evolution of the moment bounds $M_{2p}$ with $p$. Therefore, the dependence of $C_p$ in $p \in \mathbb{N}$ is unknown in general.
## F.6 Quantification of the discretization Error - Proof of Theorem 2

Based on the proof of Theorem 1 we can recover the discretization error for discrete Unadjusted Langevin Algorithm (ULA). For that purpose, we introduce the ULA Markov Chain

$$X_{k+1} = X_k - \gamma \nabla V(X_k) + \sqrt{2\gamma} Z_{k+1}, \quad (72)$$

with $Z_{k+1}$ a normalized Gaussian random variable, as well as the continuous Langevin Markov Process

$$dx_t^1 = -\gamma \nabla V(x_t^1) + \sqrt{2\gamma} dw_t^1, \quad (73)$$

with $dw_t^1$ a Wiener process. Under Assumption 1, Theorem 1 shows that $X_k$ has an invariant law $\pi_\gamma$ and is geometrically ergodic. Moreover, if $-\nabla V$ verifies Assumption 1, [22, Corollary 22] proves that $x_t^1$ has $\pi \propto e^{-V}$ as an invariant law and is geometrically ergodic.

**Theorem 7.** If the function $\nabla V$ verifies Assumption 1 and $\pi \propto e^{-V}$. Then, for $\gamma_0 = \frac{m}{L^2}$, for $p \in \mathbb{N}^*$ there exist $C_p, D_p, \tilde{C} \ge 0$ such that $\forall \gamma \in (0, \gamma_0]$, the Markov Chain is geometrically ergodic with an invariant law $\pi_\gamma$ and the discretization error is bounded by

$$\begin{aligned} \mathbf{W}_p (\pi_\gamma, \pi) &\le C_p \gamma^{\frac{1}{2p}} \\ \|\pi_\gamma - \pi\|_{TV} &\le \tilde{C} \gamma^{\frac{1}{2}} \\ \|\pi_\gamma - \pi\|_{V_p} &\le D_p \gamma^{\frac{1}{2}}, \end{aligned}$$

with $\pi = e^{-V}$ and $V_p = 1 + \|\cdot\|^p$.

*Proof.* This proof is based on the same reasoning that the proof of Theorem 5 with continuous Langevin processes instead of discrete-time Markov Chain. The key point is to write the discrete-time ULA Markov Chain as a continuous Langevin process in equation (74). Then, we verify that the step of the reasoning of the proof of Theorem 5 are still true in this continuous context.

We define the ULA continuized continuous process by

$$dx_t^2 = -b_\gamma((x_t^2)_{t \in \mathbb{R}^+}, t)dt + \sqrt{2\gamma} dw_t^2, \quad (74)$$

with $b_\gamma((x_t)_{t \in \mathbb{R}^+}, t) = \sum_{k=0}^{+\infty} \mathbb{1}_{[k\gamma, (k+1)\gamma)}(t) \nabla V(x_{k\gamma}) = \nabla V(x_{\lfloor \frac{t}{\gamma} \rfloor \gamma})$. Following the same reasoning as the one given after equation (63), the process $x_t^2$ is a continuous counterpart of the Markov Chain (72). More precisely, $x_{k\gamma}^2$ initialized with the probability law $\delta_{X_0}$ has the same law than $X_k$, the ULA Markov Chain defined in equation (72).

We denote by $P_t$ the semi-group (see definition in equation (39)) of the stochastic equation (73) and $R_\gamma$ the Markov kernel (see definition in equation (42)) of the Markov Chain defined in equation (72).

We follow a similar sketch of the proof of Theorem 5. We first demonstrate a relation similar to equation (56), by replacing [22, Corollary 2] with [22, Corollary 2] for the process $x_t^1$. We thus obtain that for $x \in \mathbb{R}^d$ and $p \in \mathbb{N}^*$, there exist $A_p \ge 0$ and $\rho_p \in (0, 1)$ such that, for $\mu$ a probability distribution on $\mathbb{R}^d$, we have

$$\begin{aligned} \|\pi_\gamma - \mu R_\gamma^{km}\|_{V_p} &\le A_p \mu (V_p^2) \rho_p^{km\gamma} \\ \|\mu P_{km\gamma} - \pi\|_{V_p} &\le A_p \mu (V_p^2) \rho_p^{km\gamma}. \end{aligned} \quad (75)$$

For $x \in \mathbb{R}^d$ and $k \in \mathbb{N}$, $m = \lfloor \frac{1}{\gamma} \rfloor$, by the triangle inequality and equation (75) applied for $\mu = \delta_x$, we have

$$\begin{aligned} \|\pi_\gamma - \pi\|_{V_p} &\le \|\pi_\gamma - \delta_x R_\gamma^{km}\|_{V_p} + \|\delta_x R_\gamma^{km} - \delta_x P_{km\gamma}\|_{V_p} + \|\delta_x P_{km\gamma} - \pi\|_{V_p} \\ &\le 2A_p V_p^2(x) \rho_p^{km\gamma} + \|\delta_x R_\gamma^{km} - \delta_x P_{km\gamma}\|_{V_p} \\ &\le 2A_p V_p^2(x) \rho_p^{km\gamma} + \sum_{j=0}^{k-1} \|\delta_x R_\gamma^{(k-j-1)m} R_\gamma^m P_{jm\gamma} - \delta_x R_\gamma^{(k-j-1)m} P_{m\gamma} P_{jm\gamma}\|_{V_p}. \end{aligned}$$
Similarly to the reasoning that led us to equation (61), we get

$$ \begin{aligned} & \|\delta_x R_{\gamma}^{(k-j-1)m} R_{\gamma}^m P_{jm\gamma} - \delta_x R_{\gamma}^{(k-j-1)m} P_{m\gamma} P_{jm\gamma}\|_{V_p} \\ & \le 2A_p \rho_p^{jm\gamma} \|\delta_x R_{\gamma}^{(k-j-1)m} R_{\gamma}^m - \delta_x R_{\gamma}^{(k-j-1)m} P_{m\gamma}\|_{V_{2p}}. \end{aligned} $$

By combining the two previous equations, we obtain

$$ \begin{aligned} & \|\pi_{\gamma} - \pi\|_{V_p} \\ & \le 2A_p V_p^2(x) \rho_p^{km\gamma} + 2A_p \sum_{j=0}^{k-1} \rho_p^{jm\gamma} \|\delta_x R_{\gamma}^{(k-j-1)m} R_{\gamma}^m - \delta_x R_{\gamma}^{(k-j-1)m} P_{m\gamma}\|_{V_{2p}}. \end{aligned} \quad (76) $$

Then, by Lemma 20 and the fact that the drift of $x_t^1$ is $\nabla V$ and the drift of $x_t^2$ is $b_\gamma$, we get

$$ \begin{aligned} & \|\delta_x R_{\gamma}^{(k-j-1)m} R_{\gamma}^m - \delta_x R_{\gamma}^{(k-j-1)m} P_{m\gamma}\|_{V_{2p}} \\ & \le \left( \delta_x R_{\gamma}^{(k-j-1)m} (V_{2p}^2) + \delta_x R_{\gamma}^{(k-j-1)m} (V_{2p}^2) \right)^{\frac{1}{2}} \\ & \times \left( \int_0^{m\gamma} \mathbb{E} \left( \|\nabla V(y_{t,j}) - b_\gamma((y_{t,j})_{t \in \mathbb{R}_+}, t)\|^2 \right) dt \right)^{\frac{1}{2}}, \end{aligned} $$

with $b_\gamma$ defined in equation (74) and $y_{t,j}$, the continuous process starting with the law $\delta_x R_{\gamma}^{(k-j-1)m}$ and following equation (73).

Using the moment bound given by [56, Lemma 17], we get that there exists $F_p(x) \ge 0$ such that $\forall k \in \mathbb{N}^*, j \in \{0, \dots, k-1\}$,

$$ \left( \delta_x R_{\gamma}^{(k-j-1)m} (V_{2p}^2) + \delta_x R_{\gamma}^{(k-j-1)m} (V_{2p}^2) \right)^{\frac{1}{2}} \le F_p(x). $$

Combining previous equations thus gives

$$ \|\delta_x R_{\gamma}^{(k-j-1)m} R_{\gamma}^m - \delta_x R_{\gamma}^{(k-j-1)m} P_{m\gamma}\|_{V_{2p}} \quad (77) $$

$$ \le F_p(x) \left( \int_0^{m\gamma} \mathbb{E} \left( \|\nabla V(y_{t,j}) - b_\gamma((y_{t,j})_{t \in \mathbb{R}_+}, t)\|^2 \right) dt \right)^{\frac{1}{2}}. \quad (78) $$

Since $\nabla V$ and $b_\gamma$ are L-Lipschitz from Assumption 1, we get, with $m = \lfloor \frac{1}{\gamma} \rfloor$, that

$$ \begin{aligned} & \int_0^{m\gamma} \mathbb{E} \left( \|\nabla V(y_{t,j}) - b_\gamma((y_{t,j})_{t \in \mathbb{R}_+}, t)\|^2 \right) dt \\ & \le L \left( \sum_{k=0}^{m-1} \int_{k\gamma}^{(k+1)\gamma} \mathbb{E} \left( \|y_{t,j} - y_{k\gamma,j}\|^2 \right) dt \right) \end{aligned} \quad (79) $$

Because, $y_{t,j}$ is following equation (73), we have $y_{t,j} - y_{k\gamma,j} = \int_{k\gamma}^t -\nabla V(y_{u,j}) du + \sqrt{2} \int_{k\gamma}^t dw_u^1$. Combined with the Itô isometry, the fact that $\|\nabla V(x)\| \le \|\nabla V(0)\| + L\|x\|$, the Cauchy-Schwarz inequality, the moment bound in equation (53) and the fact that $t - k\gamma \le \gamma$, we get

$$ \begin{aligned} \mathbb{E} (\|y_{t,j} - y_{k\gamma,j}\|^2) &= \mathbb{E} \left( \left\| \int_{k\gamma}^{t} - \nabla V(y_{u,j}) du + \sqrt{2} \int_{k\gamma}^{t} dw_u^1 \right\|^2 \right) \quad (80) \\ &\le 2\mathbb{E} \left( \left\| \int_{k\gamma}^{t} - \nabla V(y_{u,j}) du \right\|^2 \right) + 4\mathbb{E} \left( \left\| \int_{k\gamma}^{t} dw_u^1 \right\|^2 \right) \\ &\le 2\mathbb{E} \left( \left( \gamma \|\nabla V(0)\| + L \int_{k\gamma}^{t} \|y_{u,j}\| du \right)^2 \right) + 4\gamma d \\ &\le 4\gamma^2 \|\nabla V(0)\|^2 + 4L^2 \mathbb{E} \left( \left( \int_{k\gamma}^{t} \|y_{u,j}\| du \right)^2 \right) + 4\gamma d \\ &\le 4\gamma^2 \|\nabla V(0)\|^2 + 4L^2 \gamma^2 M_2 V_2(x) + 4\gamma d \\ &\le (4\gamma_0 \|\nabla V(0)\|^2 + 4L^2 \gamma_0 M_2 V_2(x) + 4d) \gamma. \end{aligned} \quad (81) $$
Combining equations (81), (79) and (78), with $B_6 = 4\gamma_0\|\nabla V(0)\|^2 + 4L^2\gamma_0 M_2 V_2(x) + 4d$ and the fact that $m\gamma \le 1$, we obtain

$$ \|\delta_x R_{\gamma}^{(k-j-1)m} R_{\gamma}^m - \delta_x R_{\gamma}^{(k-j-1)m} P_{m\gamma}\|_{V_{2p}} \le F_p(x) L^{\frac{1}{2}}(B_6 m \gamma)^{\frac{1}{2}} \gamma^{\frac{1}{2}} \le F_p(x) L^{\frac{1}{2}}(B_6)^{\frac{1}{2}} \gamma^{\frac{1}{2}}, $$

since $m = \lfloor \frac{1}{\gamma} \rfloor$ so that $m\gamma \le 1$. By injecting the previous inequality into the equation (76), we get

$$ \begin{aligned} \|\pi_{\gamma} - \pi\|_{V_p} &\le 2A_p V_p^2(x) \rho_p^{km\gamma} + 2A_p \sum_{j=0}^{k-1} \rho_p^{jm\gamma} F_p(x) (LB_6)^{\frac{1}{2}} \gamma^{\frac{1}{2}} \\ &\le 2A_p V_p^2(x) \rho_p^{km\gamma} + 2A_p \frac{1}{1 - \rho_p^{m\gamma}} (LB_6)^{\frac{1}{2}} \gamma^{\frac{1}{2}}. \end{aligned} $$

Taking $k \to +\infty$ in the previous inequality, we get the desired result in $V_p$-norm. We get the result in $W_p$ and $TV$ norms following the same reasoning as in the proof of Theorem 1 in section F.4. $\square$

# G Proofs for PSGLA

In this part, we provide proofs for the analysis of the PSGLA algorithm defined in equation (3). First, we demonstrate regularity results on the drift $b^\gamma$ defined in equation (11), then based on our IULA Markov Chain (definition in equation (2)) analysis of section 2, we derive the convergence theory for the PSGLA algorithm.

## G.1 The drift of PSGLA is verifying Assumption 1

**Lemma 21.** Under Assumptions 2-3, for all $\gamma \in (0, \min(\frac{1}{2L_g}, \frac{1}{\rho}, \gamma_1)]$, $b^\gamma$ verifies

(i) $b^\gamma$ is $L$-smooth on $\mathbb{R}^d$, with $L = 2L_f + 2L_g$.

(ii) $\forall x, y \in \mathbb{R}^d$ such that $\|x - y\| \ge 4R_0$, we have

$$ \langle b(x) - b(y), x - y \rangle \ge \frac{\mu}{4} \|x - y\|^2, $$

with the constants $L_f, L_g, \rho, \gamma_1, R_0$ and $\mu$ defined in Assumptions 2-3.

Then, for $\gamma \in (0, \min(\frac{1}{2L_g}, \frac{1}{\rho}, \gamma_1)]$, $b^\gamma$ verifies Assumption 1 with $L = 2L_f + 2L_g, R = 4R_0$ and $m = \frac{\mu}{4}$.

Note that the constant that appears in Lemma 21 are independent of $\gamma$.

*Proof.* We will demonstrate each point separately for the drift $b^\gamma(x) = \nabla f(x - \gamma \nabla g^\gamma(x)) + \nabla g^\gamma(x)$ defined in equation (11). We need that $\gamma \le \frac{1}{\rho}$ for the Moreau envelope $g^\gamma$ of $g$ to be well defined.

(i) Under Assumption 3(i), Lemma 15 of Appendix D.7 gives that for $\gamma \le \frac{1}{2L_g}$, $\nabla g^\gamma$ is $2L_g$-Lipschitz. Combined with Assumption 2(i), we get, for $x, y \in \mathbb{R}^d$

$$ \begin{aligned} & \|b^\gamma(x) - b^\gamma(y)\| \\ &\le \|\nabla f(x - \gamma \nabla g^\gamma(x)) - \nabla f(y - \gamma \nabla g^\gamma(y))\| + \|\nabla g^\gamma(x) - \nabla g^\gamma(y)\| \\ &\le L_f \|x - y - \gamma \nabla g^\gamma(x) + \gamma \nabla g^\gamma(y)\| + 2L_g \|x - y\| \\ &\le L_f (\|x - y\| + \gamma 2L_g \|x - y\|) + 2L_g \|x - y\| \\ &\le (2L_f + 2L_g) \|x - y\|, \end{aligned} $$

which proves point (i) of Lemma 21.
(ii) We want to have a lower bound of the quantity $\langle b^{\gamma}(x) - b^{\gamma}(y), x - y \rangle$, for $x, y \in \mathbb{R}^d$, when $\|x - y\|$ is sufficiently large.

For $0 < \gamma \le \frac{1}{2L_g}$ and $x, y \in \mathbb{R}^d$, we have

$$ \langle b^{\gamma}(x) - b^{\gamma}(y), x - y \rangle \ge -2L_f \|x - y\|^2 + \langle \nabla g^{\gamma}(x) - \nabla g^{\gamma}(y), x - y \rangle, $$

where we used that $\nabla f(x - \gamma \nabla g^{\gamma}(x))$ if $2L_f$-Lipschitz thanks to the computation of point (i).

We suppose that $\|x - y\| \ge 4R_0$, then there is at most one point between $x$ and $y$ that belongs to the ball $B(0, R_0)$.

We distinguish three cases.

• $[x, y] \cap B(0, R_0) = \emptyset$. From Assumption 3(ii), if $\gamma \le \gamma_1$, $\nabla^2 g^{\gamma} \ge \mu I_d$ on $\mathbb{R}^d \setminus B(0, R_0)$. Then we immediately have

$$ \langle b^{\gamma}(x) - b^{\gamma}(y), x - y \rangle \ge (\mu - 2L_f) \|x - y\|^2. $$

• If $x \in B(0, R_0)$, then $\|y\| = \|y - x + x\| \ge \|y - x\| - \|x\| \ge 3R_0$. We take $z \in [x, y]$, such that $\|z\| = R_0$, then $|z, y] \subset \mathbb{R}^d \setminus B(0, R_0)$. We define $t_0 \in [0, 1]$, such that $z = x + t_0(y - x)$. It holds that $\|z - x\| \le 2R_0$ and $\|y - z\| = \|y - x + x - z\| \ge \|y - x\| - \|x - z\| \ge 2R_0$ and $t_0 = \frac{\|z - x\|}{\|y - x\|} \le \frac{1}{2}$.

Under Assumption 3(i), by Lemma 15 of Appendix D.7, $g^{\gamma}$ is $2L_g$-smooth, we get

$$ \begin{aligned} \langle \nabla g^{\gamma}(y) - \nabla g^{\gamma}(x), y - x \rangle &= \int_{0}^{1} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle dt \\ &= \int_{0}^{t_0} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle dt \\ &+ \int_{t_0}^{1} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle dt \\ &\ge -2L_g t_0 \|y - x\|^2 + (1 - t_0) \mu \|y - x\|^2 \\ &\ge -L_g \|y - x\|^2 + \frac{\mu}{2} \|y - x\|^2, \end{aligned} $$

and we deduce that

$$ \langle b^{\gamma}(x) - b^{\gamma}(y), x - y \rangle \ge \left( \frac{\mu}{2} - 2L_f - L_g \right) \|y - x\|^2. $$

• $x, y \notin B(0, R_0)$, we look at the affine segment $t \in [0, 1] \mapsto x + t(y - x)$. If this segment does not intersect with $B(0, R_0)$, then we are in the first case. Thus we suppose that the segment crosses the ball $B(0, R_0)$. We look at the norm of the vectors $u(\alpha)$ defined, for $\alpha \in \mathbb{R}$, along the affine line that passes through points $x$ and $y$, as

$$ u(\alpha) = \|x + \alpha(y - x)\|^2 = \|x\|^2 + 2\alpha \langle x, y - x \rangle + \alpha^2 \|y - x\|^2. $$

Notice that $u$ is quadratic in $\alpha$. We assumed that $u(0), u(1) > R_0^2$ and there exists $t \in [0, 1]$ such that $u(t) < R_0^2$, then there necessarily exist $0 < t_0 < t_1 < 1$, such that $u(t_0) = u(t_1) = R_0^2$. As a consequence, $\forall t \in (t_0, t_1)$, $u(t) < R_0^2$ and $\forall t \in [0, 1] \setminus [t_0, t_1]$, $u(t) > R_0^2$. We denote

$$ z_0 = x + t_0(y - x) $$

$$ z_1 = x + t_1(y - x). $$

By definition $z_0, z_1 \in B(0, R_0)$, therefore $\|z_0 - z_1\| \le 2R_0$. Then

$$ t_1 - t_0 = \frac{\|z_1 - x\|}{\|y - x\|} - \frac{\|z_0 - x\|}{\|y - x\|} = \frac{\|z_1 - z_0\|}{\|y - x\|} \le \frac{2R_0}{\|y - x\|} \le \frac{1}{2}, \quad (82) $$
because $\|y - x\| \ge 4R_0$. We can now estimate the quantity of interest

$$
\begin{align*}
\langle \nabla g^{\gamma}(y) - \nabla g^{\gamma}(x), y - x \rangle &= \int_{0}^{1} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle \mathrm{d}t \\
&= \int_{0}^{t_0} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle \mathrm{d}t \\
&\quad + \int_{t_0}^{t_1} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle \mathrm{d}t \\
&\quad + \int_{t_1}^{1} \langle \nabla^2 g^{\gamma}(x + t(y - x))(y - x), y - x \rangle \mathrm{d}t \\
&\geq (1 + t_0 - t_1)\mu\|y - x\|^2 - 2(t_1 - t_0)L_g\|y - x\|^2 \\
&\geq (1 + t_0 - t_1)\mu\|y - x\|^2 - L_g\|y - x\|^2 \\
&\geq \frac{\mu}{2}\|y - x\|^2 - L_g\|y - x\|^2,
\end{align*}
$$

where we used equation (82) in the last relation. So, we get

$$
\langle b^{\gamma}(y) - b^{\gamma}(x), y - x \rangle \geq \left(\frac{\mu}{2} - 2L_f - L_g\right) \|y - x\|^2.
$$

To conclude, in the three cases, for $\|y - x\| \ge 4R_0$ and by Assumption 3(ii) $8L_f + 4L_g \le \mu$, we get

$$
\begin{align*}
\langle b(y) - b(x), y - x \rangle &\geq \left(\frac{\mu}{2} - 2L_f - L_g\right) \|y - x\|^2. \\
&\geq \frac{\mu}{4} \|y - x\|^2,
\end{align*}
$$

which proves the point (ii) of Lemma 21.

□

G.2 Proof of Theorem 3

First, we recall the statement of Theorem 3.

**Theorem 8.** Under Assumptions 2-3, there exist $r \in (0, 1)$, $C_1, C_2 \in \mathbb{R}_+$ such that $\forall \gamma \in (0, \bar{\gamma}]$, with $\bar{\gamma} = \min(\frac{1}{2L_g}, \frac{\mu}{32(L_f+L_g)^2}, \frac{1}{2\rho}, \gamma_1)$, where $L_g, L_f, \rho, \gamma_1$ are defined in Assumptions 2-3, and $\forall k \in \mathbb{N}$, we have

$$
\mathbf{W}_p(p_{Y_k}, \mu_\gamma) \le C_1 r^{k\gamma} + C_2 \gamma^{\frac{1}{2p}}, \quad (83)
$$

with $p_{Y_k}$ the distribution of $Y_k$ and $\mu_\gamma \propto e^{-f-g^\gamma}$.

Moreover there exist $C_3, C_4 \in \mathbb{R}_+$ such that $\forall \gamma \in (0, \bar{\gamma}]$ and $\forall k \in \mathbb{N}$, we have

$$
\mathbf{W}_p(p_{X_k}, \nu_\gamma) \le C_3 r^{k\gamma} + C_4 \gamma^{\frac{1}{2p}}, \quad (84)
$$

with $p_{X_k}$ the distribution of $X_k$ and $\nu_\gamma \propto \text{Prox}_{\gamma g} \# e^{-f-g^\gamma}$.

*Proof.* Thanks to Lemma 21, under Assumptions 2-3, the drift of the Markov Chain $Y_k$ $b^\gamma$ verifies Assumption 1. Therefore, using Theorem 5, we deduce that this Markov Chain $Y_k$ is geometrically ergodic, so there exist $C_1 \ge 0$, $r \in (0, 1)$ and an invariant law $p_\infty$ such that

$$
\mathbf{W}_1(p_{Y_k}, p_\infty) \leq C_1 r^{k\gamma}, \quad (85)
$$

with $p_{Y_k}$ the distribution of $Y_k$.

Next we introduce the Markov Chain $\tilde{Y}_k$ such that $\tilde{Y}_0 = Y_0$ and

$$
\tilde{Y}_{k+1} = \tilde{Y}_k + \gamma b_0 (\tilde{Y}_k) + \sqrt{2\gamma} \tilde{Z}_{k+1},
$$
with $b_0 = -\nabla g^\gamma - \nabla f$ and $\tilde{Z}_{k+1} \sim \mathcal{N}(0, I_d)$. Following the same proof than Lemma 21, we deduce that $b_0$ verifies Assumption 1 for $\gamma \le \min(\frac{1}{\rho}, \frac{1}{2L_g}, \gamma_1)$ with the same parameters than $b^\gamma$, $L = 2L_f + 2L_g$, $R = 4R_0$ and $m = \frac{\mu}{4}$. The parameters $L_f, L_g, R_0, \mu, \gamma_1, \rho$ are defined in Assumptions 2-3. Theorem 1 gives that for $\gamma_0 = \frac{\mu}{32(L_f+L_g)^2}$ and $\gamma \in (0, \gamma_0]$, $\tilde{Y}_k$ geometrically converges to its invariant law $p_\gamma$ and there exists $C_2 \ge 0$ such that

$$ \mathbf{W}_p(p_\infty, p_\gamma) \le C_2 \left( \mathbb{E}_{X \sim p_\infty} \left( \|b(X) - b_0(X)\|^2 \right) \right)^{\frac{1}{2p}}. $$

We define $\tilde{\gamma} = \min\left(\frac{1}{\rho}, \frac{1}{2L_g}, \gamma_1, \frac{\mu}{32(L_f+L_g)^2}\right)$.

Then, using the definition of $b$ in equation (11) and $b_0$, Assumption 1, Assumption 2(i), Assumption 3(i), the fact that $\nabla g^\gamma(x) = \nabla g(\text{Prox}_{\gamma g}(x))$ (see Lemma 12 of Appendix D.6) and Lemma 11 of Appendix D.5, we get

$$ \begin{aligned} \mathbf{W}_p(p_\infty, p_\gamma) &\le C_2 \left( \mathbb{E}_{X \sim p_\infty} \left( \|\nabla f(X - \gamma \nabla g^\gamma(X)) - \nabla f(X)\|^2 \right) \right)^{\frac{1}{2p}} \\ &\le C_2 L_f^{\frac{1}{p}} \gamma^{\frac{1}{p}} \left( \mathbb{E}_{X \sim p_\infty} \left( \|\nabla g^\gamma(X)\|^2 \right) \right)^{\frac{1}{2p}} \\ &\le C_2 L_f^{\frac{1}{p}} \gamma^{\frac{1}{p}} \left( \mathbb{E}_{X \sim p_\infty} \left( \|\nabla g(\text{Prox}_{\gamma g}(x))\|^2 \right) \right)^{\frac{1}{2p}} \\ &\le C_2 L_f^{\frac{1}{p}} \gamma^{\frac{1}{p}} \left( 2L_g^2 \mathbb{E}_{X \sim p_\infty} \left( \|\text{Prox}_{\gamma g}(x) - \text{Prox}_{\gamma g}(0)\|^2 \right) + 2\|\text{Prox}_{\gamma g}(x)\|^2 \right)^{\frac{1}{2p}} \\ &\le C_2 L_f^{\frac{1}{p}} \gamma^{\frac{1}{p}} \left( \frac{2L_g^2}{(1-\gamma\rho)^2} \mathbb{E}_{X \sim p_\infty} \left( \|X\|^2 \right) + 2\|\text{Prox}_{\gamma g}(0)\|^2 \right)^{\frac{1}{2p}}. \end{aligned} $$

By equation (53), we know that for $\gamma \le \tilde{\gamma}$, $\mathbb{E}_{X \sim p_\infty} (\|X\|^2)$ is bounded by a constant independent of $\gamma$. Then, for $\gamma \le \tilde{\gamma} \le \frac{1}{2\rho}$, we have $\frac{1}{1-\gamma\rho} \le 2$. Moreover $\text{Prox}_{\gamma g}(0) \to 0$ when $\gamma \to 0$. So there exists a constant $C_3 \ge 0$, independent of $\gamma$, such that, for $\gamma \in (0, \tilde{\gamma}]$

$$ \mathbf{W}_1(p_\infty, p_\gamma) \le C_3 \gamma^{\frac{1}{p}}. \quad (86) $$

Note that $b_0 = -\nabla(g^\gamma + f)$, so by Theorem 2, there exists $C_4 \ge 0$ such that $\forall \gamma \in (0, \tilde{\gamma}]$ we have

$$ \mathbf{W}_p(p_\gamma, \mu_\gamma) \le C_4 \gamma^{\frac{1}{2p}}, \quad (87) $$

with $\mu_\gamma \propto e^{-f-g^\gamma}$.

Combining equations (85), (86) and (87) and the triangle inequality, we get

$$ \begin{aligned} \mathbf{W}_p(p_{Y_k}, \mu_\gamma) &\le C_1 r^{k\gamma} + C_3 \gamma^{\frac{1}{p}} + C_4 \gamma^{\frac{1}{2p}} \\ &\le C_1 r^{k\gamma} + C_5 \gamma^{\frac{1}{2p}}, \end{aligned} $$

with $C_5 = C_3 \tilde{\gamma}^{\frac{1}{2p}} + C_4$. This proves the first part of Theorem 3.

Due to $X_k = \text{Prox}_{\gamma g}(Y_k)$, the distribution $p_{X_k}$ of $X_k$ is $p_{X_k} = \text{Prox}_{\gamma g} \# p_{Y_k}$, with $p_{Y_k}$ the distribution of $Y_k$ and "#" the push forward operator defined for a measurable set $A \subset \mathbb{R}^d$ by $p_{X_k}(A) = p_{Y_k}(\text{Prox}_{\gamma g}^{-1}(A))$. With $\nu_\gamma = \text{Prox}_{\gamma g} \# \mu_\gamma$, we get

$$ \mathbf{W}_p(p_{X_k}, \nu_\gamma) = \min_{\beta \in \Lambda} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x_1 - x_2\|^p d\beta(x_1, x_2), $$

with $\Lambda$ the set of all possible transport plans $\beta$ between $p_{X_k}$ and $\nu_\gamma$. By the push-forward definition and Lemma 11 of Appendix D.5, we obtain

$$ \begin{aligned} \mathbf{W}_p(p_{X_k}, \nu_\gamma) &= \min_{\beta \in \Lambda} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|\text{Prox}_{\gamma g}(x_1) - \text{Prox}_{\gamma g}(x_2)\|^p d\beta(x_1, x_2) \\ &\le \frac{1}{(1-\gamma\rho)^{2p}} \min_{\beta \in \Lambda} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x_1 - x_2\| d\beta(x_1, x_2) \\ &\le \frac{1}{(1-\tilde{\gamma}\rho)^{2p}} \mathbf{W}_p(p_{Y_k}, \mu_\gamma). \end{aligned} $$

The last equation demonstrates the second part of Theorem 3. $\square$
### G.3 On the exact law approximation with PSGLA

In this section, we provide a proof for Proposition 1. First, we recall this proposition

**Proposition 4.** With $\pi \propto e^{-f-g}$, $\mu_\gamma \propto e^{-f-g^\gamma}$ and $\nu_\gamma = \text{Prox}_{\gamma g}\#\mu_\gamma$, we have, for $V: \mathbb{R}^d \to [1, +\infty)$ and $p \ge 1$

$$ \lim_{\gamma \to 0} \|\mu_\gamma - \pi\|_V = 0 \qquad (88) $$

$$ \lim_{\gamma \to 0} \mathbf{W}_p(\mu_\gamma, \pi) = 0 \qquad (89) $$

$$ \lim_{\gamma \to 0} \mathbf{W}_p(\nu_\gamma, \pi) = 0 \qquad (90) $$

Moreover, if $g$ is L-Lipschitz, there exists $E_p \in \mathbb{R}_+$ such that $\forall \gamma \in [0, \frac{2}{L^2}]$

$$ \|\mu_\gamma - \pi\|_V \le 2\pi(V)L\gamma \qquad (91) $$

$$ \mathbf{W}_p(\mu_\gamma, \pi) \le E_p(L^2\gamma)^{\frac{1}{p}} \qquad (92) $$

$$ \mathbf{W}_p(\nu_\gamma, \pi) \le E_p(L^2\gamma)^{\frac{1}{p}} + L\gamma. \qquad (93) $$

*Proof.* This proof is a generalization of the strategy proposed in the proof of [27, Proposition 3.1]. By the definition of the $V$-norm in equation (43), we have

$$ \|\mu_\gamma - \pi\|_V = \sup_{|\phi| \le V} \int_{\mathbb{R}^d} \phi(x) (\mu_\gamma(x) - \pi(x)) dx, $$

with $\mu_\gamma(x) = \frac{e^{-f(x)-g^\gamma(x)}}{\int e^{-f-g^\gamma}}$ and $\pi(x) = \frac{e^{-f(x)-g(x)}}{\int e^{-f-g}}$. Then

$$ \begin{aligned} \|\mu_\gamma - \pi\|_V &\le \sup_{|\phi| \le V} \int_{\mathbb{R}^d} |\phi(x)| |\mu_\gamma(x) - \pi(x)| dx \\ &\le \int_{\mathbb{R}^d} V(x)\pi(x) \left|1 - e^{g(x)-g^\gamma(x)}\right| \frac{e^{-f-g}}{\int e^{-f-g^\gamma}} dx \\ &\le \int_{\mathbb{R}^d} V(x)\pi(x) \left|1 - e^{g(x)-g^\gamma(x)}\right| dx + \int_{\mathbb{R}^d} V(x)\pi(x) e^{g(x)-g^\gamma(x)} \left|1 - \frac{e^{-f-g}}{\int e^{-f-g^\gamma}}\right| dx. \end{aligned} $$

However, we know that $\forall x \in \mathbb{R}^d$, $g^\gamma(x) \le g(x)$, so we get

$$ \|\mu_\gamma - \pi\|_V \le \int_{\mathbb{R}^d} V(x)\pi(x) (e^{g(x)-g^\gamma(x)} - 1) dx + \int_{\mathbb{R}^d} V(x)\pi(x) e^{g(x)-g^\gamma(x)} \left(1 - \frac{e^{-f-g}}{\int e^{-f-g^\gamma}}\right) dx. \quad (94) $$

For $x \in \mathbb{R}^d$, $\lim_{\gamma \to 0} g^\gamma(x) = g(x)$ and the function $\gamma \to g^\gamma(x)$ is decreasing. So by the monotone convergence theorem, we have

$$ \lim_{\gamma \to 0} \|\mu_\gamma - \pi\|_V = 0. $$

Lemma 18 proves that $\lim_{\gamma \to 0} \mathbf{W}_p(\mu_\gamma, \pi) = 0$. This concludes the first part of Proposition 1.

Next, by the triangle inequality, Lemma 11 of Appendix D.5 and taking the particular coupling $(I_d, \text{Prox}_{\gamma g})\#\pi$ between $\pi$ and $\text{Prox}_{\gamma g}\#\pi$, we get

$$ \begin{aligned} \mathbf{W}_p(\nu_\gamma, \pi) &\le \mathbf{W}_p(\text{Prox}_{\gamma g}\#\mu_\gamma, \text{Prox}_{\gamma g}\#\pi) + \mathbf{W}_p(\text{Prox}_{\gamma g}\#\pi, \pi) \\ &\le \frac{1}{1 - \gamma \rho} \mathbf{W}_p(\mu_\gamma, \pi) + \left( \inf_{\beta} \int \|x_1 - x_2\|^p d\beta(x_1, x_2) \right)^{\frac{1}{p}} \\ &\le \frac{1}{1 - \gamma \rho} \mathbf{W}_p(\mu_\gamma, \pi) + \left( \int \|x - \text{Prox}_{\gamma g}(x)\|^p d\pi(x) \right)^{\frac{1}{p}}. \end{aligned} $$

By the definition of the Moreau envelope, we have $g(\text{Prox}_{\gamma g}(x)) + \frac{1}{2\gamma} \|x - \text{Prox}_{\gamma g}(x)\|^2 \le g(x)$, so we get

$$ \mathbf{W}_p(\nu_\gamma, \pi) \le \frac{1}{1 - \gamma \rho} \mathbf{W}_p(\mu_\gamma, \pi) + \sqrt{2\gamma} \left( \int |g(x) - g(\text{Prox}_{\gamma g}(x))|^\frac{p}{2} d\pi(x) \right)^{\frac{1}{p}}. $$
However, by Lemma 2 of Appendix D.2, we have that $|g(x) - g(\text{Prox}_{\gamma g}(x))|$ has a monotone convergence to 0. By the monotone convergence theorem, we obtain

$$ \lim_{\gamma \to 0} \mathbf{W}_p(\nu_\gamma, \pi) = 0. $$

If $g$ is L-Lipschitz, then, for $x \in \mathbb{R}^d$, we have

$$ \begin{aligned} g(x) - g^{\gamma}(x) &= g(x) - \inf_{y \in \mathbb{R}^d} \left( \frac{1}{2\gamma} \|x - y\|^2 + g(y) \right) = \sup_{y \in \mathbb{R}^d} g(x) - g(y) - \frac{1}{2\gamma} \|x - y\|^2 \\ &\le \sup_{y \in \mathbb{R}^d} L \|x - y\| - \frac{1}{2\gamma} \|x - y\|^2 = \gamma \frac{L^2}{2}. \end{aligned} $$

Since $\gamma(x) \ge g^{\gamma}(x)$ for all $x \in \mathbb{R}^d$, we get $\|g - g^{\gamma}\|_{\infty} \le \gamma \frac{L^2}{2}$. By injecting this inequality into equation (94), we get

$$ \begin{aligned} & \|\mu_{\gamma} - \pi\|_{V} \\ & \le \int_{\mathbb{R}^d} V(x)\pi(x) \left(e^{\gamma \frac{L^2}{2}} - 1\right) dx + \int_{\mathbb{R}^d} V(x)\pi(x) e^{\gamma \frac{L^2}{2}} \left(1 - \frac{\int e^{-f-g}}{\int e^{-f-g^{\gamma}}}\right) dx. \end{aligned} $$

Thanks to $\int \frac{e^{-f-g}}{e^{-f-g^{\gamma}}} \ge e^{-\gamma \frac{L^2}{2}}$, we get

$$ \begin{aligned} \|\mu_{\gamma} - \pi\|_{V} &\le \pi(V) \left(e^{\gamma \frac{L^2}{2}} - 1\right) + \pi(V) e^{\gamma \frac{L^2}{2}} \left(1 - e^{-\gamma \frac{L^2}{2}}\right) \\ &\le 2\pi(V) \left(e^{\gamma \frac{L^2}{2}} - 1\right). \end{aligned} $$

For $u \in [0, 1]$, $e^u - 1 \le 2u$, so for $\gamma \le \frac{2}{L^2}$, we have

$$ \|\mu_{\gamma} - \pi\|_{V} \le 2\pi(V)\gamma L^2. $$

Using Lemma 18 of Appendix F.2, we obtain

$$ \mathbf{W}_p(\mu_\gamma, \pi) \le 2 \left( \pi(V_p) \gamma L^2 \right)^{\frac{1}{p}}. $$

which concludes for the inequality in $\mathbf{W}_p$ distance between $\mu_\gamma$ and $\pi$, with $E_p = 2 (\pi(V_p))^{\frac{1}{p}}$.

Using Lemma 12, the fact that $\|\nabla g\|_{\infty} \le L$ if $g$ is L-Lipschitz, the previous bound on $\mathbf{W}_p(\mu_\gamma, \pi)$ and the particular coupling $(I_d, \text{Prox}_{\gamma g})\#\mu_\gamma$ between $\mu_\gamma$ and $\nu_\gamma$, we get

$$ \begin{aligned} \mathbf{W}_p(\nu_\gamma, \pi) &\le \mathbf{W}_p(\nu_\gamma, \mu_\gamma) + \mathbf{W}_p(\mu_\gamma, \pi) \le \left( \int \|\gamma \nabla g(\text{Prox}_{\gamma g}(x))\|^p d\mu_\gamma(x) \right)^{\frac{1}{p}} + \mathbf{W}_p(\mu_\gamma, \pi) \\ &\le \gamma L + 2 \left(\pi(V_p)\right)^{\frac{1}{p}} \left(L^2\gamma\right)^{\frac{1}{p}}. \end{aligned} $$

This demonstrates Proposition 1 with $E_p = 2(\pi(V_p))^{\frac{1}{p}}$. $\square$

## G.4 Proof of Theorem 4

The Inexact PSGLA is defined with $S_\gamma$ an inexact approximation of $\text{Prox}_{\gamma g}$ by

$$ \hat{X}_{k+1} = S_{\gamma} \left( \hat{X}_{k} - \gamma \nabla f(\hat{X}_{k}) + \sqrt{2\gamma} \hat{Z}_{k+1} \right), $$

with $\hat{Z}_{k+1} \sim \mathcal{N}(0, I_d)$ i.i.d. We can write this algorithm as a two points algorithm

$$ \begin{aligned} \hat{Y}_{k+1} &= \hat{X}_{k} - \gamma \nabla f(\hat{X}_{k}) + \sqrt{2\gamma} \hat{Z}_{k+1} \\ \hat{X}_{k+1} &= S_{\gamma} (\hat{Y}_{k+1}). \end{aligned} $$