By introducing $G_{\gamma} = \gamma^{-1}(I_d - S_{\gamma})$, we have that

$$
\begin{align*}
\hat{Y}_{k+1} &= S_{\gamma}(\hat{Y}_k) - \gamma \nabla f(S_{\gamma}(\hat{Y}_k)) + \sqrt{2\gamma} \hat{Z}_{k+1} \\
&= \hat{Y}_k - \gamma G_{\gamma}(\hat{Y}_k) - \gamma \nabla f(\hat{Y}_k - \gamma G_{\gamma}(\hat{Y}_k)) + \sqrt{2\gamma} \hat{Z}_{k+1} \\
&= \hat{Y}_k - \gamma \hat{b}^{\gamma}(\hat{Y}_k) + \sqrt{2\gamma} \hat{Z}_{k+1},
\end{align*}
$$

with the drift $\hat{b}^{\gamma}(y) = \nabla f(y - G_{\gamma}(y)) + G_{\gamma}(y)$.

We recall Theorem 4.

**Theorem 9.** 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}], \hat{Y}_k$ of equation (14) has an invariant law $\hat{\mu}_{\gamma}$ and $\forall k \in \mathbb{N}$, we have

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

Theorem 4 is a generalization of [32, Theorem 3.14] as it holds for non-convex potentials. Moreover, in our study, we do not assume a priori that the inexactness of the proximal operator is bounded.

*Proof.* We assume that $\hat{b}^{\gamma}$ verifies Assumption 1, so $\hat{Y}_k$ is geometrically ergodic with an invariant law $\hat{\mu}_{\gamma}$ and there exist $B_0 \ge 0$ and $r \in (0, 1)$, such that

$$
\mathbf{W}_p(p_{\hat{Y}_k}, \hat{\mu}_\gamma) \le B_0 r^{k\gamma}, \quad (95)
$$

with $p_{\hat{Y}_k}$ the distribution of $\hat{Y}_k$

We define the sequence $Y_{k+1} = Y_k - \gamma b^\gamma(Y_k) + \sqrt{2\gamma} Z_{k+1}$, with $b^\gamma(y) = \nabla f(y - \gamma \nabla g^\gamma(y)) + \nabla g^\gamma(y)$ defined in equation (11). Thanks to Theorem 3, under Assumption 2-3, $b^\gamma$ verifies Assumption 1 and $Y_k$ has $p_\infty$ as its invariant law. By Theorem 1, we get that there exist $\gamma_2 > 0$ and $B_1 \ge 0$ such that $\forall \gamma \in (0, \gamma_2]$

$$
\begin{align*}
\mathbf{W}_p(\hat{\mu}_{\gamma}, p_{\infty}) &\le B_1 \left( \mathbb{E}_{Y \sim \hat{\mu}_{\gamma}} (\|\hat{b}^{\gamma}(Y) - b^{\gamma}(Y)\|^2) \right)^{\frac{1}{2p}} \\
&\le 2^{\frac{1}{2p}} B_1 \left( \mathbb{E}_{Y \sim \hat{\mu}_{\gamma}} (\|G_{\gamma}(Y) - \nabla g^{\gamma}(Y)\|^2 + \|\nabla f(Y - \gamma G_{\gamma}(Y)) - \nabla f(Y - \gamma \nabla g^{\gamma}(Y))\|^2) \right)^{\frac{1}{2p}} \\
&\le 2^{\frac{1}{2p}} B_1 \left( \mathbb{E}_{Y \sim \hat{\mu}_{\gamma}} (\|G_{\gamma}(Y) - \nabla g^{\gamma}(Y)\|^2 + L_f^2 \gamma_2 \|G_{\gamma}(Y) - \nabla g^{\gamma}(Y)\|^2) \right)^{\frac{1}{2p}} \\
&\le 2^{\frac{1}{2p}} B_1 (1 + L_f^2 \gamma_2)^{\frac{1}{2p}} \left( \mathbb{E}_{Y \sim \hat{\mu}_{\gamma}} (\|\gamma^{-1} (\text{Prox}_{\gamma g}(Y) - S_{\gamma}(Y))\|^2) \right)^{\frac{1}{2p}} \\
&\le 2^{\frac{1}{2p}} B_1 (1 + L_f^2 \gamma_2)^{\frac{1}{2p}} \gamma^{-\frac{1}{p}} \|\text{Prox}_{\gamma g} - S_{\gamma}\|_{\ell_2(\hat{\mu}_{\gamma})}^{\frac{1}{2p}}. \tag{96}
\end{align*}
$$

By combining equation (86), equation (87), equation (95) and equation (96) and the triangle inequality, we get that

$$
\begin{align*}
\mathbf{W}_p(p_{\hat{Y}_k}, \mu_\gamma) &\le B_0 r^{k\gamma} + 2^{\frac{1}{2p}} B_1 (1 + L_f^2 \gamma_2)^{\frac{1}{2p}} \gamma^{-\frac{1}{p}} \|\mathrm{Prox}_{\gamma g} - S_\gamma\|_{\ell_2(\hat{\mu}_\gamma)}^{\frac{1}{2p}} + C_3 \gamma^{\frac{1}{p}} + C_4 \gamma^{\frac{1}{2p}} \\
&\le B_0 r^{k\gamma} + 2^{\frac{1}{2p}} B_1 (1 + L_f^2 \gamma_2)^{\frac{1}{2p}} \gamma^{-\frac{1}{p}} \|\mathrm{Prox}_{\gamma g} - S_\gamma\|_{\ell_2(\hat{\mu}_\gamma)}^{\frac{1}{2p}} + C_3 \gamma_2^{\frac{1}{2p}} \gamma^{\frac{1}{2p}} + C_4 \gamma^{\frac{1}{2p}} \\
&\le B_0 r^{k\gamma} + B_2 \gamma^{-\frac{1}{p}} \|\mathrm{Prox}_{\gamma g} - S_\gamma\|_{\ell_2(\hat{\mu}_\gamma)}^{\frac{1}{2p}} + B_3 \gamma^{\frac{1}{2p}},
\end{align*}
$$

$$
\text{with } B_2 = 2^{\frac{1}{2p}} B_1 (1 + L_f^2 \gamma_2)^{\frac{1}{2p}} \text{ and } B_3 = C_3 \gamma_2^{\frac{1}{2p}} + C_4.
$$

Now, we prove that for $\gamma$ small enough, $G_\gamma$ is Lipschitz. By Assumption 1(i), $\hat{b}^\gamma$ is $L$-Lipschitz and by Assumption 2(i) $\nabla f$ is $L_f$-Lipschitz, so for $\gamma \in (0, \gamma_2]$ and $x, y \in \mathbb{R}^d$, we have

$$
\begin{align*}
& \|\hat{b}^{\gamma}(x) - \hat{b}^{\gamma}(y)\| \le L \|x - y\| \\
& \|G_{\gamma}(x) - G_{\gamma}(y)\| + \nabla f(x - \gamma G_{\gamma}(x)) - \nabla f(y - \gamma G_{\gamma}(y)) \| \le L \|x - y\| \\
& \|G_{\gamma}(x) - G_{\gamma}(y)\| - \|\nabla f(x - \gamma G_{\gamma}(x)) - \nabla f(y - \gamma G_{\gamma}(y))\| \le L \|x - y\| \\
& \|G_{\gamma}(x) - G_{\gamma}(y)\| - L_f (\|x - y\| + \gamma \|G_{\gamma}(x)G_{\gamma}(y)\|) \le L \|x - y\| \\
& (1 - \gamma L_f) \|G_{\gamma}(x) - G_{\gamma}(y)\| \le (L + L_f) \|x - y\|.
\end{align*}
$$
For $\gamma \le \frac{1}{2L_f}$, we get

$$ \|G_{\gamma}(x) - G_{\gamma}(y)\| \le 2(L + L_f)\|x - y\|. \quad (97) $$

So $G_{\gamma}$ is $L_G = 2(L + L_f)$-Lipschitz. We denote $\gamma_3 = \min(\gamma_2, \frac{1}{2L_f})$.

We have that $p_{\hat{X}_k} = S_{\gamma}\#p_{\hat{Y}_k}$ and $S_{\gamma} = I_d - \gamma G_{\gamma}$ is $(1 + \gamma_3 L_G)$ Lipschitz. Hence we obtain

$$ \begin{aligned} \mathbf{W}_p(p_{\hat{X}_k}, \nu_{\gamma}) &= \min_{\beta \in \Lambda} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|S_{\gamma}(x_1) - S_{\gamma}(x_2)\|^p d\beta(x_1, x_2) \\ &\le (1 + \gamma_3 L_G) \min_{\beta \in \Lambda} \int_{\mathbb{R}^d \times \mathbb{R}^d} \|x_1 - x_2\| d\beta(x_1, x_2) \\ &\le (1 + \gamma_3 L_G) \mathbf{W}_p(p_{\hat{Y}_k}, \mu_{\gamma}), \end{aligned} $$

with $\Lambda$ the set of distribution on $\mathbb{R}^d \times \mathbb{R}^d$ with marginals $p_{\hat{Y}_k}$ and $\mu_{\gamma}$. This proves the desired result with $\hat{\gamma} = \gamma_3$, $C_5 = (1 + \gamma_3 L_G)B_0$, $C_6 = (1 + \gamma_3 L_G)B_2$ and $C_6 = (1 + \gamma_3 L_G)B_3$. $\square$