# Coupling Analysis of the Asymptotic Behaviour of a Primal-Dual Langevin
Algorithm

Martin Burger  Helmholtz Imaging, Deutsches Elektronen-Synchrotron DESY,
Notkestr. 85, 22607 Hamburg, Germany. Fachbereich Mathematik, Universität
Hamburg, Bundesstrasse 55, Hamburg, 20146, Germany.  Matthias J. Ehrhardt
Department of Mathematical Sciences, University of Bath, Claverton Down, Bath
BA2 7AY, United Kingdom.  Lorenz Kuger  Helmholtz Imaging, Deutsches
Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany. Department
Mathematik, Friedrich-Alexander-Universität Erlangen-Nürnberg, Cauerstr. 11,
91058 Erlangen, Germany.  Corresponding author –<EMAIL_ADDRESS>Lukas
Weigand  Helmholtz Imaging, Deutsches Elektronen-Synchrotron DESY, Notkestr.
85, 22607 Hamburg, Germany.

###### Abstract

We analyze a recently proposed algorithm for the problem of sampling from
probability distributions $\mu^{\ast}$ in $\mathbb{R}^{d}$ with a Lebesgue
density of the form $\mu^{\ast}(x)\propto\exp(-f(Kx)-g(x))$, where $K$ is a
linear operator and $f,g$ convex and non-smooth. The algorithm is a
generalization of the primal-dual hybrid gradient (PDHG) convex optimization
algorithm to a sampling scheme. We analyze the method’s continuous time limit,
an SDE in the joint primal-dual variable. We give mild conditions under which
the corresponding Fokker-Planck equation converges to a unique stationary
state, which however does not concentrate in the dual variable and
consequently does not have $\mu^{\ast}$ as its primal marginal. Under a
smoothness assumption on $f$, we show that the scheme converges to the purely
primal overdamped Langevin diffusion in the limit of small primal and large
dual step sizes. We further prove that the target can never be the primal
marginal of the invariant solution for any modification of the SDE with space-
homogeneous diffusion coefficient. A correction with inhomogeneous diffusion
coefficient and the correct invariant solution is proposed, but the scheme
requires the same smoothness assumptions on $f$ and is numerically inferior to
overdamped Langevin diffusion. We demonstrate our findings numerically first
on small-scale examples in which we can exactly verify the theoretical
results, and subsequently on typical examples of larger scale from Bayesian
imaging inverse problems.

##  1 Introduction

We analyze a recently proposed algorithm for sampling from a probability
distribution on $\mathcal{X}\coloneqq\mathbb{R}^{d}$ with log-concave, non-
smooth density. In particular, we are interested in sampling from target
distributions with Lebesgue densities of the form

$\mu^{\ast}(x)=\frac{\exp(-f(Kx)-g(x))}{\int\exp(-f(Kw)-g(w))\,\mathrm{d}w},$
(1.1)

where we assume that the integral in the denominator exists and is finite. The
operator $K:\mathcal{X}\to\mathbb{R}^{m}\eqqcolon\mathcal{Y}$ is assumed to be
linear and $f:\mathcal{Y}\to\mathbb{R}\cup\\{\infty\\}$,
$g:\mathcal{X}\to\mathbb{R}\cup\\{\infty\\}$ will be called potential terms
and are both assumed convex and lower semi-continuous but not necessarily
differentiable. We will denote the full potential by $f\circ K+g\eqqcolon
h:\mathcal{X}\to\mathbb{R}\cup\\{\infty\\}$.

A standard class of algorithms to draw samples in Bayesian inference are
Markov chain Monte Carlo (MCMC) methods. Compared to other types of MCMC
methods, sampling schemes based on Langevin diffusion processes have proved
efficient in inverse problems and imaging applications due to their preferable
scaling behaviour to high dimensions [33, 7, 17]. MCMC sampling algorithms
based on Langevin diffusion typically perform a time discretisation of the
stochastic differential equation (SDE)

$\mathrm{X}_{t}=\nabla\log\mu^{\ast}(X_{t})\,\mathrm{d}t+\sqrt{2}\,\mathrm{d}W_{t},$
(1.2)

where $W_{t}$ is standard $\mathbb{R}^{d}$-valued Brownian motion. Assuming
$\mu^{\ast}$ satisfies a log-Sobolev inequality, it is the unique invariant
distribution of the SDE and the law of $X_{t}$ will converge to $\mu^{\ast}$
as $t$ grows. The standard Euler-Maruyama discretization of (1.2) takes the
form

$\left\\{\begin{array}[]{l}\xi^{n+1}\sim\mathrm{N}(0,I_{d})\\\
X^{n+1}=X^{n}+\tau\nabla\log\mu^{\ast}(X_{k})+\sqrt{2\tau}\xi^{n+1},\end{array}\right.$
(1.3)

this iteration is also known as unadjusted Langevin algorithm (ULA). ULA and
also similar implicit or semi-implicit time discretizations introduce a bias
between the law of the samples and the target measure, hence many works have
been concerned with the characterization of the invariant distribution of the
Markov chain $(X^{n})$ and its bias [33, 23, 21, 15]. Further questions
typically concern preconditioning and convergence speed of the Markov chain
and the behaviour of Metropolis-Hastings correction steps to overcome the bias
and draw unbiased samples from $\mu^{\ast}$.

In this work, we are concerned with a new type of Langevin sampling algorithm
which has been proposed as Unadjusted Langevin Primal-Dual Algorithm (ULPDA)
in [10] and has since been mentioned in [5, 3]. Unlike ULA or its variants,
ULPDA can not be interpreted as a time discretization of (1.2). The algorithm
consists in carrying out steps of a primal-dual optimization method, where the
same type of stochasticity as in unadjusted Langevin algorithms is added to
the primal variable in each step. Its motivation is hence inspired by the
optimization view: Just as (1.3) can be interpreted as a gradient descent with
an additive stochastic component, one might be able to obtain more general
sampling algorithms by adding the right amount of stochasticity to more
general classes of first-order optimization methods. This is, of course, not a
rigorous argument for convergence, and until now, to the best of the authors’
knowledge, no convergence theory for ULPDA has been provided. Indeed, the
situation seems to be more intricate: When we run ULPDA on a simple toy
example, the distribution of the samples seems to depend more critically on
the choice of step size hyperparameters $\tau,\sigma$ of the algorithm than it
is known from ULA, see Figure 1.

Figure 1: In a simple toy example, we choose $f,g:\mathbb{R}\to\mathbb{R}$ to
be quadratic functions, $K>0$ a real number, so that the target $\mu^{\ast}$
is Gaussian with the density depicted in orange. We run ULPDA with very small
primal step size $\tau$ and sufficient burn-in period and repeat the
experiment for three different ratios $\sigma/\tau$. The gray histograms show
the distribution of samples produced by ULPDA. For comparison, the right
figure is produced by ULA with the same primal step size $\tau$.

###  1.1 Related Works

There is a vast literature on unadjusted Langevin algorithms, where we use the
term to refer not only to (1.3), but any time discretization of (1.2). The
analysis is based on the fact that the target distribution is the stationary
distribution of the SDE (1.2) and that the dynamics converge to the target in
continuous time [33]. This is based on the fact that the dynamics can be
understood as a gradient flow of relative entropy with respect to the target
in Wasserstein space [32, 21]. Standard convergence results on ULA typically
require a Lipschitz continuous gradient of the potential in order to bound the
error introduced by time discretization [15, 17, 16].

Since non-smooth densities do not satisfy this gradient Lipschitz condition,
many works have treated the extension to the non-smooth case. Generally, the
approaches can be divided into two classes. The first regularizes non-smooth
terms, resulting in a modified, smooth target and applies typical ULA to this
modified distribution. The challenge then typically consists in bounding the
terms between the target and the modified distribution, and numerically in
trading off the introduced bias and convergence speed via the smoothing
parameter [25, 7, 20, 8, 1]. The second approach consists in finding different
discretizations of (1.2), i.e. implicit or semi-implicit schemes, that are
implementable and provably convergent. These schemes are usually referred to
as proximal or implicit Langevin methods, see e.g. [25, 18, 14, 4, 2]. A work
of particular relevance for our setting is [3], since the authors assume the
same type of splitting of the potential $h=f\circ K+g$ as we do. Their method
‘Prox-Sub’ can be understood as a first-order forward-backward discretization
of (1.2), with the forward step in $f\circ K$ realized as a subgradient since
$f$ is assumed non-smooth. We will see that this method is intimately
connected to the primal-dual sampling algorithm since it forms the limiting
case in which the primal-dual stationary distribution can be related to the
target. In [10], the authors proposed the primal-dual scheme that we consider
here, in order to generalize from the class of methods that employ purely
primal discretizations of (1.2). A partial dualization of the potential is
also considered in [1], although their analysis still requires the relevant
part $f$ of the potential to be smooth. We also refer to [5] for a recent
comparison of current methods applicable in the non-smooth case.

Constraints on the state space can be viewed as a special case of non-smooth
potentials, where the constraint is enforced by adding an indicator function
to the potential. Several works have considered this special case usually by
incorporating a proximal (which amounts to a projection in the case of an
indicator) or a reflection step [22, 19, 6]. Mirroring the Langevin dynamics
can be another approach to precondition the dynamics as well as to enforce
geometrical constraints [30, 11, 9, 8].

###  1.2 Contributions

  * •

For the continuous time limit of the iteration, we prove a contraction result
using a coupling argument, which makes the dynamics asymptotically stable and
guarantees convergence to the corresponding stationary solution.

  * •

We show that the marginal with respect to the primal variable of the
stationary solution in continuous time is not the target distribution. This
can be interpreted as being a consequence of the dynamics not concentrating in
the dual variable, although the diffusion coefficient is degenerate with the
stochastic term only added to the primal variable.

  * •

We consider two potential frameworks for “correcting” the bias between
continuous time stationary and target distribution. A proposed modification of
the dynamics with non-homogeneous diffusion coefficient can be employed to
force the target to become the stationary solution. Secondly, we consider the
unmodified scheme, but with the ratio of primal and dual step sizes in the
algorithm growing to infinity. In the limit, the stationary distribution of
the unmodified scheme is shown to converge to the target under regularity
assumptions on the dualized potential.

  * •

In discrete time, we show stability of the algorithm in the sense that
Wasserstein distances between two different instances are bounded if the
potentials are convex. Under stronger regularity assumptions on $f^{\ast},g$,
there exists a stationary solution to which the distribution of the iterates
converges.

The manuscript is organized as follows. In Section 2, we give the considered
sampling procedure based on primal-dual optimization and derive its continuous
time limit SDE. Since this limiting SDE differs from (1.2), we devote Section
3 to the analysis of the continuous time setting. We derive the mean-field
limit Fokker-Planck equation and show contraction properties in Section 3.1.
In Section 3.2, we prove a central ‘negative’ result, i.e. that the target can
not be the stationary distribution for any homogeneous diffusion coefficient
and illustrate the (non-)concentration of the mean-field limit, which is
crucial to understand how its stationary state relates to the target
distribution. This relation is proven in Section 3.3, which allows to bound
the bias due to the dualization in Wasserstein distance. We further show that
the PDE can be modified to force the target to be its stationary distribution
by employing a non-homogeneous diffusion coefficient. By transferring some of
the time continuous arguments to the discretized algorithm and using typical
computations known from primal-dual optimization methods, we show stability of
the algorithm and convergence to a stationary state in Section 4. We present
numerical illustrations of our most important findings in Section 5

##  2 The Primal-Dual Sampling Scheme

Although the algorithmic approach is not limited to this setting, in the
context of Bayesian inference, the target distribution $\mu^{\ast}$ is
typically a posterior distribution. Such a posterior target is defined through
the likelihood $p(z|x)$ of observing $z$ given the true parameter $x$,
combined with a prior distribution $p(x)$. It is given by

$p(x|z)=\frac{1}{Z}p(z|x)p(x),$

where the normalization constant $Z=\int p(z|w)p(w)\,\mathrm{d}w$ is typically
intractable in high dimensions. As an example class of problems that we
consider, in Bayesian inference for imaging inverse problems the distributions
are typically high-dimensional and involve non-smooth terms. This is because
model-based priors on image distributions often assume sparse coefficients in
some basis, the recovery of which corresponds to minimization of non-smooth
functions [31]. Therefore, many relevant posteriors $p(x|z)$ fall into the
class of targets given in (1.1). We start by explaining the motivation for the
primal-dual sampling from the point of view of convex, non-smooth
optimization.

###  2.1 Motivation for ULPDA

A common task in Bayesian inverse problems is the computation of the maximum a
posteriori estimate (MAP), i.e. the mode of the distribution
$\mu^{\ast}=p(\cdot|z)$. This task forms the connection between Bayesian
formulations of inverse problems and variational regularization approaches,
since it amounts to the problem of computing the solution to the following
minimization problem

$\operatorname*{arg\,min}_{x}f(Kx)+g(x).$ (2.1)

In order to solve (2.1), depending on the structure of $f$, $K$ and $g$ it can
be helpful to employ techniques from convex duality on parts of the potential
$h=f\circ K+g$. A popular approach consists in reformulating the variational
problem (2.1) to the saddle point problem

$\displaystyle\min_{x}h(x)$ $\displaystyle=\min_{x}(g(x)+\max_{y}(\langle
Kx,y\rangle-f^{*}(y))$ $\displaystyle=\min_{x}\max_{y}s(x,y)$ (2.2)

with $s(x,y):=g(x)+\langle Kx,y\rangle-f^{*}(y)$. Here $f^{\ast}$ denotes the
convex conjugate of $f$ defined by

$f^{\ast}(y):=\sup_{x}\left\\{\langle x,y\rangle-f(x)\right\\}.$

Under mild assumptions on $g$ and $f^{\ast}$, one can prove existence of a
saddle point $(x^{\ast},y^{\ast})$ (see section 19 in [28] for details) the
primal part $x^{\ast}$ of which is a solution to the original problem (2.1).
Instead of minimizing $h$, one can then employ primal-dual type algorithms to
solve the saddle point problem (2.2). For the context of imaging inverse
problems, we refer to [24] for an overview.

The primal-dual optimization method (primal descent and dual ascent with
extrapolation) for (2.2) takes the form

$\left\\{\begin{array}[]{l}Y^{n+1}=\operatorname{prox}_{\sigma
f^{\ast}}\left(Y^{n}+\sigma KX_{\theta}^{n}\right)\\\
X^{n+1}=\operatorname{prox}_{\tau g}(X^{n}-\tau K^{T}Y^{n})\\\
X_{\theta}^{n+1}=X^{n+1}+\theta(X^{n+1}-X^{n})\end{array}\right.$ (2.3)

with step sizes $\tau$ and $\sigma$ and a relaxation parameter $\theta$ chosen
accordingly. Herein, the proximal operator for a convex functional $\varphi$
and a step size $\gamma>0$ is defined by

$\operatorname{prox}_{\gamma\varphi}(z)\coloneqq(id+\gamma\partial\varphi)^{-1}(z)=\operatorname*{arg\,min}_{w}\left\\{\varphi(w)+\frac{1}{2\gamma}\lVert
w-z\rVert_{2}^{2}\right\\}.$

The minimum is attained by a unique point, and the application of the proximal
operator can be viewed as an implicit time step in the gradient flow of the
functional $\varphi$. For convergence theory of the scheme (2.3) and possible
accelerations we refer to [29, 24].

In order to draw the parallels between computing the MAP of $\mu^{\ast}$ and
drawing samples from $\mu^{\ast}$, we want to recall the standard unadjusted
Langevin algorithm (1.3). Under appropriate conditions on $h$ and the step
size $\tau$, ULA converges to a stationary distribution whose distance to the
target $\mu^{\ast}$ can be bounded by a factor depending on $h$ and the
considered metric on the space of probability distributions [15]. ULA can be
interpreted as a scheme alternating a gradient descent step in $h$ and a
diffusion step (or, equivalently, a time step along the gradient flow of
entropy in the Wasserstein space). In particular, the additional diffusion
step ensures the samples have the correct covariance (up to a bias scaling
linearly in the step size $\tau$).

It appears natural to take a similar perspective for primal-dual optimization
methods and add a diffusion step at some point in the iteration. Such a scheme
was proposed (without further theoretical analysis) under the name unadjusted
Langevin primal-dual algorithm (ULPDA) in [10]. The iteration takes the form

$\left\\{\begin{array}[]{l}Y^{n+1}=\operatorname{prox}_{\sigma
f^{\ast}}\left(Y^{n}+\sigma KX_{\theta}^{n}\right)\\\
\xi^{n+1}\sim\mathrm{N}(0,I_{d})\\\ X^{n+1}=\operatorname{prox}_{\tau
g}(X^{n}-\tau K^{T}Y^{n})+\sqrt{2\tau}\xi^{n+1}\\\
X_{\theta}^{n+1}=X^{n+1}+\theta(X^{n+1}-X^{n}),\end{array}\right.$ (2.4)

where the diffusion step is added in the primal variable $X_{k}$ and the
overrelaxation parameter $\theta$ was set to $\theta=1$ in [10] in the spirit
of primal-dual optimization methods. Note that the stochastic term in the
update for $X^{n+1}$ may be added in the argument of the proximal operator (as
in [10]) or after applying the proximal operator (as in 2.4). For the
considerations in continuous time, this makes no difference, however it will
simplify some of our analysis in discrete time based on couplings in Section
4.

###  2.2 Continuous time perspective

In the analysis of ULA’s convergence behaviour, it is helpful to view the
update step as an Euler-Maruyama time discretization of the overdamped
Langevin diffusion (1.2). Convergence results for ULA are based on the fact
that the target measure $\mu^{\ast}$ is the unique invariant distribution of
the Markov process $X_{t}$.

In order to analyze ULPDA, we want to define a similar time-continuous limit
of (2.4). For that, we assume that the step sizes $\tau,\sigma$ are coupled
via a parameter $\lim_{\tau,\sigma\to 0}\frac{\sigma}{\tau}=\lambda$. Taking
both step sizes to zero at this ratio, the formal limiting Markov process
satisfies the diffusion equation

$\left\\{\begin{array}[]{l}\,\mathrm{d}X_{t}\in-(\partial
g(X_{t})+K^{*}Y_{t})\,\mathrm{d}t+\sqrt{2}\,\mathrm{d}W_{t}\\\
\,\mathrm{d}Y_{t}\in-\lambda(\partial
f^{*}(Y_{t})-KX_{t})\,\mathrm{d}t,\end{array}\right.$ (2.5)

which we can write in a joint variable $Z_{t}=(X_{t},Y_{t})$ as

$\,\mathrm{d}Z_{t}\in-
A_{\text{cPD}}(Z_{t})\,\mathrm{d}t+B_{\text{cPD}}\,\mathrm{d}W_{t}.$ (2.6)

where we introduced the continuous-time primal-dual coefficients

$A_{\text{cPD}}(Z_{t})=\begin{pmatrix}\partial g&K^{T}\\\ -\lambda
K&\lambda\partial f^{\ast}\end{pmatrix},\qquad
B_{\text{cPD}}=\begin{pmatrix}\sqrt{2}I\\\ 0\end{pmatrix},$ (2.7)

with $I$ the identity in $\mathbb{R}^{d}$. Note in particular that the
diffusion coefficient
$\Sigma_{\text{cPD}}=\frac{1}{2}B_{\text{PD}}B_{\text{cPD}}^{T}$ is only
positive semi-definite (or degenerate) with an $m$-dimensional null space,
where $m$ is the dimension of the dual variable $Y$. As we will see in more
detail, the corresponding mean-field limit PDE is not elliptic anymore, which
complicates analysis of existence and uniqueness of solutions.

###  2.3 Spaces and notation

We denote by $\mathcal{P}(\mathbb{R}^{d})$ the space of probability measures
on $\mathbb{R}^{d}$, and by $\mathcal{P}_{p}(\mathbb{R}^{d})$ the subspace of
distributions with finite $p$-th moment. Given a norm $\lVert\cdot\rVert$ on
$\mathbb{R}^{d}$, a transport distance on $\mathcal{P}_{2}(\mathbb{R}^{d})$
can be defined by

$\mathcal{T}(\mu,\tilde{\mu}):=\inf_{\mathcal{M}\in\Gamma(\mu,\tilde{\mu})}\left(\mathbb{E}_{(x,\tilde{x})\sim\mathcal{M}}\left[\lVert
x-\tilde{x}\rVert^{2}\right]\right)^{1/2}.$

where $\Gamma(\mu,\tilde{\mu})$ denotes the set of all couplings with
marginals $\mu$ and $\tilde{\mu}$

$\Gamma(\mu,\tilde{\mu})=\left\\{\mathcal{M}\in\mathcal{P}(\mathbb{R}^{n}\times\mathbb{R}^{n})~{}:~{}\mu=\int_{\mathbb{R}^{n}}\mathcal{M}(\cdot,\tilde{x})\,d\tilde{x}~{}\textrm{and}~{}\tilde{\mu}=\int_{\mathcal{X}}\mathcal{M}(x,\cdot)\,dx\right\\}.$

$\mathcal{T}$ defines a metric on $\mathcal{P}_{2}(\mathbb{R}^{d})$. If
$\lVert\cdot\rVert$ is the Euclidean norm, then $\mathcal{T}$ is the
2-Wasserstein distance $\mathcal{W}_{2}$.

In the following, we will denote the state spaces $\mathcal{X}=\mathbb{R}^{d}$
for the primal and $\mathcal{Y}=\mathbb{R}^{m}$ for the dual variable. In
order to simplify notation, we introduce a weighted inner product and norm on
the product space $\mathcal{Z}:=\mathcal{X}\times\mathcal{Y}$. For a scalar
$a>0$, we define $\langle\cdot,\cdot\rangle_{a}$ on $\mathcal{X}$ (and
analogously on $\mathcal{Y}$) by $\langle x,\tilde{x}\rangle_{a}:=\langle
ax,\tilde{x}\rangle$ for all $x,\tilde{x}\in\mathcal{X}$, where
$\langle\cdot,\cdot\rangle$ is the Euclidean inner product on $\mathcal{X}$.
Further, denote by $\lVert\cdot\rVert_{a}$ the weighted norm $\lVert
x\rVert_{a}=\sqrt{\langle x,x\rangle_{a}}$. Given $a,b>0$, we define a
weighted inner product on the product space $\cal Z$ by

$\langle z,\tilde{z}\rangle_{a,b}=\langle x,\tilde{x}\rangle_{a}+\langle
y,\tilde{y}\rangle_{b},\quad
z=(x,y),\tilde{z}=(\tilde{x},\tilde{y})\in\mathcal{Z}.$

Accordingly, we define a weighted norm by

$\lVert z\rVert_{a,b}:=\sqrt{\langle z,z\rangle_{a,b}}=\sqrt{\lVert
x\rVert_{a}^{2}+\lVert y\rVert_{b}^{2}},\quad z=(x,y)\in\mathcal{Z}.$

On the space $\mathcal{P}(\mathcal{Z})$, we will also denote
$\mathcal{T}_{a,b}$ for the transport distance induced by the norm
$\lVert\cdot\rVert_{a,b}$. Note that $\mathcal{T}_{a,b}$ is a metric on
$\mathcal{P}_{2}(\mathcal{Z})$ and $\mathcal{T}_{1,1}=\mathcal{W}_{2}$. For
every $a,b>0$, $\mathcal{T}_{a,b}$ is equivalent to $\mathcal{W}_{2}$, since
the fact $\sqrt{\min\\{a,b\\}}\,d_{1,1}(z,\tilde{z})\leq
d_{a,b}(z,\tilde{z})\leq\sqrt{\max\\{a,b\\}}\,d_{1,1}(z,\tilde{z})$ directly
implies

$\sqrt{\min\\{a,b\\}}\,\mathcal{W}_{2}\leq\mathcal{T}_{a,b}\leq\sqrt{\max\\{a,b\\}}\,\mathcal{W}_{2}.$

##  3 The Stationary Solution in Continuous Time

We start by proving the existence of a unique stationary solution of the time-
continuous limit (2.6) of the ULPDA algorithm. We further show that the
corresponding primal-dual Fokker-Planck equation is contractive, implying
convergence of any solution path to the stationary distribution. The
contraction result is formulated for a general class of diffusion SDEs. It
follows a standard calculation from the literature for overdamped Langevin
diffusion and is generalized to the setting of an SDE with any strongly
monotone operator as non-homogeneous drift coefficient.

###  3.1 Fokker-Planck Equation and Existence of a Stationary Solution

####  3.1.1 Contraction of Overdamped Langevin Diffusion via Coupling Argument

Suppose for now that the potential $h:\mathcal{X}\to\mathbb{R}$ is
differentiable, that $\nabla h$ is $L$-Lipschitz-continuous and
$\omega_{h}$-strongly convex. We repeat here an instructive calculation that
has been carried out, e.g., in [13, 17] to show stability of the overdamped
Langevin diffusion process

$dX_{t}=-\nabla h(X_{t})dt+\sqrt{2}dW_{t}.$ (3.1)

The distribution $\mu_{t}$ of the process $X_{t}$ in (3.1), denoted
$\mu_{t}=\operatorname{Law}(X_{t})$, satisfies the Fokker-Planck equation (or
forward Kolmogorov equation) [27] given by

$\frac{\partial\mu_{t}}{\partial t}=\nabla_{x}\cdot(\mu_{t}\nabla
h(x))+\Delta_{x}\mu_{t}.$

The target distribution $\mu^{\ast}$ is the unique stationary solution of the
Fokker-Planck equation. Furthermore, the PDE is contractive in Wasserstein
distance, which can be shown using the following coupling argument. We will
call two processes coupled if both solve (3.1) with the same Brownian motion.
Let $X_{t},\tilde{X}_{t}$ be coupled with $X_{0}\sim\mu_{0}$ and
$\tilde{X}_{0}\sim\tilde{\mu}_{0}$. Then the distribution $\mathcal{M}_{t}$ of
the joint variable $(X_{t},\tilde{X}_{t})$ satisfies the Fokker-Planck
equation

$\frac{\partial\mathcal{M}_{t}}{\partial
t}=\nabla_{x}\cdot(\mathcal{M}_{t}\nabla
h(x))+\nabla_{\tilde{x}}\cdot(\mathcal{M}_{t}\nabla
h(\tilde{x}))+(\nabla_{x}+\nabla_{\tilde{x}})\cdot(\nabla_{x}+\nabla_{\tilde{x}})\mathcal{M}_{t}.$
(3.2)

Note that we can rewrite this using
$(\nabla_{x}+\nabla_{\tilde{x}})\cdot(\nabla_{x}+\nabla_{\tilde{x}})=\Delta_{x+\tilde{x}}$
by a coordinate change $(x,\tilde{x})\mapsto(x+\tilde{x},x-\tilde{x})$. For
any distance function of the form $d(x,\tilde{x})=\rho(\lVert
x-\tilde{x}\rVert)$, we directly see $\Delta_{x+\tilde{x}}d(x,\tilde{x})=0$.
Using this and the property (3.2) of the joint distribution $\mathcal{M}_{t}$,
a contraction in Wasserstein distance can hence be shown in the following way.
We start by estimating

$\displaystyle\frac{1}{2}\frac{d}{dt}\mathbb{E}_{(x,\tilde{x})\sim\mathcal{M}_{t}}\left[\lVert
x-\tilde{x}\rVert^{2}\right]$ $\displaystyle=\frac{1}{2}\frac{d}{dt}\int\lVert
x-\tilde{x}\rVert^{2}\mathcal{M}_{t}(\mathrm{d}x,\mathrm{d}\tilde{x})$
$\displaystyle=\frac{1}{2}\int\lVert
x-\tilde{x}\rVert^{2}\left[\nabla_{x}\cdot(\mathcal{M}_{t}\nabla
h(x))+\nabla_{\tilde{x}}\cdot(\mathcal{M}_{t}\nabla
h(\tilde{x}))+\Delta_{x+\tilde{x}}\mathcal{M}_{t}\right]\,\mathrm{d}x\,\mathrm{d}\tilde{x}$
$\displaystyle=-\int(x-\tilde{x})\cdot(\nabla h(x)-\nabla
h(\tilde{x}))\mathcal{M}_{t}(\mathrm{d}x,\mathrm{d}\tilde{x})$
$\displaystyle\leq-\omega_{h}\mathbb{E}_{(x,\tilde{x})\sim\mathcal{M}_{t}}\left[\lVert
x-\tilde{x}\rVert^{2}\right],$

where we used that due to strong convexity of $h$, we have
$(x-\tilde{x})\cdot(\nabla h(x)-\nabla h(\tilde{x}))\geq\omega_{h}\lVert
x-\tilde{x}\rVert^{2}$. Applying now Gronwall’s lemma, we obtain

$\mathbb{E}_{(x,\tilde{x})\sim\mathcal{M}_{t}}\left[\lVert
x-\tilde{x}\rVert^{2}\right]\leq
e^{-2\omega_{h}t}\mathbb{E}_{(x,\tilde{x})\sim\mathcal{M}_{0}}\left[\lVert
x-\tilde{x}\rVert^{2}\right].$

The left side is bounded from below by
$\mathcal{W}_{2}^{2}(\mu_{t},\tilde{\mu}_{t})$ by definition. By taking the
infimum over all couplings $\mathcal{M}_{0}$ (with marginals $\mu_{0}$ and
$\tilde{\mu}_{0}$) on the right hand side, we obtain

$\mathcal{W}_{2}(\mu_{t},\tilde{\mu}_{t})\leq
e^{-\omega_{h}t}\mathcal{W}_{2}(\mu_{0},\tilde{\mu}_{0}).$

In particular, if we choose one of the processes to be initialized at the
target distribution as $\tilde{\mu}_{0}=\mu_{\ast}$, then the distribution of
$\tilde{X}_{t}$ is stationary at $\tilde{\mu}_{t}=\mu^{\ast}$ and by the above
we obtain

$\mathcal{W}_{2}(\mu_{t},\mu^{\ast})\leq
e^{-\omega_{h}t}\mathcal{W}_{2}(\mu_{0},\mu^{\ast}).$

####  3.1.2 Contraction and Stationary Solution of the Primal-Dual Diffusion
Equation

We can generalize the above argument allowing an analysis of the time-
continuous limit (2.6) of the primal dual sampling scheme (2.4). For that, we
consider $Z_{t}$ to satisfy a general stochastic differential inclusion with
non-homogeneous coefficients

$\mathrm{d}Z_{t}\in-A(Z_{t})\,\mathrm{d}t+B(Z_{t})\,\mathrm{d}W_{t}.$ (3.3)

Here $W_{t}$ is Brownian motion in $\mathbb{R}^{K}$ and
$B:\mathbb{R}^{d}\to\mathbb{R}^{d\times K}$. The drift $A:\mathbb{R}^{n}\to
2^{\mathbb{R}^{d}}$ is a set-valued monotone operator in the norm
$\lVert\cdot\rVert$, meaning

$\langle w-\tilde{w},z-\tilde{z}\rangle\geq 0,\qquad\forall w\in
A(z),\tilde{w}\in A(\tilde{z}).$

We say that $A$ is $\omega$-strongly monotone if the following stronger
condition is satisfied:

$\langle w-\tilde{w},z-\tilde{z}\rangle\geq\omega\lVert
z-\tilde{z}\rVert^{2},\qquad\forall w\in A(z),\tilde{w}\in A(\tilde{z}).$

######  Proposition 3.1.

Let $Z_{t}$ be a solution to (3.3) and denote by $\pi_{t}$ the density w.r.t.
the Lebesgue measure of the distribution of $Z_{t}$. Then $\pi_{t}$ satisfies
the Fokker-Planck equation

$\frac{\partial\pi_{t}}{\partial
t}=\nabla_{z}\cdot\left(w\pi_{t}+\Sigma(z)\nabla_{z}\pi_{t}\right),\qquad w\in
A(z).$ (3.4)

with $\Sigma=\frac{1}{2}BB^{T}\in\mathbb{R}^{n\times n}$.

Since we are repeatedly going to make use of coupling arguments, note that by
applying Proposition 3.1 to a coupled diffusion

$\begin{pmatrix}\mathrm{d}Z_{t}\\\
\mathrm{d}\tilde{Z}_{t}\end{pmatrix}\in-\begin{pmatrix}A(Z_{t})\\\
A(\tilde{Z}_{t})\end{pmatrix}\mathrm{d}t+\begin{pmatrix}B\\\
B\end{pmatrix}\mathrm{d}W_{t},$ (3.5)

we can also obtain a Fokker-Planck equation for the joint density of the
coupled processes.

######  Proposition 3.2.

Let $Z_{t},\tilde{Z}_{t}$ be two coupled processes that both solve (3.3) with
the same Brownian motion, i.e. $(Z_{t},\tilde{Z}_{t})$ solves (3.5). Then the
joint distribution $\Pi_{t}(z,\tilde{z})$ of $(Z_{t},\tilde{Z}_{t})$ satisfies
the Fokker-Planck equation

$\frac{\partial\Pi_{t}}{\partial
t}=\nabla_{z}\cdot\left(w\pi_{t}\right)+\nabla_{\tilde{z}}\cdot\left(\tilde{w}\Pi_{t}\right)+(\nabla_{z}+\nabla_{\tilde{z}})\cdot(\Sigma(\nabla_{z}+\nabla_{\tilde{z}})\Pi_{t}),\qquad
w\in A(z),\tilde{w}\in A(\tilde{z}).$ (3.6)

For showing the contractivity of overdamped Langevin dynamics in Section
3.1.1, the relevant requirement was the strong convexity of $h$, implying that
$\nabla h$ is a strongly monotone operator. The computation can be directly
generalized under a strong monotonicity assumption on $A$.

######  Lemma 3.3.

Suppose $A$ is $\omega$-strongly monotone. Then (3.3) is contractive in the
sense that, if $Z_{t},\tilde{Z}_{t}$ are two coupled processes that solve
(3.3) with the same Brownian motion, then for $t\geq 0$ their associated
distributions $\pi_{t},\tilde{\pi}_{t}$ satisfy

$\mathcal{W}_{2}(\pi_{t},\tilde{\pi}_{t})\leq e^{-\omega
t}\mathcal{W}_{2}(\pi_{0},\tilde{\pi}_{0}).$

###### Proof.

By employing (3.6), we obtain with the same reasoning as before

$\displaystyle\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\mathbb{E}_{(z,\tilde{z})\sim\Pi_{t}}\left[\lVert
z-\tilde{z}\rVert^{2}\right]$
$\displaystyle=\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int\lVert
z-\tilde{z}\rVert^{2}\Pi_{t}(\mathrm{d}z,\mathrm{d}\tilde{z})$
$\displaystyle=\frac{1}{2}\int\lVert
z-\tilde{z}\rVert^{2}\left(\nabla_{z}\cdot\left(w\Pi_{t}\right)+\nabla_{\tilde{z}}\cdot\left(\tilde{w}\Pi_{t}\right)+(\nabla_{z}+\nabla_{\tilde{z}})\cdot(\Sigma(\nabla_{z}+\nabla_{\tilde{z}})\Pi_{t})\right)\,\mathrm{d}z\,\mathrm{d}\tilde{z}$
$\displaystyle=-\int\left((z-\tilde{z})\cdot
w+(\tilde{z}-z)\cdot\tilde{w}\right)\Pi_{t}(\mathrm{d}z,\mathrm{d}\tilde{z})$
$\displaystyle\leq-\omega\mathbb{E}_{(z,\tilde{z})\sim\Pi_{t}}\left[\lVert
z-\tilde{z}\rVert^{2}\right]$

In the second line we used Proposition 3.2 and in the third line partial
integration together with the fact that
$(\nabla_{z}+\nabla_{\tilde{z}})(\lVert z-\tilde{z}\rVert^{2})=0$. We now
apply Gronwall’s Lemma and obtain

$\mathbb{E}_{(z,\tilde{z})\sim\Pi_{t}}\left[\lVert
z-\tilde{z}\rVert^{2}\right]\leq e^{-2\omega
t}\mathbb{E}_{(z,\tilde{z})\sim\Pi_{0}}\left[\lVert
z-\tilde{z}\rVert^{2}\right]$

which implies

$\mathcal{W}_{2}(\pi_{t},\tilde{\pi}_{t})\leq e^{-\omega
t}\mathbb{E}_{(z,\tilde{z})\sim\Pi_{0}}\left[\lVert
z-\tilde{z}\rVert^{2}\right].$

Taking the infimum over all couplings of $\pi_{0}$ and $\tilde{\pi}_{0}$ on
the right hand side concludes the proof. ∎

Consider now again the primal-dual type diffusion (2.6) and denote $\pi_{t}$
for the law of $Z_{t}$. By Proposition 3.1, $\pi_{t}$ solves the Fokker-Planck
equation

$\frac{\partial\pi_{t}}{\partial
t}=\nabla_{z}\cdot\left(w\pi_{t}+\Sigma_{\text{{cPD}}}\nabla_{z}\pi_{t}\right),\qquad
w\in A_{\text{{cPD}}}(z).$ (3.7)

with $\Sigma_{\text{{cPD}}}=\frac{1}{2}B_{\text{{cPD}}}B_{\text{{cPD}}}^{T}$
and $A_{\text{{cPD}}},B_{\text{{cPD}}}$ as in (2.7).

Since the drift coefficient $A_{\text{{cPD}}}$ is a monotone operator, the
primal-dual type diffusion falls into the class of general diffusion equations
considered in Lemma 3.3. Under strong convexity of $g,f^{\ast}$, the drift
coefficient is strongly monotone and we obtain the following stability result.

######  Theorem 3.4.

Assume that $g,f^{\ast}$ are $\omega_{g}$-strongly convex and
$\omega_{f^{\ast}}$-strongly convex, respectively, with
$\omega_{g},\omega_{f^{\ast}}>0$. Let $Z_{t},\tilde{Z}_{t}$ be two coupled
processes that solve the primal-dual diffusion equation (2.5) with
corresponding densities $\pi_{t}$ and $\tilde{\pi}_{t}$. Then it holds

$\mathcal{T}_{1,\lambda^{-1}}(\pi_{t},\tilde{\pi}_{t})\leq e^{-\omega
t}\mathcal{T}_{1,\lambda^{-1}}(\pi_{0},\tilde{\pi}_{0}),$

where $\omega=\min\\{\omega_{g},\lambda\omega_{f^{\ast}}\\}>0$. In particular,
the dynamics are contractive and there exists a unique stationary distribution
of the continuous time primal-dual diffusion process, denoted
$\pi^{\text{{stat}}}_{\text{{cPD}}}$, which solves the primal-dual Fokker-
Planck equation (3.7).

###### Proof.

We employ Lemma 3.3, hence the only part we need to show is the strong
monotonicity of the negative drift coefficient $A_{\text{{cPD}}}$ defined by

$A_{\text{{cPD}}}=\begin{pmatrix}\partial g&K^{T}\\\ -\lambda
K&\lambda\partial f^{\ast}\end{pmatrix}.$

We write $z=(x,y)\in\mathcal{X}\times\mathcal{X}=\mathcal{Z}$ for the joint
primal and dual variable and equip $\mathcal{Z}$ with the weighted inner
product $\langle\cdot,\cdot\rangle_{1,\lambda^{-1}}$. Let now $w\in
A_{\text{{cPD}}}(z)$ and $\tilde{w}\in A_{\text{{cPD}}}(\tilde{z})$. Then
there are $u\in\partial g(x)$ and $v\in\partial f^{\ast}(y)$ such that
$w=(u+K^{T}y,\lambda(v-Kx))$. Equivalently, if $\tilde{w}\in A(\tilde{z})$
then there exist $\tilde{u}\in\partial g(\tilde{x})$ and $\tilde{v}\in\partial
f^{\ast}(\tilde{y})$ with
$\tilde{w}=(\tilde{u}+K^{T}\tilde{y},\lambda(\tilde{v}-K\tilde{x}))$. It
follows that

$\displaystyle\langle w-\tilde{w},z-\tilde{z}\rangle_{1,\lambda^{-1}}$
$\displaystyle=\langle u-\tilde{u},x-\tilde{x}\rangle+\lambda\langle
v-\tilde{v},y-\tilde{y}\rangle_{\lambda^{-1}}$ $\displaystyle\qquad+\langle
K^{T}(y-\tilde{y}),x-\tilde{x}\rangle-\lambda\langle
K(x-\tilde{x}),y-\tilde{y}\rangle_{\lambda^{-1}}$ $\displaystyle=\langle
u-\tilde{u},x-\tilde{x}\rangle+\langle v-\tilde{v},y-\tilde{y}\rangle$
$\displaystyle\geq\omega_{g}\lVert
x-\tilde{x}\rVert^{2}+\omega_{f^{\ast}}\lVert y-\tilde{y}\rVert^{2}$
$\displaystyle\geq\omega\lVert z-\tilde{z}\rVert_{1,\lambda^{-1}}^{2}.$

∎

The stationary solution of (3.7) will be denoted by
$\pi^{\text{{stat}}}_{\text{{cPD}}}$ from here. We will also frequently refer
to its $x$-marginal as $\mu^{\text{{stat}}}_{\text{{cPD}}}$.

###  3.2 Non-Concentration of the Stationary Solution

Having established the existence of $\pi^{\text{{stat}}}_{\text{{cPD}}}$, we
can try to relate the algorithmic, time discretized scheme to it. However,
while overdamped Langevin dynamics is designed to have the target $\mu^{\ast}$
as its stationary solution, this is not the case in the primal-dual setting.
In general, $\mu^{\ast}$ is not the marginal distribution of $X$ when
$(X,Y)\sim\pi^{\text{{stat}}}_{\text{{cPD}}}$, i.e.
$\mu^{\ast}\neq\mu^{\text{{stat}}}_{\text{{cPD}}}$. We demonstrate this in a
simple example of a Gaussian target, where we can solve the Fokker-Planck
equation (3.7) in closed form.

######  Theorem 3.5.

Consider a one-dimensional setting with the target
$\mu^{\ast}(x)\propto\exp(-h(x))$ given by

$f(y)\coloneqq\frac{1}{2c_{f}}y^{2},\ g(x)\coloneqq\frac{1}{2c_{g}}x^{2},\
h(x)\coloneqq f(kx)+g(x)=\frac{k^{2}c_{g}+c_{f}}{2c_{f}c_{g}}x^{2},\qquad
x,y\in\mathbb{R},$

where $k\neq 0$ and $c_{g},c_{f}>0$. It holds
$f^{\ast}(y)=\frac{c_{f}}{2}y^{2}$ for all $y\in\mathbb{R}$. The primal-dual
flow (2.6) of the joint variable $Z_{t}=(X_{t},Y_{t})^{T}\in\mathbb{R}^{2}$
reads

$\mathrm{d}Z_{t}=-A_{\text{{cPD}}}Z_{t}\mathrm{d}t+B_{\text{{cPD}}}\mathrm{d}W_{t},\quad\text{with
}A_{\text{{cPD}}}=\begin{pmatrix}c_{g}^{-1}&k\\\ -\lambda k&\lambda
c_{f}\end{pmatrix},B_{\text{{cPD}}}=\begin{pmatrix}\sqrt{2}\\\ 0\end{pmatrix}$

where $W_{t}$ is one-dimensional Brownian motion and the weight $\lambda>0$
relates primal and dual time scales. Denoting
$E(t):=(e_{ij}(t))_{i,j}:=\exp(-tA_{\text{{cPD}}}^{T})$, the following
statements hold:

  1. (i)

For any (possibly singular) initialization $\pi_{0}$ such that
$Z_{0}\sim\pi_{0}$, the distribution
$\pi_{t}\coloneqq\operatorname{Law}(Z_{t})$ has a $C^{\infty}$-smooth Lebesgue
density for $t>0$.

  2. (ii)

If $\pi_{0}$ is Gaussian (including degenerate Gaussians with singular
covariance matrix), then $\pi_{t}$ is (non-degenerate) Gaussian for all $t>0$.
As $t\to\infty$, the solution $\pi_{t}$ converges to the Gaussian stationary
$\pi^{\text{{stat}}}$ with mean zero and covariance matrix

$\frac{1}{(c_{f}+k^{2}c_{g})(1+\lambda
c_{f}c_{g})}\begin{pmatrix}c_{g}(c_{f}+k^{2}c_{g})+\lambda
c_{f}^{2}c_{g}^{2}&k\lambda c_{f}c_{g}^{2}\\\ k\lambda
c_{f}c_{g}^{2}&k^{2}\lambda c_{g}^{2}\end{pmatrix}.$

The $x$-marginal $\mu^{\text{{stat}}}_{\text{{cPD}}}$ is not equal to the
target $\mu^{\ast}$. However, the limiting joint distribution as
$\lambda\to\infty$ is the degenerate Gaussian
$(id,f^{\prime}(k\,\cdot))_{\\#}\mu^{\ast}$ whose $x$-marginal is the target
$\mu^{\ast}$.

###### Proof.

From the diffusion equation, it follows that the distribution
$\pi_{t}(z)=\pi_{t}(x,y)$ of $Z_{t}$ satisfies the Fokker-Planck equation

$\partial_{t}\pi_{t}=\nabla_{z}\cdot(\pi_{t}A_{\text{{cPD}}}z)+\Delta_{x}\pi_{t},$

with an initial value $\pi_{0}$. The Cauchy problem can be solved in closed
form, note that a similar computation has been carried out in [34]. By taking
the Fourier transform in the variable $z$ (denoted $\mathcal{F}$), we obtain
the first order initial value problem for
$p_{t}(\zeta)=\mathcal{F}(\pi_{t})(\zeta)$ with $\zeta=(\xi,\omega)^{T}$:

$\left\\{\begin{array}[]{l}\partial_{t}p_{t}=-\zeta^{T}A_{\text{{cPD}}}\nabla_{\zeta}p_{t}-\xi^{2}p_{t}\\\
p_{0}(\xi,\omega)=\mathcal{F}\pi_{0}.\end{array}\right.$

The corresponding characteristic equations are

$\displaystyle\dot{\zeta}(t)$ $\displaystyle=A_{\text{{cPD}}}^{T}\zeta(t)$
$\displaystyle\dot{p}(t)$ $\displaystyle=-\xi^{2}(t)p(t),$

and they are solved by

$\displaystyle\zeta(t)$ $\displaystyle=\exp(tA_{\text{{cPD}}}^{T})\zeta_{0}$
$\displaystyle p(t)$
$\displaystyle=p_{0}(\zeta_{0})\exp\left(-\int_{0}^{t}\xi^{2}(s)\,\mathrm{d}s\right)$

By eliminating the initial value
$\zeta_{0}=\exp(-tA_{\text{{cPD}}}^{T})\zeta(t)$, the solution for the initial
value problem is

$p_{t}(\zeta)=p_{0}(\exp(-tA_{\text{{cPD}}}^{T})\zeta)\exp\left(-\frac{1}{2}\int_{0}^{t}(B_{\text{{cPD}}}^{T}\exp(-sA_{\text{{cPD}}}^{T})\zeta)^{2}\,\mathrm{d}s\right).$

The quadratic form inside the integral is positive semidefinite in the
variable $\zeta$. For $t>0$, it is also positive definite: Otherwise, there
would need to be a $\zeta$ such that
$B_{\text{{cPD}}}^{T}\exp(-sA_{\text{{cPD}}}^{T})\zeta=0$ for all $s\in[0,t]$,
which implies $B_{\text{{cPD}}}^{T}(A_{\text{{cPD}}}^{T})^{k}\zeta=0$ for all
$k\geq 0$. This is equivalent to $(0,1)^{T}$ being an element of an invariant
subspace of $A_{\text{{cPD}}}^{T}$, which is impossible since the eigenvectors
of $A_{\text{{cPD}}}^{T}$ are of the form

$\left(\frac{c_{g}^{-1}-\lambda
c_{f}}{2k}\pm\frac{1}{2k}\sqrt{(c_{g}^{-1}-\lambda c_{f})^{2}-4\lambda
k^{2}},1\right)^{T},$

neither of which is a multiple of $(0,1)^{T}$. From the positive definiteness
of the quadratic form it follows that $p_{t}$ is a Schwartz function for every
$t>0$. In particular, its inverse Fourier transform $\pi_{t}$ is
$C^{\infty}$-smooth in the variable $z$ for $t>0$, which proves (i).

If the initialization $\pi_{0}$ is Gaussian with mean $z_{0}=(x_{0},y_{0})$
and (possibly singular) covariance $\Sigma_{0}$, we can simply take the
inverse Fourier transform to obtain the solution $\pi_{t}$. It is worth
noting, however, that we do not even have to compute the full distribution
$\pi_{t}$ in order to analyze the $x$-marginal $\mu_{t}$ of $\pi_{t}$. Since
marginalizing out the dual variable $y$ in $\pi_{t}$ corresponds to evaluating
$p_{t}$ at $\zeta=(\xi,0)$ in Fourier space, we get, denoting by
$\mathcal{F}^{-1}_{\xi}$ the inverse Fourier transform in $\xi$ and
$E(t):=(e_{ij}(t))_{i,j}:=\exp(-tA_{\text{{cPD}}}^{T})$,

$\displaystyle\mu_{t}(x)$ $\displaystyle:=\int\pi_{t}(x,y)\,\mathrm{d}y$
$\displaystyle=\mathcal{F}^{-1}_{\xi}(p_{t}(\xi,0))(x)$
$\displaystyle\propto\int\exp\left(i\xi(x-z_{0}\cdot
e_{1}(t))-\frac{\xi^{2}}{2}\left(\lVert(e_{1}(t)\rVert_{\Sigma_{0}}^{2}+2\int_{0}^{t}e^{2}_{11}(s)\,\mathrm{d}s\right)\right)\,\mathrm{d}\xi$
$\displaystyle=\exp\left(-\frac{1}{2}(x-z_{0}\cdot
e_{1}(t))^{2}\left(\lVert(e_{1}(t)\rVert_{\Sigma_{0}}^{2}+2\int_{0}^{t}e^{2}_{11}(s)\,\mathrm{d}s\right)^{-1}\right),$

where we abbreviated $e_{1}(t)$ for the first column of $E(t)$. Hence, the
primal variable follows a Gaussian distribution at every time $t>0$, with mean
$z_{0}^{T}e_{1}(t)$ and variance
$\lVert(e_{1}(t)\rVert_{\Sigma_{0}}^{2}+2\int_{0}^{t}e^{2}_{11}(s)\,\mathrm{d}s$.
Due to the contraction result Theorem 3.4, we know that $\mu_{t}$ converges to
the stationary distribution $\mu^{\text{{stat}}}_{\text{{cPD}}}$ as
$t\to\infty$. Indeed, this holds here since $\lim_{t\to\infty}E(t)=0$, so
$\mathbb{E}_{X\sim\mu^{\text{{stat}}}_{\text{{cPD}}}}[X]=0$. The variance of
$\mu^{\text{{stat}}}_{\text{{cPD}}}$ is

$\operatorname{Var}_{X\sim\mu^{\text{{stat}}}_{\text{{cPD}}}}[X]=2\int_{0}^{\infty}e^{2}_{11}(s)\,\mathrm{d}s=\frac{c_{g}(c_{f}+k^{2}c_{g})+\lambda
c_{f}^{2}c_{g}^{2}}{(c_{f}+k^{2}c_{g})(1+\lambda c_{f}c_{g})},$

which does not coincide with the variance of the target distribution
$\mu^{\ast}$. The target is however approximated as $\lambda$ becomes large,
since
$\lim_{\lambda\to\infty}\operatorname{Var}[\mu^{\text{{stat}}}_{\text{{cPD}}}]=c_{g}c_{f}/(c_{f}+k^{2}c_{g})$.

Analogously, the full solution $\pi_{t}$ of the Fokker-Planck equation can be
computed by taking the inverse Fourier transform in $\zeta$ above. We obtain
that $\pi_{t}$ is Gaussian at every time $t>0$ with mean $E(t)^{T}z_{0}$ and
positive definite covariance matrix

$\operatorname{Var}_{Z\sim\pi_{t}}[Z]=E(t)^{T}\Sigma_{0}E(t)+2\begin{pmatrix}\int_{0}^{t}e_{11}^{2}\,\mathrm{d}s&\int_{0}^{t}e_{11}e_{12}\,\mathrm{d}s\\\
\int_{0}^{t}e_{11}e_{12}\,\mathrm{d}s&\int_{0}^{t}e_{12}^{2}\,\mathrm{d}s\end{pmatrix}.$

In the limit $t\to\infty$, we obtain the stationary distribution which has
mean zero and covariance matrix

$\operatorname{Var}_{Z\sim\pi^{\text{{stat}}}_{\text{{cPD}}}}[Z]=\frac{1}{(c_{f}+k^{2}c_{g})(1+\lambda
c_{f}c_{g})}\begin{pmatrix}c_{g}(c_{f}+k^{2}c_{g})+\lambda
c_{f}^{2}c_{g}^{2}&k\lambda c_{f}c_{g}^{2}\\\ k\lambda
c_{f}c_{g}^{2}&k^{2}\lambda c_{g}^{2}\end{pmatrix}$

In the limit $\lambda\to\infty$, the stationary distribution approaches a
degenerate Gaussian (in the sense that the covariance matrix is singular),
with mean zero and covariance matrix

$\frac{1}{(c_{f}+k^{2}c_{g})}\begin{pmatrix}c_{f}c_{g}&kc_{g}\\\
kc_{g}&k^{2}c_{f}^{-1}c_{g}\end{pmatrix}.$

This is identical to $(id,f^{\prime}(k\,\cdot))_{\\#}\mu^{\ast}$, proving
(ii). ∎

The negative result Theorem 3.5 naturally raises the question whether the
diffusion equation can be modified such that the target becomes its stationary
distribution. While we will consider the more general case of non-homogeneous,
problem-dependent diffusion coefficients in Section 3.4, the (algorithmically)
most straightforward approach would be to add homogeneous diffusion not only
in the primal variable $X_{k}$ of the algorithm, but also in the dual $Y_{k}$.
Unfortunately, the next theorem shows this approach can not be successful.

######  Theorem 3.6.

Consider the same target $\mu^{\ast}$ as in Theorem 3.5, but a more general
diffusion process

$\mathrm{d}Z_{t}=-A_{\text{{cPD}}}Z_{t}\mathrm{d}t+B_{\text{{cPD}}}\mathrm{d}W_{t},\quad\text{with
}A_{\text{{cPD}}}=\begin{pmatrix}c_{g}^{-1}&k\\\ -\lambda k&\lambda
c_{f}\end{pmatrix},B_{\text{{cPD}}}=\begin{pmatrix}b_{11}&b_{12}\\\
b_{21}&b_{22}\end{pmatrix}$ (3.8)

where $W_{t}$ is Brownian motion in $\mathbb{R}^{2}$. Then there exists no
choice for the entries of $B_{\text{{cPD}}}\in\mathbb{R}^{2\times 2}$
independent of the constants $c_{g},c_{f}$ and $k$ such that the stationary
solution $\pi^{\text{{stat}}}_{\text{{cPD}}}$ of (3.8) has $x$-marginal
$\mu^{\ast}$.

###### Proof.

We start by noting that the contraction result Theorem 3.4 can easily be
generalized to the more general
$\Sigma_{\text{{cPD}}}=\frac{1}{2}B_{\text{{cPD}}}B_{\text{{cPD}}}^{T}$ with
any constant coefficient $B_{\text{{cPD}}}$, since the relevant Lemma 3.3
holds for any constant diffusion coefficient. In view of this contraction
property, we can assume the initialization $Z_{0}\sim\pi_{0}$ to be Gaussian
with mean $z_{0}$, covariance $\Sigma_{0}$ and obtain
$\pi^{\text{{stat}}}_{\text{{cPD}}}$ by computing $\lim_{t\to\infty}\pi_{t}$.
The reason for initializing at a Gaussian is that we can again obtain the
closed-form solution $\pi_{t}$ of the Cauchy problem for the Fokker-Planck
equation. We carry out the same computations as in the proof of Theorem 3.5,
showing that $\pi_{t}$ is Gaussian with mean $E(t)^{T}z_{0}$ and covariance
matrix

$E(t)^{T}\Sigma_{0}E(t)+2\int_{0}^{t}E(s)^{T}B_{\text{{cPD}}}B_{\text{{cPD}}}^{T}E(s)\,\mathrm{d}s.$

The entries of $E$ can be computed explicitly (we omit the terms here) and by
integrating and taking the limit $t\to\infty$, we find that the $x$-marginal
of the stationary distribution $\pi^{\text{{stat}}}_{\text{{cPD}}}$ has
variance

$\operatorname{Var}[\mu^{\text{{stat}}}_{\text{{cPD}}}]=\frac{(b_{11}^{2}+b_{12}^{2})(\lambda
c_{f}^{2}c_{g}^{2}+c_{g}(c_{f}+k^{2}c_{g}))-2(b_{11}b_{21}+b_{12}b_{22})kc_{f}c_{g}^{2}+(b_{21}^{2}+b_{22}^{2})\lambda^{-1}k^{2}c_{g}^{2}}{2(1+\lambda
c_{f}c_{g})(c_{f}+k^{2}c_{g})}.$

We want to show that there is no choice of $b_{ij}$ (independent of
$c_{f},c_{g}$ and $k$) such that the target $\mu^{\ast}$ coincides with the
$x$-marginal $\mu^{\text{{stat}}}_{\text{{cPD}}}$ of the stationary solution.
The variance of both distributions is equal if and only if

$2c_{g}c_{f}(1+\lambda c_{g}c_{f})=(b_{11}^{2}+b_{12}^{2})(\lambda
c_{f}^{2}c_{g}^{2}+c_{g}(c_{f}+k^{2}c_{g}))-2(b_{11}b_{21}+b_{12}b_{22})kc_{f}c_{g}^{2}+(b_{21}^{2}+b_{22}^{2})\lambda^{-1}k^{2}c_{g}^{2}.$

Isolating the terms including $k^{2}$ shows that this implies

$0=(b_{11}^{2}+b_{12}^{2})c_{g}^{2}+(b_{21}^{2}+b_{22}^{2})\lambda^{-1}c_{g}^{2},$

which is only satisfied if all $b_{ij}$ are zero since $\lambda>0$. Hence, a
problem-independent choice for $B_{\text{{cPD}}}$ that results in the target
being invariant with respect to the diffusion does not exist. ∎

We want to elaborate a little more on the concentration properties of the
stationary solution $\pi_{\text{{cPD}}}^{\text{{stat}}}$. Theorem 3.5 showed
that in the limiting case $\lambda\to\infty$, the stationary distribution does
not have a Lebesgue density anymore since the mass concentrates on the
optimality condition $\partial f^{\ast}(y)=Kx$. In that case, the primal
marginal actually converges to the target distribution, which is intuitively
clear since plugging in $Y_{t}\in\partial f(KX_{t})$ into the equation for
$X_{t}$ reduces the equation to overdamped Langevin diffusion. We will now
generalize these properties of the toy example in Theorem 3.5 in two ways.
Specifically, we will show in Section 3.3 under which conditions the
stationary solution converges in the limit $\lambda\to\infty$ to a singular
push-forward of the target $\mu^{\ast}$ to the space of joint distributions in
$(x,y)$. Also, we now show that the concentration at the optimality generally
does not occur for finite $\lambda$, even though the diffusion coefficient is
degenerate. In the case of smooth coefficients, the Fokker-Planck operator can
be shown to be hypoelliptic, implying that any solution to the dynamics has a
smooth density with respect to the Lebesgue measure for all times $t>0$, even
in the case of a concentrated initialization.

######  Definition 3.7.

Let $\mathcal{L}$ be a partial differential operator with smooth coefficients.
If, for every open set $\Omega$ and every distribution $u$ on $\Omega$,
$\mathcal{L}u\in C^{\infty}(\Omega)$ implies $u\in C^{\infty}(\Omega)$, then
$\mathcal{L}$ is called hypoelliptic.

Every elliptic operator is hypoelliptic, and equally any operator of the form
$-\partial_{t}+\mathcal{L}$ where $\mathcal{L}$ is the generator of a
diffusion process with uniformly positive definite diffusion coefficient. A
special case is overdamped Langevin diffusion with smooth drift, which
naturally has a smooth transition density.

The degenerate diffusion coefficient in (2.5) results in a constant but only
positive semidefinite diffusion coefficient $\Sigma_{\text{{cPD}}}$ in the
Fokker-Planck equation (3.7). In order to prove that the Fokker-Planck
equation is hypoelliptic, we hence have to employ a weaker condition for
hypoellipticity, usually termed ‘Hörmander’s condition’ [34]. In order to
state it, we need the definition of the commutator operation and Lie algebras
spanned by smooth vector fields.

######  Definition 3.8.

For two vector fields $A,B$ on $\mathbb{R}^{n}$, the Lie bracket or commutator
$[A,B]$ is the vector field on $\mathbb{R}^{n}$ defined by

$[A,B](x)=DB(x)A(x)-DA(x)B(x),$

with the derivative matrices $(DA)_{i,j}(x)=\partial_{j}A_{i}(x)$.

######  Lemma 3.9 ([34]).

Assume that the generator $\mathcal{L}$ of a diffusion process with state
space $\mathbb{R}^{n}$ has $C^{\infty}$-smooth coefficients and can be written
in the form

$\mathcal{L}=A_{0}+\sum_{j=1}^{n}A_{j}^{2}$

with first-order differential operators $A_{0},\dots,A_{n}$. Interpreting the
differential operators as vector fields, define a collection of vector fields
$\mathcal{A}_{k}$ by

$\mathcal{A}_{0}\coloneqq\\{A_{j}\,:\,j=1,\dots,n\\},\quad\mathcal{A}_{k+1}=\mathcal{A}_{k}\cup\\{[A,A_{j}]\,:\,A\in\mathcal{A}_{k},j=0,\dots,n\\}.$

If the vector spaces
$\mathcal{A}_{k}(x)\coloneqq\operatorname{span}\left(\\{A(x)\,:\,A\in\mathcal{A}_{k}\\}\right)$
satisfy $\bigcup_{k\geq 0}\mathcal{A}_{k}(x)=\mathbb{R}^{n}$ for every
$x\in\mathbb{R}^{n}$, then $\mathcal{L}$ is said to satisfy the parabolic
Hörmander condition. In that case, $\partial_{t}-\mathcal{L}$ and
$\partial_{t}-\mathcal{L}^{\ast}$ are both hypoelliptic, hence any process
generated by $\cal L$ has $C^{\infty}$-smooth transition probability.In
particular, any solution to the Fokker-Planck equation
$\partial_{t}p=\mathcal{L}^{\ast}p$ is
$C^{\infty}((0,\infty)\times\mathbb{R}^{n})$-smooth.

We refer to Chapter 6 of the book [27] for a discussion in the context of
diffusion SDEs and an application to the Langevin diffusion SDE with friction.
For the diffusion equation arising in our primal-dual setting, the generator
is given by

$\mathcal{L}_{\text{{cPD}}}u=-A_{\text{{cPD}}}^{T}(z)\nabla_{z}u+\Delta_{x}u.$

The following result gives a sufficient condition under which the primal-dual
dynamics considered here satisfy Hörmander’s condition, so that the
corresponding Fokker-Planck operator is hypoelliptic.

######  Theorem 3.10.

Assume that $g\in C^{\infty}(\mathcal{X})$ and $f^{\ast}\in
C^{\infty}(\mathcal{Y})$, so that the diffusion (2.6) has smooth coefficients.
If $\dim(\mathcal{Y})=m\leq d=\dim(\mathcal{X})$ and $K$ has full rank $m$,
then $\mathcal{L}_{\text{{cPD}}}$ satisfies Hörmander’s condition. As a
consequence, the stationary solution $\pi^{\text{{stat}}}_{\text{{cPD}}}$ has
a $C^{\infty}(\mathcal{Z})$-smooth Lebesgue density.

###### Proof.

We apply Lemma 3.9. $\mathcal{L}_{\text{{cPD}}}$ has the correct form with
vector fields $A_{j}=\partial_{x_{j}}$, $j=1,\dots,d$ and

$A_{0}=-A_{\text{{cPD}}}^{T}(z)\nabla_{z}=\sum_{j=1}^{d}\left(-\frac{\partial
g}{\partial
x_{j}}(x)-k_{j}^{T}y\right)\partial_{x_{j}}+\lambda\sum_{i=1}^{m}\left((k^{i})^{T}x-\frac{\partial
f^{\ast}}{\partial y_{i}}(y)\right)\partial_{y_{i}},$

where we denoted the columns of $K$ by $k_{1},\dots,k_{d}$ and the rows by
$k^{1},\dots,k^{m}$. Clearly, the vector fields $A_{1},\dots,A_{k}$ span the
first $d$ dimensions $\partial_{x_{j}}$, $j=1,\dots,d$, i.e.
$\mathcal{A}_{0}(z)=\mathcal{X}$ for all $z\in\mathcal{Z}$. The Lie brackets
in $\mathcal{A}_{1}$ are given by

$[A_{j},A_{0}]=-(H_{g}(x))_{j}^{T}\nabla_{x}-\lambda k_{j}^{T}\nabla_{y},$

where $H_{g}$ is the Hessian of $g$. Due to $\lambda>0$ and the full rank
assumption, the columns $k_{j}$ of $K$ span $\mathcal{Y}$ and we have
$\mathcal{A}_{1}(z)=\mathcal{Z}$ for all $z\in\mathcal{Z}$. ∎

###  3.3 Convergence to Overdamped Langevin Diffusion

We now continue to show that in the limiting case $\lambda\to\infty$ and under
some smoothness assumption on $f$, the stationary solution
$\pi^{\text{{stat}}}_{\text{{cPD}}}$ actually becomes singular, and further
that the primal marginal distribution coincides in the limit with the target.
Note that the formal case $\lambda=\infty$ was the subject of a recent paper
[3], where the authors showed convergence in the case of perfect concentration
$y\in\partial f(Kx)$.

Consider the rescaled model $(\lambda=\frac{1}{\varepsilon})$ taking the form

$\left\\{\begin{array}[]{l}dX_{t}\in-(\partial
g(X_{t})+K^{*}Y_{t})dt+\sqrt{2}dW_{t}\\\ \varepsilon dY_{t}\in-(\partial
f^{*}(Y_{t})-KX_{t})dt,\end{array}\right.$ (3.9)

where we again denote the distribution of $(X_{t},Y_{t})$ by $\pi_{t}$, with
$x$-marginal $\mu_{t}$. We are now interested in the behaviour of $\mu_{t}$ in
the limit $\varepsilon\rightarrow 0$. Intuitively, as $\varepsilon$ becomes
small, in this formulation it becomes clear that $Y_{t}$ is forced to
concentrate according to the optimality condition $Y_{t}\in\partial
f(KX_{t})$. In this concentrated case, plugging in the dual variable into the
primal inclusion shows that the dynamics actually coincide with overdamped
Langevin dynamics, which has the correct invariant distribution. Time-
discretizations of the limiting case $\varepsilon=0$ have been considered in
[3], where an intermediate iterate $Y_{t}\in\partial f(KX_{t})$ is computed
explicitly. Indeed, we can use a coupling argument between (3.9) and
overdamped Langevin dynamics to bound the Wasserstein distance between the
stationary state of (3.9) and the target.

######  Theorem 3.11.

Assume that $f\in C^{3}$, that $g,f^{\ast}$ are $\omega_{g}$\- respectively
$\omega_{f^{\ast}}$-strongly convex, respectively, and that the target
distribution $\mu^{\ast}$ satisfies

$\displaystyle\mathbb{E}_{X\sim\mu^{\ast}}\left[\lVert\nabla(\nabla f\circ
K)(X)\rVert_{F}^{2}\right]=:C_{1}<\infty,$
$\displaystyle\mathbb{E}_{X\sim\mu^{\ast}}\left[\lVert\nabla\cdot\nabla(\nabla
f\circ K)(X)\rVert^{2}\right]=:C_{2}<\infty,$
$\displaystyle\mathbb{E}_{X\sim\mu^{\ast}}\left[\lVert\nabla(\nabla f\circ
K)(X)\left(\nabla g(X)+K^{T}\nabla
f(KX)\right)\rVert^{2}\right]=:C_{3}<\infty,$

where in the second line the Laplacian is understood acting component-wise on
the $\mathcal{Y}$-valued function $\nabla f\circ K$. For
$0<\varepsilon<\omega_{f^{\ast}}/\omega_{g}$ and a constant
$C(\varepsilon)=\sqrt{\varepsilon
C_{1}/\omega_{g}}+\mathcal{O}(\varepsilon^{3/2})$, the following two
statements hold:

  1. 1.

If the primal-dual dynamics (3.9) are initialized with $X_{0}\sim\mu_{0}$ and
$Y_{0}$ concentrated at the optimality condition $Y_{0}=\nabla f(KX_{0})$,
then for every $t>0$ it holds

$\mathcal{W}_{2}(\mu_{t},\mu^{\ast})\leq\left((1+\varepsilon\omega_{f^{\ast}}\lVert
K\rVert)^{1/2}\mathcal{W}_{2}(\mu_{0},\mu^{\ast})+C(\varepsilon)\right)e^{-\omega_{g}t}+C(\varepsilon).$

  2. 2.

The $x$-marginal $\mu^{\text{{stat}}}_{\text{{cPD}}}$ of the stationary
solution $\pi^{\text{{stat}}}_{\text{{cPD}}}$ of (3.9) satisfies

$\mathcal{W}_{2}(\mu^{\text{{stat}}}_{\text{{cPD}}},\mu^{\ast})\leq
C(\varepsilon).$

###### Proof.

Assume $(X_{t},Y_{t})$ with distribution $\pi_{t}$ solves (3.9). We construct
a coupling with a process $\tilde{X}_{t}$ which is initialized as
$\tilde{X}_{0}\sim\mu^{\ast}$ and solves the usual overdamped Langevin
dynamics

$d\tilde{X}_{t}\in-(\partial g(\tilde{X}_{t})+K^{T}\nabla
f(K\tilde{X}_{t}))dt+\sqrt{2}dW_{t}$ (3.10)

with the same Brownian motion $W_{t}$ as in (3.9). Since $\mu^{\ast}$ is the
invariant distribution of (3.10), it follows $\tilde{X}_{t}\sim\mu^{\ast}$ for
all $t>0$. The joint measure $\Pi_{t}$ of $(X_{t},Y_{t},\tilde{X}_{t})$ solves
the Fokker-Planck equation

$\displaystyle\frac{\partial}{\partial t}\Pi_{t}$
$\displaystyle=\nabla_{x}\cdot((p+K^{*}y)\Pi_{t})+\frac{1}{\varepsilon}\nabla_{y}\cdot((q-Kx)\Pi_{t})+\nabla_{\tilde{x}}\cdot((\tilde{p}+K^{T}\tilde{y})\Pi_{t})$
(3.11)
$\displaystyle\qquad+(\nabla_{x}+\nabla_{\tilde{x}})\cdot(\nabla_{x}+\nabla_{\tilde{x}})\Pi_{t}.$

where as before $p\in\partial g(x),q\in\partial
f^{\ast}(y),\tilde{p}\in\partial g(\tilde{x})$. With the artificial,
concentrated dual variable $\tilde{Y}_{t}=\nabla f(\tilde{X}_{t})$, we denote
the singular joint law of $(\tilde{X}_{t},\tilde{Y}_{t})$ by
$\pi^{\ast}=(I,\nabla f\circ K)_{\\#}\mu^{\ast}$.

We now consider the transport distance $\mathcal{T}_{1,\varepsilon}$ on
$\mathcal{P}(\mathcal{Z})$ in order to bound the distance between $\pi_{t}$
and $\pi^{\ast}$, which due to the concentration of $\pi^{\ast}$ can be
written as

$\mathcal{T}_{1,\varepsilon}(\pi_{t},\pi^{\ast})=\inf_{\Pi\in\Gamma(\pi,\mu^{\ast})}\left(\mathbb{E}_{(x,y,\tilde{x})\sim\Pi}\left[d_{1,\varepsilon}(x,y,\tilde{x})^{2}\right]\right)^{1/2},$

with an augmented distance function

$d_{1,\varepsilon}(x,y,\tilde{x}):=\left(\lVert
x-\tilde{x}\rVert^{2}+\varepsilon\lVert y-\nabla
f(K\tilde{x})\rVert^{2}\right)^{1/2}.$

By using the joint Fokker-Planck equation (3.11) and the notations
$\tilde{y}(\tilde{x})\coloneqq\nabla f(K\tilde{x})$ as well as
$M\coloneqq\nabla\tilde{y}$ and $m\coloneqq\nabla\cdot\nabla\tilde{y}$, we
obtain

$\displaystyle\frac{1}{2}\frac{\,\mathrm{d}}{\,\mathrm{d}t}\mathbb{E}_{(x,y,\tilde{x})\sim\Pi_{t}}\left[d_{1,\varepsilon}(x,y,\tilde{x})^{2}\right]$
$\displaystyle=\frac{1}{2}\int
d_{1,\varepsilon}(x,y,\tilde{x})^{2}\left(\nabla_{x}\cdot((p+K^{T}y)\Pi_{t})+\frac{1}{\varepsilon}\nabla_{y}\cdot((q-Kx)\Pi_{t})+\nabla_{\tilde{x}}\cdot((\tilde{p}+K^{T}\tilde{y}(\tilde{x}))\Pi_{t})\right)\,\mathrm{d}(x,y,\tilde{x})$
$\displaystyle\qquad+\frac{1}{2}\int
d_{1,\varepsilon}(x,y,\tilde{x})^{2}\left((\nabla_{x}+\nabla_{\tilde{x}})\cdot(\nabla_{x}+\nabla_{\tilde{x}})\Pi_{t}\right)\,\mathrm{d}(x,y,\tilde{x})$
$\displaystyle=-\int\left((x-\tilde{x})\cdot(p+K^{T}y)+(y-\tilde{y}(\tilde{x}))\cdot(q-Kx)+(\tilde{x}-x)\cdot(\tilde{p}+K^{T}\tilde{y}(\tilde{x})\right)\,\mathrm{d}\Pi_{t}$
$\displaystyle\qquad-\varepsilon\int\left(M(\tilde{x})(\tilde{y}(\tilde{x})-y)\right)\cdot(\tilde{p}+K^{T}\tilde{y}(\tilde{x}))\,\mathrm{d}\Pi_{t}-\varepsilon\int\left(M(\tilde{x})(\tilde{y}(\tilde{x})-y)\right)\cdot(\nabla_{x}+\nabla_{\tilde{x}})\Pi_{t}d(x,y,\tilde{x})$
$\displaystyle=-\int(x-\tilde{x})\cdot(p-\tilde{p})\,\mathrm{d}\Pi_{t}-\int(y-\tilde{y}(\tilde{x}))\cdot(q-K\tilde{x})\,\mathrm{d}\Pi_{t}$
$\displaystyle\qquad+\varepsilon\int\left[(\tilde{y}(\tilde{x})-y)\cdot\left(m(\tilde{x})-M^{T}(\tilde{x})(\tilde{p}+K^{T}\tilde{y}(\tilde{x}))\right)+\lVert
M(\tilde{x})\rVert_{F}^{2}\right]\,\mathrm{d}\Pi_{t}$

For the terms in the second to last line, we can use strong convexity of $g$
and $f^{\ast}$ to estimate

$\displaystyle\int(x-\tilde{x})\cdot(p-\tilde{p})\,\mathrm{d}\Pi_{t}\geq\omega_{g}\mathbb{E}_{(x,y,\tilde{x})\sim\Pi_{t}}\left[\lVert
x-\tilde{x}\rVert^{2}\right]$
$\displaystyle\int(y-\tilde{y}(\tilde{x}))\cdot(q-K\tilde{x})\,\mathrm{d}\Pi_{t}\geq\frac{\omega_{f^{\ast}}}{\varepsilon}\mathbb{E}_{(x,y,\tilde{x})\sim\Pi_{t}}\left[\lVert
y-\tilde{y}(\tilde{x})\rVert_{\varepsilon}^{2}\right]$

Since the $\tilde{x}$-marginal of $\Pi_{t}$ is the target $\mu^{\ast}$, by
assumption we have

$\varepsilon\int\lVert
M(\tilde{x})\rVert_{F}^{2}\,\mathrm{d}\Pi_{t}\leq\varepsilon C_{1}.$

In a similar way, using the Cauchy-Schwarz inequality and Young’s inequality,
for any $\delta>0$ we can bound the final term by

$\displaystyle\varepsilon\int(\tilde{y}(\tilde{x})-y)$
$\displaystyle\cdot\left(m(\tilde{x})-M^{T}(\tilde{x})(\tilde{p}+K^{T}\tilde{y}(\tilde{x}))\right)\,\mathrm{d}\Pi_{t}$
$\displaystyle\leq\varepsilon\int\lVert\tilde{y}(\tilde{x})-y\rVert\cdot\lVert
m(\tilde{x})-M^{T}(\tilde{x})(\tilde{p}+K^{T}\tilde{y}(\tilde{x}))\rVert\,\mathrm{d}\Pi_{t}$
$\displaystyle\leq\frac{\delta}{2}\mathbb{E}_{(x,y,\tilde{x})\sim\Pi_{t}}\left[\lVert
y-\tilde{y}(\tilde{x})\rVert_{\varepsilon}^{2}\right]+\frac{\varepsilon}{2\delta}\mathbb{E}_{\tilde{x}\sim\mu^{\ast}}\left[\lVert
m(\tilde{x})-M^{T}(\tilde{x})(\tilde{p}+K^{T}\tilde{y}(\tilde{x}))\rVert^{2}\right]$
$\displaystyle\leq\frac{\delta}{2}\mathbb{E}_{(x,y,\tilde{x})\sim\Pi_{t}}\left[\lVert
y-\tilde{y}(\tilde{x})\rVert_{\varepsilon}^{2}\right]+\frac{\varepsilon(C_{2}+C_{3})}{2\delta}.$

Denoting
$\varphi(t):=\mathbb{E}_{(x,y,\tilde{x})\sim\Pi_{t}}\left[d_{1,\varepsilon}(x,y,\tilde{x})^{2}\right]$,
we therefore have

$\varphi^{\prime}(t)\leq-2\omega\varphi(t)+2\varepsilon
C_{1}+\frac{\varepsilon(C_{2}+C_{3})}{\delta},$

where we abbreviate
$\omega\coloneqq\min\\{\omega_{g},\frac{\omega_{f^{\ast}}}{\varepsilon}-\frac{\delta}{2}\\}$.
Using Gronwall’s lemma, this implies

$\displaystyle\varphi(t)$ $\displaystyle\leq e^{-2\omega
t}\varphi(0)+\int_{0}^{t}e^{-2\omega(t-s)}\left(2\varepsilon
C_{1}+\frac{\varepsilon(C_{2}+C_{3})}{\delta}\right)\,\mathrm{d}s$
$\displaystyle=e^{-2\omega t}\varphi(0)+\left(2\varepsilon
C_{1}+\frac{\varepsilon(C_{2}+C_{3})}{\delta}\right)\frac{1-e^{-2\omega
t}}{2\omega}$

By the definition of $\mathcal{T}_{1,\varepsilon}$ and taking the infimum over
all possible initializations $\Pi_{0}$ with $(x,y)$-marginal $\pi_{0}$ and
$(\tilde{x},\tilde{y})$-marginal $\pi^{\ast}=(id,\nabla f\circ
K)_{\\#}\mu^{\ast}$ on the right side, we get

$\mathcal{T}_{1,\varepsilon}^{2}(\pi_{t},\pi^{\ast})\leq\varphi(t)\leq
e^{-2\omega
t}\mathcal{T}_{1,\varepsilon}^{2}(\pi_{0},\pi^{\ast})+\left(2\varepsilon
C_{1}+\frac{\varepsilon(C_{2}+C_{3})}{\delta}\right)\frac{1-e^{-2\omega
t}}{2\omega}.$ (3.12)

For $\varepsilon<\omega_{f^{\ast}}/\omega_{g}$, the bound on the right is
minimized by $\delta=2(\frac{\omega_{f^{\ast}}}{\varepsilon}-\omega_{g})$ for
any finite $t$ as well as in the limit $t\to\infty$. This choice implies
$\omega=\omega_{g}$ and we denote the resulting constant terms by

$C(\varepsilon)\coloneqq\left(\frac{\varepsilon
C_{1}}{\omega_{g}}+\frac{\varepsilon(C_{2}+C_{3})}{4\omega_{g}(\omega_{f^{\ast}}-\varepsilon\omega_{g})}\right)^{1/2}=\sqrt{\varepsilon
C_{1}/\omega_{g}}+\mathcal{O}(\varepsilon^{3/2}).$

The first desired inequality for $t>0$ follows by taking the square root of
(3.12) and using that
$\mathcal{W}_{2}(\mu_{t},\mu^{\ast})\leq\mathcal{T}_{1,\varepsilon}(\pi_{t},\pi^{\ast})$
and, assuming a concentrated initialization $Y_{0}=\nabla f(KX_{0})$, that

$\mathcal{T}_{1,\varepsilon}(\pi_{0},\pi^{\ast})\leq(1+\varepsilon\omega_{f^{\ast}}\lVert
K\rVert)^{1/2}\mathcal{W}_{2}(\mu_{0},\mu^{\ast}).$

Since $\pi_{t}$ converges to the stationary distribution of interest
$\pi^{\text{{stat}}}_{\text{{cPD}}}$, taking the limit $t\to\infty$ in (3.12)
and using
$\mathcal{W}_{2}(\mu^{\text{{stat}}}_{\text{{cPD}}},\mu^{\ast})\leq\mathcal{T}_{1,\varepsilon}(\pi^{\text{{stat}}}_{\text{{cPD}}},\pi^{\ast})$
proves the second statement. ∎

###  3.4 A Bias-Corrected Modification of Primal-Dual Diffusion

In this section, we want to follow up on the negative result of Theorem 3.6,
where we saw that there is no way to construct a diffusion SDE with the
primal-dual type drift coefficient and homogeneous diffusion that has the
correct stationary distribution. Naturally, we can also consider modifications
of the diffusion SDE with problem-dependent and/or non-homogeneous diffusion
coefficient, so that the samples become (in continuous time) unbiased.

We aim to correct the Fokker-Planck equation (3.7) by adding terms such that
$\pi^{\ast}=(I,\nabla f\circ K)_{\\#}\mu^{\ast}$ becomes the stationary
solution of the PDE. Throughout this subsection, we will assume that $g\in
C^{1}(\mathbb{R}^{d})$ and $f\in C^{3}(\mathbb{R}^{m})$, we will denote as
before $M:\mathbb{R}^{d}\to\mathbb{R}^{d\times m}$ for
$M(x)\coloneqq\nabla(\nabla f\circ K)^{T}(x)=K^{T}H_{f}(Kx)$. Instead of
(2.5), consider the modified stochastic differential equation

$\,\mathrm{d}Z_{t}=-A_{\text{{mod}}}(Z_{t})\,\mathrm{d}t+B_{\text{{mod}}}(Z_{t})\,\mathrm{d}W_{t},$
(3.13)

where the drift and diffusion coefficients are given by

$\displaystyle A_{\text{{mod}}}(z)$
$\displaystyle\coloneqq\begin{pmatrix}\nabla g(x)+K^{T}y\\\ \lambda(\nabla
f^{*}(y)-Kx)+M(x)^{T}(\nabla g(x)+K^{T}\nabla f(Kx))\end{pmatrix}$
$\displaystyle B_{\text{{mod}}}(z)$
$\displaystyle\coloneqq\sqrt{2}\,\nabla(id,\nabla f\circ
K)(x)=\sqrt{2}\begin{pmatrix}I\\\ M(x)^{T}\end{pmatrix}.$

Note that $B_{\text{{mod}}}$ differs from $B_{\text{{cPD}}}$ in a non-
homogeneous coefficient in the dual part $y$. The drift $A_{\text{{mod}}}$
contains the same terms as $A_{\text{{cPD}}}$, and an additional term
$-M(x)^{T}\nabla\log\mu^{\ast}(x)$.

######  Remark 3.12.

As an intuitive explanation why (3.13) has the desired stationary
distribution, note that we want the stationary solution to concentrate its
mass in $\mathcal{Z}=\mathbb{R}^{d+m}$ on the $d$-dimensional manifold
$\mathcal{U}$ described by $y=\nabla f(Kx)$. The tangent space to
$\mathcal{U}$ at any point $(x,\nabla f(Kx))$, denoted in the following by
$\mathcal{T}_{x}\mathcal{U}$, is spanned by the columns of $\nabla(id,\nabla
f\circ K)=(I,M(x))^{T}$. This coincides with the column space of the modified
diffusion coefficient
$\Sigma_{\text{{mod}}}\coloneqq\frac{1}{2}B_{\text{{mod}}}B_{\text{{mod}}}^{T}$,
hence locally on $\mathcal{U}$, the diffusion is confined to the manifold. The
term $\nabla f^{\ast}(y)-Kx$ in the drift vanishes upon concentration of the
mass on $\mathcal{U}$, the remaining terms are just the gradient of the
negative log-density $-\nabla\log\mu^{\ast}(x)=\nabla g(x)+K^{T}\nabla f(Kx)$
in the primal variable $x$ mapped to the tangent space of $\mathcal{U}$ by the
map $(I,M(x))^{T}$. The drift hence does not transport mass away from
$\mathcal{U}$ and locally on $\mathcal{U}$, (3.13) coincides with overdamped
Langevin diffusion.

In order to show formally that the partially singular distribution
$\pi^{\ast}$ is a stationary solution of (3.13), we have to consider the weak
formulation of the corresponding Fokker-Planck equation, given by

$\frac{\,\mathrm{d}}{\,\mathrm{d}t}\int\varphi
d\pi_{t}=\int\nabla_{z}^{T}\Sigma_{\text{{mod}}}(z)\nabla_{z}\varphi\,\mathrm{d}\pi_{t}-\int\nabla\varphi\cdot
A_{\text{{mod}}}(z)\,\mathrm{d}\pi_{t}\quad\forall\varphi\in
C^{\infty}(\mathcal{Z}).$ (3.14)

######  Theorem 3.13.

The distribution $\pi^{\ast}\coloneqq(I,\nabla f\circ K)_{\\#}\mu^{*}$ is a
stationary solution of (3.14).

In principle, this theorem might suggest that a time discretization of the
modified (3.13) may be a better idea than the primal-dual sampling algorithm
(2.4). Note, however, apart from the higher-order derivative terms in $f$
(which counter the motivation to use a primal-dual scheme in the first place),
the term $A_{\text{{mod}}}$ contains just the gradient of the target log-
density, i.e. we explicitly add overdamped Langevin diffusion in the dual
variable. This suggests that algorithmically, the dualization is somewhat
redundant and any sampler would be inferior to a standard unadjusted Langevin
algorithm. In order to prove Theorem 3.13, we need the following vector
calculus identity.

######  Lemma 3.14.

Let $\psi:\mathbb{R}^{d}\times\mathbb{R}^{m}\rightarrow\mathbb{R}^{n}$,
$\eta:\mathbb{R}^{d}\rightarrow\mathbb{R}^{m}$ and
$S:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times n}$, each with $C^{1}$
entries. Define $r:\mathbb{R}^{d}\to\mathbb{R}^{n}$ by
$r_{j}(x)\coloneqq\sum_{i}\frac{\partial S_{ij}}{\partial x_{i}}(x)$ and
denote the partial Jacobians of $\psi$ by
$D_{x}\psi:\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{n\times d}$ with
$(D_{x}\psi)_{i,j}(x,y)\coloneqq\frac{\partial\psi_{i}}{\partial x_{j}}(x,y)$
and $D_{y}\psi:\mathbb{R}^{d}\times\mathbb{R}^{m}\to\mathbb{R}^{n\times m}$
where $(D_{y}\psi)_{i,k}(x,y)\coloneqq\frac{\partial\psi_{i}}{\partial
y_{k}}(x,y)$. Then for all $x\in\mathbb{R}^{d}$ it holds

$\nabla_{x}\cdot(S(x)\psi(x,\eta(x)))=r(x)\cdot\psi(x,\eta(x))+\operatorname{tr}(S(x)D_{x}\psi(x,\eta(x))+\operatorname{tr}(S(x)D_{y}\psi(x,\eta(x))D_{x}\eta(x)).$
(3.15)

###### Proof.

The proof is a straightforward application of product and chain rule.

$\displaystyle\nabla_{x}\cdot($ $\displaystyle S(x)\psi(x,\eta(x)))$
$\displaystyle=\sum_{i=1}^{d}\frac{\partial}{\partial
x_{i}}(\sum_{j=1}^{d}S_{ij}(x)\psi_{j}(x,\eta(x)))$
$\displaystyle=\sum_{ij}\frac{\partial S_{ij}}{\partial
x_{i}}(x)\psi_{j}(x,\eta(x))+S_{ij}(x)\frac{\partial\psi_{j}}{\partial
x_{i}}(x,\eta(x))+S_{ij}(x)\sum_{k=1}^{m}\frac{\partial\psi_{j}}{\partial
y_{k}}(x,\eta(x))\frac{\partial h_{k}}{\partial x_{i}}(x)$
$\displaystyle=r(x)\cdot\psi(x,\eta(x))+\operatorname{tr}(S(x)D_{x}\psi(x,\eta(x))+\operatorname{tr}(S(x)D_{y}\psi(x,\eta(x))D_{x}\eta(x)).$

∎

###### Proof of Theorem 3.13.

It suffices to check that the right side of (3.14) evaluates to zero for any
$\varphi\in C^{\infty}(\mathbb{R}^{d}\times\mathbb{R}^{d})$ when we insert
$\pi^{\ast}$:

$\displaystyle\int\begin{pmatrix}\nabla_{x}\\\
\nabla_{y}\end{pmatrix}^{T}\begin{pmatrix}I&M(x)\\\
M(x)^{T}&M(x)M(x)^{T}\end{pmatrix}\begin{pmatrix}\nabla_{x}\\\
\nabla_{y}\end{pmatrix}\varphi\,\mathrm{d}\pi^{\ast}-\int\nabla_{x}\varphi\cdot(\nabla
g(x)+K^{*}y)\,\mathrm{d}\pi^{\ast}$ (3.16a)
$\displaystyle-\int\nabla_{y}\varphi\cdot(\lambda\nabla f^{*}(y)-\lambda
Kx)\,\mathrm{d}\pi^{\ast}-\int\nabla_{y}\varphi\cdot M(x)^{T}(\nabla
g(x)+K^{T}\nabla f(Kx))\,\mathrm{d}\pi^{\ast}$ (3.16b)

Using $(\nabla f)^{-1}=\nabla f^{\ast}$ and the fact that $\pi^{\ast}$ is
partially concentrated on $y=\nabla f(Kx)$, it follows that

$\int\nabla_{y}\varphi\cdot(\lambda\nabla f^{\ast}(y)-\lambda
Kx)\,\mathrm{d}\pi^{\ast}=0.$

Using partial integration and Lemma 3.14 with $S\equiv I$,
$\psi=\nabla_{x}\varphi$ and $\eta=\nabla f\circ K$, the second integral in
(3.16a) can be rewritten as

$\displaystyle-\int\nabla_{x}$ $\displaystyle\varphi\cdot(\nabla
g(x)+K^{*}y)\,\mathrm{d}\pi^{\ast}$
$\displaystyle=\int\nabla_{x}\varphi(x,\nabla f(Kx))\cdot(-\nabla
g(x)-K^{*}\nabla f(Kx))\exp(-f(Kx)-g(x))\,\mathrm{d}x$
$\displaystyle=\int\nabla_{x}\varphi(x,\nabla
f(Kx))\cdot\nabla_{x}(\exp(-f(Kx)-g(x)))\,\mathrm{d}x$
$\displaystyle=-\int\nabla_{x}\cdot(\nabla_{x}\varphi(x,\nabla
f(Kx)))\,\mathrm{d}\mu^{\ast}(x)$
$\displaystyle=-\int\operatorname{tr}(D_{x}\nabla_{x}\varphi(x,\nabla
f(Kx)))\,\mathrm{d}\mu^{\ast}(x)-\int\operatorname{tr}(D_{y}\nabla_{x}\varphi(x,\nabla
f(Kx))M(x)^{T})\,\mathrm{d}\mu^{\ast}(x)$
$\displaystyle=-\int\Delta_{x}\varphi(x,\nabla
f(Kx))\,\mathrm{d}\mu^{\ast}(x)-\int\nabla_{y}\cdot(M(x)^{T}\nabla_{x}\varphi(x,\nabla
f(Kx)))\,\mathrm{d}\mu^{\ast}(x)$

where we used that $M(x)^{T}=D_{x}\eta(x)$ for $\eta=\nabla f\circ K$. These
two terms cancel with the terms from the first column of the diffusion matrix
in (3.16a). The second integral in (3.16b) can also be reformulated using
partial integration and then Lemma 3.14, this time with $S=M$,
$\psi=\nabla_{y}\varphi$ and $\eta=\nabla f\circ K$. In the notation of Lemma
3.14, we get $r:\mathbb{R}^{d}\to\mathbb{R}^{m}$ with
$r_{j}(x)\coloneqq\sum_{i}\frac{\partial M_{ij}}{\partial x_{i}}(x)$ and
obtain

$\displaystyle-\int$ $\displaystyle\nabla_{y}\varphi(x,y)\cdot M(x)^{T}(\nabla
g(x)+K^{T}\nabla f(Kx))\,\mathrm{d}\pi^{\ast}$
$\displaystyle=\int\nabla_{y}\varphi(x,\nabla f(Kx))\cdot M(x)^{T}(-\nabla
g(x)-K^{T}\nabla f(Kx))\exp(-f(Kx)-g(x))\,\mathrm{d}x$
$\displaystyle=\int\nabla_{y}\varphi(x,\nabla f(Kx))\cdot
M(x)^{T}\nabla_{x}(\exp(-f(Kx)-g(x)))\,\mathrm{d}x$
$\displaystyle=-\int\nabla_{x}\cdot(M(x)\nabla_{y}\varphi(x,\nabla
f(Kx)))\,\mathrm{d}\mu^{\ast}(x)$ $\displaystyle=-\int
r(x)\cdot\nabla_{y}\varphi(x,\nabla
f(Kx))\,\mathrm{d}\mu^{\ast}(x)-\int\operatorname{tr}(M(x)D_{x}\nabla_{y}\varphi(x,\nabla
f(Kx)))\,\mathrm{d}\mu^{\ast}(x)$
$\displaystyle\qquad\qquad-\int\operatorname{tr}(M(x)D_{y}\nabla_{y}\varphi(x,\nabla
f(Kx))M(x)^{T})\,\mathrm{d}\mu^{\ast}(x)$

Since $\operatorname{tr}(MD_{y}\nabla_{y}\varphi
M^{T})=\operatorname{tr}(M^{T}MD_{y}\nabla_{y}\varphi)=\nabla_{y}\cdot(M^{T}M\nabla_{y}\varphi)$,
the term in the last line cancels with the term from the lower right entry in
the diffusion matrix. By the product rule,
$\nabla_{x}\cdot(M\nabla_{y}\varphi)=\operatorname{tr}(MD_{x}\nabla_{y}\varphi)+r\cdot\nabla_{y}\varphi$,
hence the second to last line cancels with the term from the upper right entry
in the diffusion matrix. It follows that (3.16) is equal to zero. ∎

######  Remark 3.15.

In general, it is hard to generalize a contraction result like Lemma 3.3 to
the case of non-homogeneous diffusion coefficients. In section 2 of [13], this
question is addressed in more detail for the heat equation with non-
homogeneous diffusion. The authors show that it becomes possible to prove
stability in a Wasserstein distance with a modified underlying metric, where
the metric depends on the diffusion coefficient via an elliptic PDE.
Since the modified equation (3.13) is of limited numerical interest, we only
sketch some main ideas and omit a detailed proof. At every point
$x\in\mathcal{X}$ one can define an inner product $\langle
u,v\rangle_{x}\coloneqq u^{T}V^{-1}(x)v$ where $V(x)=W(x)W(x)^{T}$ is
symmetric positive definite. The first $d$ columns of $W$ consist of a basis
of the tangent space $\mathcal{T}_{x}\mathcal{U}$ of $\mathcal{U}$ at $x$, and
the remaining $m$ columns span the normal space $\mathcal{N}_{x}\mathcal{U}$.
We then need to equip $\mathcal{Z}$ with the induced Riemannian metric

$d_{\text{{mod}}}(z,\tilde{z})=\inf_{\zeta(0)=z,\,\zeta(1)=\tilde{z}}\int_{0}^{1}\left(\dot{\zeta}(t)^{T}W^{-1}(\zeta(t))\dot{\zeta}(t)\right)^{1/2}\,\mathrm{d}t.$

The difficulty then usually consists in finding a more convenient form of
$d_{\text{{mod}}}(z,\tilde{z})$, which allows to integrate the Fokker-Planck
equation against $d_{\text{{mod}}}^{2}(z,\tilde{z})$ in order to show
contraction. We leave this to future work.

######  Example 3.16.

We can check that in the one-dimensional quadratic case that was considered as
an illustrative counterexample in Theorems 3.5 and 3.6, the modified equation
now has the desired invariant distribution $\pi^{\ast}$. Furthermore, every
solution $\pi_{t}$ of the Fokker-Planck equation converges to $\pi^{\ast}$ as
$t\to\infty$. This can be seen as follows: We carry out the same computations
as in the proof of Theorem 3.5 to obtain the Fourier transform $p_{t}$ of the
solution $\pi_{t}$, given by

$\displaystyle
p_{t}(\zeta)=p_{0}(\exp(-tA_{\text{{mod}}}^{T})\zeta)\exp\left(-\frac{1}{2}\int_{0}^{t}(B_{\text{{mod}}}^{T}\exp(-sA_{\text{{mod}}}^{T})\zeta)^{2}\,\mathrm{d}s\right),$
$\displaystyle A_{\text{{mod}}}=\begin{pmatrix}c_{g}^{-1}&k\\\ -\lambda
k+kc_{f}^{-1}c_{g}^{-1}+k^{3}c_{f}^{-2}&\lambda c_{f}\end{pmatrix},\quad
B_{\text{{mod}}}=\begin{pmatrix}1\\\ kc_{f}^{-1}\end{pmatrix},$

This time, the quadratic form under the integral is not positive definite in
$\zeta$ since $B_{\text{{mod}}}$ is an eigenvector of $A_{\text{{mod}}}$ with
eigenvalue $v_{h}\coloneqq c_{g}^{-1}+k^{2}c_{f}^{-1}$. Instead, it holds

$\int_{0}^{t}(B_{\text{{mod}}}^{T}\exp(-sA_{\text{{mod}}}^{T})\zeta)^{2}\,\mathrm{d}s=v_{h}^{-1}\left(1-e^{-tv_{h}}\right)(B_{\text{{mod}}}^{T}\zeta)^{2}.$

Taking the inverse Fourier transform shows that for every (possibly
degenerate) initialization $\pi_{0}$ the solution $\pi_{t}$ converges to the
invariant distribution, which is degenerated Gaussian with mean zero and
covariance matrix $v_{h}^{-1}B_{\text{{mod}}}B_{\text{{mod}}}^{T}$. This is
precisely the target distribution $\pi^{\ast}$.

##  4 Stability Analysis in Discrete Time

We want to apply similar techniques as in Section 3 to the time discretized
setting of the algorithm (2.4). Indeed, couplings like the ones we employed in
the continuous time setting above have been used before in the stability and
convergence analysis of Langevin algorithms. We begin this section by
recalling such a computation in the simple case of ULA as an instructive
example, and then expand to the primal-dual algorithm considered here.

Consider two Markov chains $(X_{n})_{n},(\tilde{X}_{n})_{n}$ defined by
initial values $X_{0},\tilde{X}_{0}$ and

$\displaystyle\xi^{n+1}$ $\displaystyle\sim\mathrm{N}(0,I)$ $\displaystyle
X^{n+1}$ $\displaystyle=X^{n}-\tau\nabla h(X_{n})+\sqrt{2\tau}\xi^{n+1}$
$\displaystyle\tilde{X}^{n+1}$ $\displaystyle=\tilde{X}^{n}-\tau\nabla
h(\tilde{X}^{n})+\sqrt{2\tau}\xi^{n+1}$

Then we can compute that (see [17], Prop. 3 where this calculation has been
carried out)

$\displaystyle\lVert X^{n+1}-\tilde{X}^{n+1}\rVert^{2}$ $\displaystyle=\lVert
X^{n}-\tau\nabla h(X^{n})-\tilde{X}^{n}+\tau\nabla h(\tilde{X}^{n})\rVert^{2}$
$\displaystyle=\lVert X^{n}-\tilde{X}^{n}\rVert^{2}-2\tau
D_{h,\textrm{sym}}(X^{n},\tilde{X}^{n})+\tau^{2}\lVert\nabla h(X^{n})-\nabla
h(\tilde{X}^{n})\rVert^{2}.$

We can use the strong convexity of $h$, which guarantees

$D_{h,\textrm{sym}}(X^{n},\tilde{X}^{n})\geq\frac{\mu_{h}L}{\mu_{h}+L}\lVert
X^{n}-\tilde{X}^{n}\rVert^{2}+\frac{1}{\mu_{h}+L}\lVert\nabla h(X^{n})-\nabla
h(\tilde{X}^{n})\rVert^{2},$

and we obtain

$\lVert
X^{n+1}-\tilde{X}^{n+1}\rVert^{2}\leq\left(1-\tau\frac{2\mu_{h}L}{\mu_{h}+L}\right)\lVert
X^{n}-\tilde{X}^{n}\rVert^{2}+\tau\left(\tau-\frac{2}{\mu_{h}+L}\right)\lVert\nabla
h(X^{n})-\nabla h(\tilde{X}^{n})\rVert^{2}.$

If the step size is chosen to be $\tau=2/(\mu_{h}+L)$, we obtain the
convergence rate

$\lVert
X^{n}-\tilde{X}^{n}\rVert^{2}\leq\left(\frac{L-\mu_{h}}{L+\mu_{h}}\right)^{2n}\lVert
X^{0}-\tilde{X}^{0}\rVert^{2}.$

Denoting the distribution of $X^{n}$ by $\mu^{n}$ and that of $\tilde{X}^{n}$
by $\tilde{\mu}^{n}$, we obtain by taking the infimum over all couplings of
$\mu^{0}$ and $\tilde{\mu}^{0}$ on the right side that

$\mathcal{W}_{2}(\mu_{n},\tilde{\mu}_{n})\leq\left(\frac{L-\mu_{h}}{L+\mu_{h}}\right)^{n}\mathcal{W}_{2}(\mu_{0},\tilde{\mu}_{0}).$

This shows that ULA is stable and will converge to its unique stationary
distribution. Similar computations can be carried out for other time-
discretized schemes based on Langevin diffusion, e.g. time-implicit leading to
proximal sampling schemes. Since the stationary distribution of the time-
discrete scheme is not the target distribution $\mu^{\ast}$, the missing step
to showing convergence of the Markov chain is to couple a time-discrete chain
and a time-continuous process and then finding an error bound. Typically, this
is more challenging and the introduced bias (in Wasserstein distance) between
the two settings scales as $\sqrt{\tau}$.

We now turn to the primal-dual setting of (2.4). Consider again two Markov
chains with potentially different initial values $(X^{0},Y^{0})\sim\pi^{0}$
and $(\tilde{X}^{0},\tilde{Y}^{0})\sim\tilde{\pi}^{0}$, but equal stochastic
input:

$\begin{array}[]{lcl}\xi^{n+1}\sim\mathrm{N}(0,I)&&\\\
Y^{n+1}=\operatorname{prox}_{\sigma f^{\ast}}(Y^{n}+\sigma
KX_{\theta}^{n})&&\tilde{Y}^{n+1}=\operatorname{prox}_{\sigma
f^{\ast}}(\tilde{Y}^{n}+\sigma K\tilde{X}_{\theta}^{n})\\\
X^{n+1}=\operatorname{prox}_{\tau g}(X^{n}-\tau
K^{T}Y^{n+1})+\sqrt{2\tau}\xi^{n+1}&&\tilde{X}^{n+1}=\operatorname{prox}_{\tau
g}(\tilde{X}^{n}-\tau K^{T}\tilde{Y}^{n+1})+\sqrt{2\tau}\xi^{n+1}\\\
X_{\theta}^{n+1}=X^{n+1}+\theta(X^{n+1}-X^{n})&&\tilde{X}_{\theta}^{n+1}=\tilde{X}^{n+1}+\theta(\tilde{X}^{n+1}-\tilde{X}^{n})\end{array}$
(4.1)

By convention, we set $X^{-1}=X^{0}$ and $\tilde{X}^{-1}=\tilde{X}^{0}$. In
the following, we denote the distribution of $(X^{n},Y^{n})$ by $\pi^{n}$ and
the $X^{n}$-marginal by $\mu^{n}$. Denote further $L\coloneqq\lVert K\rVert$
where $\lVert\cdot\rVert$ is the operator norm induced by the Euclidean norm
on $\mathcal{X}$ and $\mathcal{Y}$. We assume that also in the discrete
setting, the step sizes are coupled via the relation $\lambda=\sigma/\tau$.

We start by giving a stability result for the case where neither $g$ nor
$f^{\ast}$ are strongly convex.

######  Theorem 4.1.

Let $(X^{n},Y^{n})$ and $(\tilde{X}^{n},\tilde{Y}^{n})$ be given by the
coupling (4.1) with relaxation parameter $\theta=1$ and $\tau\sigma L^{2}\leq
1$. Then for the primal variable it holds

$\mathcal{W}_{2}(\mu^{N},\tilde{\mu}^{N})\leq\mathcal{T}_{1,\lambda^{-1}}(\pi^{0},\tilde{\pi}^{0}),$

where $\mathcal{T}_{1,\lambda^{-1}}$ is the transport distance over the joint
primal-dual variable defined in Section 2.3. If $\tau\sigma L^{2}<1$, then the
joint variable obeys

$\mathcal{T}_{1,\lambda^{-1}}(\pi^{N},\tilde{\pi}^{N})\leq
C\mathcal{T}_{1,\lambda^{-1}}(\pi^{0},\tilde{\pi}^{0}),$

with the constant $C=(1-\tau\sigma L^{2})^{-1/2}$.

In the case of strongly convex $g$ and $f^{\ast}$, we can prove the following
contraction result.

######  Theorem 4.2.

Assume that $g$ and $f^{\ast}$ are $\omega_{g}$\- and
$\omega_{f^{\ast}}$-strongly convex respectively. Let $(X^{n},Y^{n})$ and
$(\tilde{X}^{n},\tilde{Y}^{n})$ be given by the coupling (4.1) with
extrapolation parameter $\theta$ satisfying

$\max\\{(1+2\omega_{g}\tau)^{-1},(1+2\omega_{f^{\ast}}\sigma)^{-1}\\}\leq\theta<1,$

and further $\theta\sigma\tau L^{2}\leq 1$. Then with the constant
$\tilde{\lambda}\coloneqq\frac{1+2\omega_{g}\tau}{1+2\omega_{f^{\ast}}\sigma}\,\lambda$,
it holds

$\mathcal{W}_{2}(\mu^{N},\tilde{\mu}^{N})\leq\theta^{N/2}\mathcal{T}_{1,\tilde{\lambda}^{-1}}(\pi^{0},\tilde{\pi}^{0}).$

and if $\theta\sigma\tau L^{2}<1$, then with
$\tilde{C}\coloneqq(1-\theta\tau\sigma L^{2})^{-1/2}$ we have for the joint
variable

$\mathcal{T}_{1,\tilde{\lambda}^{-1}}(\pi^{N},\tilde{\pi}^{N})\leq\tilde{C}\theta^{N/2}\mathcal{T}_{1,\tilde{\lambda}^{-1}}(\pi^{0},\tilde{\pi}^{0}).$

It follows directly that in the latter case a stationary distribution
$\pi^{\text{{stat}}}$ exists and the distribution $\pi^{N}$ of samples
converges to $\pi^{\text{{stat}}}$ as $N\to\infty$.

In order to prove Theorems 4.1 and 4.2, we need the following auxiliary
inequality, which has been used several times in similar forms for primal-dual
optimization algorithms and is not specific to the sampling scheme here, see
e.g. [29].

######  Lemma 4.3.

Let $(X^{n},Y^{n})$ and $(\tilde{X}^{n},\tilde{Y}^{n})$, $n\geq 1$ be given by
the coupling (4.1) with arbitrary extrapolation parameter $\theta$. Using the
definition of the proximal mapping, we know that

$\displaystyle q^{n+1}$
$\displaystyle\coloneqq(Y^{n}-Y^{n+1})/\sigma+KX_{\theta}^{n}$ $\displaystyle
p^{n+1}$
$\displaystyle\coloneqq(X^{n}-X^{n+1}+\sqrt{2\tau}\xi^{n+1})/\tau-K^{T}Y^{n+1},$

define analogously $\tilde{p}^{n+1},\tilde{q}^{n+1}$ for the second chain
$(\tilde{X}^{n},\tilde{Y}^{n})$. Then it holds

$\displaystyle 2\langle
X^{n+1}-\tilde{X}^{n+1},p^{n+1}-\tilde{p}^{n+1}\rangle+\lVert
X^{n+1}-\tilde{X}^{n+1}\rVert_{\tau^{-1}}^{2}+\lVert
X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}$
$\displaystyle\quad+2\langle
Y^{n+1}-\tilde{Y}^{n+1},q^{n+1}-\tilde{q}^{n+1}\rangle+\lVert
Y^{n+1}-\tilde{Y}^{n+1}\rVert_{\sigma^{-1}}^{2}$ $\displaystyle\quad+2\langle
K(X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}),Y^{n+1}-\tilde{Y}^{n+1}\rangle+(1-\sigma\tau\theta
L^{2})\lVert
Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\leq\lVert
X^{n}-\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}+\theta\lVert
X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n}-\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\quad+2\theta\langle
K(X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}),Y^{n}-\tilde{Y}^{n}\rangle.$

###### Proof.

By plugging in the representation of the subgradient $p^{n+1}$ of $g$ at
$X^{n+1}$ (and likewise for the second chain) given in the statement of the
lemma, we can rewrite the symmetric Bregman divergence of $g$ at
$(X^{n+1},\tilde{X}^{n+1})$ as

$\displaystyle 2\langle
X^{n+1}-\tilde{X}^{n+1},p^{n+1}-\tilde{p}^{n+1}\rangle$
$\displaystyle=2\langle
X^{n+1}-\tilde{X}^{n+1},(X^{n}-X^{n+1})/\tau-K^{T}Y^{n+1}-(\tilde{X}^{n}-\tilde{X}^{n+1})/\tau+K^{T}\tilde{Y}^{n+1}\rangle$
$\displaystyle=-\frac{2}{\tau}\langle
X^{n+1}-\tilde{X}^{n+1},X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}\rangle-2\langle
X^{n+1}-\tilde{X}^{n+1},K^{T}(Y^{n+1}-\tilde{Y}^{n+1})\rangle$
$\displaystyle=\lVert X^{n}-\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}-\lVert
X^{n+1}-\tilde{X}^{n+1}\rVert_{\tau^{-1}}^{2}-\lVert
X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}$
$\displaystyle\quad-2\langle
K(X^{n+1}-\tilde{X}^{n+1}),Y^{n+1}-\tilde{Y}^{n+1}\rangle,$

where we used the identity $-\langle A-B,A\rangle=\frac{1}{2}\lVert
B\rVert^{2}-\frac{1}{2}\lVert A\rVert^{2}-\frac{1}{2}\lVert A-B\rVert^{2}$ in
the last equality. Repeating the same calculation for the subgradients
$q^{n+1}$ of $f^{\ast}$ at $Y^{n+1}$, we further get

$\displaystyle 2\langle
Y^{n+1}-\tilde{Y}^{n+1},q^{n+1}-\tilde{q}^{n+1}\rangle$
$\displaystyle=2\langle
Y^{n+1}-\tilde{Y}^{n+1},(Y^{n}-Y^{n+1})/\sigma+KX_{\theta}^{n}-(\tilde{Y}^{n}-\tilde{Y}^{n+1})/\sigma-K\tilde{X}_{\theta}^{n}\rangle$
$\displaystyle=-\frac{2}{\sigma}\langle
Y^{n+1}-\tilde{Y}^{n+1},Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rangle+2\langle
Y^{n+1}-\tilde{Y}^{n+1},K(X_{\theta}^{n}-\tilde{X}_{\theta}^{n})\rangle$
$\displaystyle=\lVert Y^{n}-\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}-\lVert
Y^{n+1}-\tilde{Y}^{n+1}\rVert_{\sigma^{-1}}^{2}-\lVert
Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\quad+2\langle
K(X_{\theta}^{n}-\tilde{X}_{\theta}^{n}),Y^{n+1}-\tilde{Y}^{n+1}\rangle.$

In order to shorten the expressions, we introduce the shorthand notation
$U^{n}:=X^{n}-\tilde{X}^{n}$,
$U_{\theta}^{n}=X_{\theta}^{n}-\tilde{X}_{\theta}^{n}$ and
$V^{n}:=Y^{n}-\tilde{Y}^{n}$. Adding up the two above equalities gives

$\displaystyle 2\langle U^{n+1},p^{n+1}-\tilde{p}^{n+1}\rangle+2\langle
V^{n+1},q^{n+1}-\tilde{q}^{n+1}\rangle$ $\displaystyle=\lVert
U^{n}\rVert_{\tau^{-1}}^{2}-\lVert U^{n+1}\rVert_{\tau^{-1}}^{2}-\lVert
U^{n+1}-U^{n}\rVert_{\tau^{-1}}^{2}$ $\displaystyle\quad+\lVert
V^{n}\rVert_{\sigma^{-1}}^{2}-\lVert V^{n+1}\rVert_{\sigma^{-1}}^{2}-\lVert
V^{n+1}-V^{n}\rVert_{\sigma^{-1}}^{2}$ $\displaystyle\quad-2\langle
KU^{n+1},V^{n+1}\rangle+2\langle KU_{\theta}^{n},V^{n+1}\rangle.$ (4.2)

By definition it holds $U_{\theta}^{n}=U^{n}+\theta(U^{n}-U^{n-1})$, so the
two inner products in the last line can be rewritten in the following way

$\displaystyle-2\langle KU^{n+1},V^{n+1}\rangle+2\langle
KU_{\theta}^{n},V^{n+1}\rangle$ $\displaystyle=-2\langle
K(U^{n+1}-U^{n}),V^{n+1}\rangle+2\theta\langle
K(U^{n}-U^{n-1}),V^{n+1}\rangle$ $\displaystyle=-2\langle
K(U^{n+1}-U^{n}),V^{n+1}\rangle+2\theta\langle
K(U^{n}-U^{n-1}),V^{n}\rangle+2\theta\langle
K(U^{n}-U^{n-1}),V^{n+1}-V^{n}\rangle$ $\displaystyle\leq-2\langle
K(U^{n+1}-U^{n}),V^{n+1}\rangle+2\theta\langle K(U^{n}-U^{n-1}),V^{n}\rangle$
$\displaystyle\quad+\theta\lVert
U^{n}-U^{n-1}\rVert_{\tau^{-1}}^{2}+\tau\sigma\theta L^{2}\lVert
V^{n+1}-V^{n}\rVert_{\sigma^{-1}}^{2}.$

The inequality follows from applying Cauchy-Schwarz and the inequality $ab\leq
a^{2}/2c+cb^{2}/2$ for any $c>0$ (Young’s inequality, with $a=\lVert
K(U^{n}-U^{n-1})\rVert$ and $b=\lVert V^{n+1}-V^{n}\rVert$ and $c=\tau L^{2}$
here) to the last inner product in the expression before. Inserting this bound
into (4) gives precisely the statement of the Lemma. ∎

###### Proof of Theorem 4.1.

In the convex case, we start by noting that symmetric Bregman distances for
convex functions are always non-negative, i.e.

$\displaystyle\langle X^{n+1}-\tilde{X}^{n+1},p^{n+1}-\tilde{p}^{n+1}\rangle$
$\displaystyle\geq 0,$ $\displaystyle\langle
Y^{n+1}-\tilde{Y}^{n+1},q^{n+1}-\tilde{q}^{n+1}\rangle$ $\displaystyle\geq 0.$

Setting $\theta=1$, the inequality in Lemma 4.3 hence gives

$\displaystyle\lVert X^{n+1}-\tilde{X}^{n+1}\rVert_{\tau^{-1}}^{2}+\lVert
X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n+1}-\tilde{Y}^{n+1}\rVert_{\sigma^{-1}}^{2}$ $\displaystyle\quad+2\langle
K(X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}),Y^{n+1}-\tilde{Y}^{n+1}\rangle+(1-\sigma\tau
L^{2})\lVert
Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\leq\lVert X^{n}-\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}+\lVert
X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n}-\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$ $\displaystyle\quad+2\langle
K(X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}),Y^{n}-\tilde{Y}^{n}\rangle.$

Iterating this inequality over $n=0,\dots,N-1$ steps and using $X^{-1}=X^{0}$,
$\tilde{X}^{-1}=X^{0}$ gives

$\displaystyle\lVert X^{N}-\tilde{X}^{N}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{N}-\tilde{Y}^{N}\rVert_{\sigma^{-1}}^{2}+\lVert
X^{N}-\tilde{X}^{N}-X^{N-1}+\tilde{X}^{N-1}\rVert_{\tau^{-1}}^{2}$
$\displaystyle\quad+2\langle
K(X^{N}-\tilde{X}^{N}-X^{N-1}+\tilde{X}^{N-1}),Y^{N}-\tilde{Y}^{N}\rangle+\left(1-\tau\sigma
L^{2}\right)\sum_{n=0}^{N-1}\lVert
Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\leq\lVert X^{0}-\tilde{X}^{0}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{0}-\tilde{Y}^{0}\rVert_{\sigma^{-1}}^{2}$

Using Cauchy-Schwarz and Young’s inequality on the remaining inner product and
non-negativity of the sum of norms implies

$\displaystyle\lVert X^{N}-\tilde{X}^{N}\rVert_{\tau^{-1}}^{2}+(1-\tau\sigma
L^{2})\lVert Y^{N}-\tilde{Y}^{N}\rVert_{\sigma^{-1}}^{2}\leq\lVert
X^{0}-\tilde{X}^{0}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{0}-\tilde{Y}^{0}\rVert_{\sigma^{-1}}^{2}.$

Taking expectations on both sides, by definition of the Wasserstein distance,
the left side is lower bounded by
$\tau^{-1}\mathcal{W}_{2}^{2}(\mu^{N},\tilde{\mu}^{N})$. Since the resulting
inequality holds for all couplings of initializations $(X^{0},Y^{0})$ and
$(\tilde{X}^{0},\tilde{Y}^{0})$, taking the infimum over all couplings on the
right side gives
$\tau^{-1}\mathcal{T}^{2}_{1,\lambda^{-1}}(\pi^{0},\tilde{\pi}^{0})$, proving
the first inequality. The second inequality in the statement of Theorem 4.1
results from lower bounding the left side not by the Wasserstein distance in
the primal variable, but instead by
$C^{-2}\tau^{-1}\mathcal{T}^{2}_{1,\lambda^{-1}}(\pi^{N},\tilde{\pi}^{N})$. ∎

###### Proof of Theorem 4.2.

We now assume strong convexity of $g,f^{\ast}$. In that case, the symmetric
Bregman distances obey the following lower bounds

$\displaystyle\langle X^{n+1}-\tilde{X}^{n+1},p^{n+1}-\tilde{p}^{n+1}\rangle$
$\displaystyle\geq\omega_{g}\lVert X^{n+1}-\tilde{X}^{n+1}\rVert^{2},$
$\displaystyle\langle Y^{n+1}-\tilde{Y}^{n+1},q^{n+1}-\tilde{q}^{n+1}\rangle$
$\displaystyle\geq\omega_{f^{\ast}}\lVert Y^{n+1}-\tilde{Y}^{n+1}\rVert^{2}.$

The inequality in Lemma 4.3 now implies

$\displaystyle\lVert
X^{n+1}-\tilde{X}^{n+1}\rVert_{\tau^{-1}+2\omega_{g}}^{2}+\lVert
X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n+1}-\tilde{Y}^{n+1}\rVert_{\sigma^{-1}+2\omega_{f^{\ast}}}^{2}$
$\displaystyle\quad+2\langle
K(X^{n+1}-\tilde{X}^{n+1}-X^{n}+\tilde{X}^{n}),Y^{n+1}-\tilde{Y}^{n+1}\rangle+(1-\sigma\tau\theta
L^{2})\lVert
Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\leq\lVert
X^{n}-\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}+\theta\lVert
X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n}-\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\quad+2\theta\langle
K(X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}),Y^{n}-\tilde{Y}^{n}\rangle.$

Define now

$\displaystyle\Delta^{n}$ $\displaystyle\coloneqq\lVert
X^{n}-\tilde{X}^{n}\rVert_{\tau^{-1}+2\omega_{g}}^{2}+\lVert
X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n}-\tilde{Y}^{n}\rVert_{\sigma^{-1}+2\omega_{f^{\ast}}}^{2}$
$\displaystyle\quad 2\langle
K(X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}),Y^{n}-\tilde{Y}^{n}\rangle.$

Then under the assumptions on $\theta$ given in Theorem 4.2, it holds

$\displaystyle\theta\Delta^{n}$ $\displaystyle\geq\lVert
X^{n}-\tilde{X}^{n}\rVert_{\tau^{-1}}^{2}+\theta\lVert
X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}\rVert_{\tau^{-1}}^{2}+\lVert
Y^{n}-\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$ $\displaystyle\quad
2\theta\langle
K(X^{n}-\tilde{X}^{n}-X^{n-1}+\tilde{X}^{n-1}),Y^{n}-\tilde{Y}^{n}\rangle$
$\displaystyle\geq\Delta^{n+1}+(1-\sigma\tau\theta L^{2})\lVert
Y^{n+1}-\tilde{Y}^{n+1}-Y^{n}+\tilde{Y}^{n}\rVert_{\sigma^{-1}}^{2}$
$\displaystyle\geq\Delta^{n+1}$

Iterating this inequality over $n=0,\dots,N-1$ we obtain

$\displaystyle\Delta^{0}$ $\displaystyle\geq\theta^{-N}\Delta^{N}.$

We use again Cauchy-Schwarz and Young’s inequality on the inner product in
$\Delta^{N}$ as well as the initial condition $X^{-1}=X^{0}$,
$\tilde{X}^{-1}=X^{0}$ and obtain

$\displaystyle\lVert
X^{N}-\tilde{X}^{N}\rVert_{\tau^{-1}+2\omega_{g}}^{2}+(1-\theta\tau\sigma
L^{2})\lVert Y^{N}-\tilde{Y}^{N}\rVert_{\sigma^{-1}+2\omega_{f^{\ast}}}^{2}$
$\displaystyle\quad\leq\theta^{N}\left(\lVert
X^{0}-\tilde{X}^{0}\rVert_{\tau^{-1}+2\omega_{g}}^{2}+\lVert
Y^{0}-\tilde{Y}^{0}\rVert_{\sigma^{-1}+2\omega_{f^{\ast}}}^{2}\right)$

From here, we proceed as in the proof of Theorem 4.1: Taking expectations on
both sides, the left side can be bounded below by
$(\tau^{-1}+2\omega_{g})W_{2}^{2}(\mu^{N},\tilde{\mu}^{N})$ or alternatively
in the joint variable by
$\tilde{C}^{-2}(\tau^{-1}+2\omega_{g})\mathcal{T}_{1,\tilde{\lambda}^{-1}}^{2}(\pi^{N},\tilde{\pi}^{N})$.
The transport distance on the right is obtained by taking the infimum over all
initialization couplings, completing the proof. ∎

######  Remark 4.4.

In continuous time, our theory currently provides us with a bound between the
stationary distribution and the target $\mu^{\ast}$ via Theorem 3.11. In
discrete time, we obtained bounds between the stationary distribution and the
samples’ law. In order to make convergence statements between the samples and
the target, we are currently working on a final coupling result to bridge the
gap between discrete and continuous time stationary solutions. It will be
reported in a future version of this manuscript.

##  5 Numerical Examples

In the following, we demonstrate the behaviour and our new results on ULPDA on
several numerical examples, first small-scale ones to demonstrate the
convergence behaviour in detail and subsequently larger imaging inverse
problems. We always simulate (2.5) in a time interval $[0,T]$ by running the
discretized algorithm (2.4) with small step sizes $\tau,\sigma$ for $T/\tau$
steps. The step sizes are coupled via $\lambda\tau=\sigma$ and are chosen
small enough so that $\tau\sigma L^{2}\leq 1$, i.e. the assumptions of Theorem
4.1 are satisfied and the continuous-time dynamics are approximated well.

###  5.1 Revisiting the Quadratic Case

We start by illustrating the toy example of two one-dimensional quadratic
potentials that was considered in Theorems 3.5 and 3.6. Recall that the target
distribution has density $\exp(-h(x))/Z$ with the potential $h$ defined as

$h(x)=f(kx)+g(x)=\frac{k^{2}}{2c_{f}}x^{2}+\frac{1}{2c_{g}}x^{2}=\frac{c_{g}k^{2}+c_{f}}{2c_{f}c_{g}}x^{2},$

hence the target distribution is Gaussian
$\mathrm{N}(0,c_{f}c_{g}(c_{g}k^{2}+c_{f})^{-1})$. For the numerical
simulation, we choose $c_{f}=1,c_{g}=2,k=1.5$. We simulate the primal-dual
Langevin dynamics (2.5) for different values of $\lambda$ by running ULPDA
(2.4). We set $\theta=1$ and set the step sizes $\tau$ and $\sigma$ via the
relations $\sigma/\tau=\lambda$ and $\sigma\tau k^{2}=c$ where $c\leq 1$ is
required for stability. We set $c=10^{-4}$ in order to approximate closely the
continuous time dynamics, resulting in step sizes
$\tau=10^{-2}\sqrt{\lambda}k$ and $\sigma=10^{-2}\sqrt{\lambda^{-1}}k$. Figure
2 shows the joint distribution of primal and dual samples $(X^{n},Y^{n})$
produced by (2.4). As $\lambda$ grows, the joint distribution of primal and
dual variable concentrates along the linear subspace $y=kx/c_{f}$ and the
primal marginal distribution converges to the target.

Figure 2: Approximation of the target $\mu^{\ast}=\exp(-h)$ for a one-
dimensional quadratic $h$ by (2.4) and its continuous time diffusion limit
(2.5). From left to right, results for $\lambda=1$, $\lambda=10$ and
$\lambda=100$ are displayed. The center scatter plot shows pairs
$(X^{n},Y^{n})$, the background color indicates the mass of the concentrated
joint target $\pi^{\ast}$ (zero is blue, the mass is concentrated along the
white line $y=kx/c_{f}$). The marginal boxes show the densities of the primal
and dual targets $\mu^{\ast}$ and $\nu^{\ast}=(\partial f\circ
K)_{\\#}\mu^{\ast}$ (dashed red), the marginals
$\mu^{\text{stat}}_{\text{{cPD}}},\nu^{\text{stat}}_{\text{{cPD}}}$ of the
invariant solution of primal-dual diffusion (black solid) and samples drawn
using (2.4) (grey histograms). For increasing $\lambda$, the joint
distribution of samples tends to concentrate and the primal marginal converges
to the target, as proven in Theorem 3.11.

###  5.2 TV Denoising from Gaussian Noise

For the potential $f\circ K$ being is dualized in the algorithm, we now
consider the total variation (TV) functional, a typical model-based prior in
imaging inverse problems. The (isotropic) TV functional is defined by

$\operatorname{TV}(x):=\sum_{i,j}\sqrt{(\nabla_{h}x)_{i,j}^{2}+(\nabla_{v}x)_{i,j}^{2}},$

where $\nabla_{h},\nabla_{v}$ denote horizontal and vertical finite
differences in the pixelated image $x$. We consider a Gibbs prior on $x$ with
density $\exp(-\alpha\operatorname{TV}(x))$ where $\alpha$ is a fixed
regularization parameter (Note that this does only defines a probability
distribution on a subspace orthogonal to constant functions since TV has a
one-dimensional null space. With a log-likelihood that grows fast enough at
infinity along this null space, the posterior is however well-defined.). The
corresponding potential term can be written in the form
$\alpha\operatorname{TV}(x)=\alpha\lVert Kx\rVert_{2,1}=f(Kx)$ where
$K=(\nabla_{h},\nabla_{v})$ and $f=\alpha\lVert\cdot\rVert_{2,1}$ denotes a
constant multiple of the nested $l_{2}$-$l_{1}$-norm (inner $l_{2}$ over
vertical/horizontal component, outer $l_{1}$ over all pixels). We further
assume that we make a noisy observation $\tilde{x}=x+\varepsilon$ of $x$,
where the noise is distributed like
$\varepsilon\sim\mathrm{N}(0,\sigma_{\varepsilon}^{2}I)$ for some observation
variance $\sigma_{\varepsilon}^{2}$. The corresponding negative log-likelihood
is

$g(x):=l_{\tilde{x}}(x)=\frac{1}{2\sigma_{\varepsilon}^{2}}\lVert
x-\tilde{x}\rVert_{2}^{2},$

which is $\sigma^{-2}$-strongly convex. We now want to sample from the
posterior $\mu^{\ast}$, which by Bayes law has density proportional to
$\exp(-\alpha\lVert Kx\rVert_{2,1}-l_{\tilde{x}}(x))=\exp(-f(Kx)-g(x))$.

####  5.2.1 Toy Example in Two Dimensions

We start with a small-scale example in which we are able to compute accurate
approximations of Wasserstein distances between the samples’ distribution and
the true posterior. Consider the artificial case where the ‘image’ consists of
$2\times 1$ pixels, i.e. $\mathcal{X}=\mathbb{R}^{2}$. Then the operator $K$
reduces to $Kx=x_{2}-x_{1}\in\mathbb{R}=\mathcal{Y}$ with $\lVert
K\rVert=\sqrt{2}$, and $f$ is given by $f(y)=\alpha\left|y\right|$. The convex
conjugate $f^{\ast}$ is the indicator function of the interval
$[-\alpha,\alpha]$ and the proximal mapping $\operatorname{prox}_{\sigma
f^{\ast}}$ is, for any $\sigma>0$, the projection onto $[-\alpha,\alpha]$.

The distribution $\mu^{n}$ of iterate $X^{n}$ in ULPDA is approximated by
running finitely many, here $10^{4}$, independent Markov chains in parallel
and using the resulting empirical distribution $\hat{\mu}^{n}$. In Figure 3,
we show scatter plots of the samples at time $T$, where $T$ is chosen large
enough that the empirical measures $\hat{\mu}^{n}$ remain stationary both
qualitatively and in Wasserstein distance to the target (see below). As
comparison, we show samples drawn using the algorithm Prox-Sub proposed in
[3]. Prox-Sub solves the optimality condition $y\in\partial f(Kx)$ exactly at
every step (by assuming that the subgradient can be evaluated), and hence
corresponds exactly to the concentrated limit $\lambda=\infty$ of ULPDA
explored in Section 3.3. The authors of [3] showed the convergence of the
distribution of Prox-Sub samples to the target in 2-Wasserstein distance,
total variation distance and Kullback-Leibler divergence, where the bias
decreases depends as usual on the step size of the algorithm.

Figure 3: Scatter plots of $\hat{\mu}^{n}$ in the two-dimensional TV-denoising
experiment. From left to right, the plots show approximations of the same
posterior with values $\lambda=10^{1},10^{2},10^{3}$ and ‘$\lambda=\infty$’
(Prox-Sub). Horizontal and vertical axes correspond to the two coordinates in
the primal variable $x\in\mathbb{R}^{2}$, the dual variable $y$ is not
displayed. Grey dots show the samples at time $T$ produced by $10^{4}$
parallel chains, the black arrows indicate the samples’ mean and the principal
components of their covariance (i.e. eigenvectors of the covariance matrix
scaled in length by the corresponding eigenvalues). For reference, the blue
arrows are the principal components of the posterior that should be
approximated. The noisy observation $\tilde{x}$ is marked by a yellow triangle
and the background color is a contour plot of the posterior density
$\mu^{\ast}(x)$.

Using numerical integration, in this low-dimensional example we are able to
compute an estimate of the normalization constant
$Z=\int_{\mathcal{X}}\exp(-f(Kx)-g(x))\,\mathrm{d}x$. Using this, we can
access an approximation of the posterior density $\mu^{\ast}(x)$. This allows
us to compute approximations of the Wasserstein distances of current iterates
to the target $\mathcal{W}_{2}(\hat{\mu}^{n},\mu^{\ast})$ without the need to
rely on a second sampler. We evaluate the distance at 20 equidistantly spaced
points in $[0,T]$ (i.e. for $\hat{\mu}^{kN},k=1,\dots,20$ where $N=T/20\tau$).
The distances are computed using the Python optimal transport module [12] and
are displayed in Figure 4. As expected, the samples drawn by the Prox-Sub
algorithm are closest to the target distribution (their bias in Wasserstein
distance is bounded by a constant multiple of the step size $\tau$). For
finite $\lambda$, the distributions $\hat{\mu}^{n}$ exhibit an additional
strong bias to the target due to the dualization, which becomes smaller as
$\lambda$ grows.

Figure 4: Wasserstein distances between the empirical measures $\hat{\mu}^{n}$
and the target $\mu^{\ast}$ in the two-dimensional TV-denoising experiment. We
compare three different values of $\lambda$ and the algorithm Prox-Sub [3]
corresponding formally to $\lambda=\infty$. The discretization via partial
dualization of $f\circ K$ introduces a bias between samples and target that is
significantly larger than the pure discretization bias of Prox-Sub. As
$\lambda$ grows, the iteration approximates the purely primal setting and the
bias is reduced.

####  5.2.2 Image Denoising

We now sample from the posterior of an image denoising task. Due to the higher
dimension, we can not compute accurate approximations of the Wasserstein
distances to the target anymore. We therefore resort to illustrating the
behaviour of ULPDA using moment estimates. A ground truth image normalized to
intensities $[0,1]$ is distorted by noise with standard deviation $0.25$ (PSNR
$12.02\,$dB in the noisy observation). We then draw $10^{6}$ samples with a
single chain and compute the sample mean and pixel-wise variance after a burn-
in phase of $1e3$ samples. The experiment is repeated three times for values
$\lambda=1,10,100$, each with the same step size $\tau$. A fourth chain is
computed using Prox-Sub with the same $\tau$. MMSE and pixel-wise variance
estimates are reported in Figure 5. The sample mean (MMSE) estimates produced
by the different chains converge quickly and are very similar for all chains
(PSNR values between $23.82\,$dB and $23.85\,$dB). The pixel-wise variance,
however, shows the same effect of overdispersion for smaller levels of
$\lambda$ as in the lower-dimensional example.

While the ‘primal’ samples $X^{n}$ are overdispersed for finite $\lambda$, the
dual iterates $Y^{n}$ characteristically exhibit underdispersion (see Figure
6). This effect was consistent throughout all our numerical experiments
(compare also Figures 2 and 3). A heuristic explanation for overdispersion in
the primal and underdispersion in the dual variable lies in the fact that the
step in the dual variable solves the optimality condition $y\in\partial f(Kx)$
only incompletely. On the particle level in continuous time, when the primal
variable $X_{t}$ is driven towards the tail regions of its marginal
distribution by the Brownian motion, the dual variable is not evolved entirely
to the corresponding position $\partial f(KX_{t})$ in the tails of its
marginal, but rather remains near the mode of its marginal. The effect
decreased with growing parameters $\lambda$ in all our experiments.

Figure 5: MMSE estimates and pixel-wise variance computed via the posterior
sample mean in the $l_{2}$-TV denoising experiment. Top row shows the MMSE
estimates, bottom row the corresponding pixel-wise (logarithmically scaled)
variances. From left to right results for $\lambda=1,10,100$ and
‘$\lambda=\infty$’ (using Prox-Sub) are shown.

Figure 6: Pixel-wise (logarithmic) standard deviation of the dual variables
$Y^{n}$, where an $l_{2}$-norm over horizontal and vertical component was
applied. From left to right, we show the results for $\lambda=1,10,100$ and
‘$\lambda=\infty$’ (Prox-Sub). For the latter case, the algorithm computes the
dual variable as $Y^{n}\in\partial f(KX^{n})$. For finite $\lambda$, the dual
samples suffer from underdispersion, the effect vanishes as $\lambda$ grows
and the dual values concentrate around the dual optimality condition.

###  5.3 Total Generalized Variation Based Image Denoising

In this experiment, we consider an example with an image prior based on second
order total generalized variation (TGV). Like the TV functional, total
generalized variation allows for discontinuous edges in images, but it
generalizes TV in that it also employs higher-order derivatives. This way, it
avoids the typical staircasing effect of TV in smooth but non-constant regions
of images. For optimization purposes, the second order TGV functional is
usually defined as

$\displaystyle\operatorname{TGV}^{2}_{\alpha}(u)=\inf_{v\in
V}J^{2}_{\alpha}(u,v),$ (5.1) $\displaystyle
J^{2}_{\alpha}(u,v)=\alpha_{1}\lVert\nabla u-v\rVert_{2,1}+\alpha_{0}\lVert
Ev\rVert_{2,1}.$

Here $u\in\mathcal{U}$ is an image while $v$ is in the corresponding space
$\mathcal{V}=\mathcal{U}\times\mathcal{U}$ of (discretized) gradients of $u$.
$\alpha_{0},\alpha_{1}>0$ are regularization parameters.
$\nabla=(\nabla_{h},\nabla_{v}):\mathcal{U}\to\mathcal{V}$ is a (forward
discretized) gradient operator as in the case of TV, and
$E:\mathcal{V}\to\mathcal{W}$ is a symmetrized (backward discretized) gradient
operator that acts as second order derivative. We refer to Section 2.1 of [26]
for details on the implementation of the gradient operators. In TGV-based
variational minimization, the state space is typically expanded from the image
space $\mathcal{U}$ to the joint space
$\mathcal{X}=\mathcal{U}\times\mathcal{V}$ and the minimization over $u$ (in
the variational problem) and $v$ (in 5.1) is carried out simultaneously. By
applying primal-dual optimization methods, the iteration can be carried out
efficiently without the need to evaluate the infimum in (5.1) explicitly [26].
In the same way, we want to avoid evaluating $\operatorname{TGV}^{2}_{\alpha}$
or its gradients explicitly in a sampling procedure. In order to recover the
image $u$, we therefore propose to sample from a posterior of the form

$\mu^{\ast}(x)\propto
p(\tilde{x}|x)\exp\left(-J^{2}_{\alpha}(u,v)\right),\qquad
x=(u,v)\in\mathcal{U}\times\mathcal{V}=\mathcal{X},$

where $p(\tilde{x}|x)$ is the likelihood, for denosing with normally
distributed noise given as before by

$-\log p(\tilde{x}|x)=\frac{1}{2\sigma_{\varepsilon}^{2}}\lVert
u-\tilde{u}\rVert_{\mathcal{U}}^{2},\qquad
x=(u,v),\tilde{x}=(\tilde{u},\tilde{v})\in\mathcal{X}$

We note that this is a relaxed version of regularization with
$\operatorname{TGV}_{\alpha}^{2}$, since we dropped the requirement that $v$
solves (5.1). Practically, the prior $\exp(-J^{2}_{\alpha}(u,v))$ seems to
yield very similar smoothing behaviour to the effect of total generalized
variation, we leave its detailed analysis to future work. Note that the chosen
target log-density is in the necessary form (1.1): $g(x)$ is, as before, the
negative log-likelihood $\lVert
u-\tilde{u}\rVert^{2}_{\mathcal{U}}/(2\sigma_{\varepsilon}^{2})$ for all
$x=(u,v)$. We define the linear operator

$K=\begin{pmatrix}\nabla&-I\\\
0&E\end{pmatrix}:\mathcal{X}\to\mathcal{V}\times\mathcal{W}\eqqcolon\mathcal{Y},$

where $I$ is the identity on $\mathcal{V}$. The functional
$f:\mathcal{Y}\to\mathbb{R}$ is given by $f(y)=\alpha_{1}\lVert
v\rVert_{2,1}+\alpha_{0}\lVert w\rVert_{2,1}$ for all $y=(v,w)$. Due to the
separable structure in $v$ and $w$, the proximal mapping of the convex
conjugate $f^{\ast}$ is simply a componentwise projection in the dual
$l_{2}$-$l_{\infty}$-norms in the respective spaces $\mathcal{V}$,
$\mathcal{W}$.

We run the experiment with empirically set values of $\lambda$ that make sure
the variance is not strongly overdispersed while keeping the computational
expense moderate – recall that as $\lambda$ grows, $\tau$ needs to be chosen
smaller, which slows down the convergence. $10^{5}$ samples are drawn in each
run, we report the resulting sample mean and pixel-wise variance in Figure 7.

Figure 7: Results from the TGV-based denoising experiment. Left: noisy
observation $\tilde{u}$, Center: MMSE estimate via sample mean, Right: Log-
scaled pixel-wise standard deviation. Top is a moderately noisy phantom
designed for TGV-regularization ($\sigma_{\varepsilon}=0.1$, PSNR of noisy
image $19.98\,$dB, MMSE PSNR $34.00\,$dB), bottom a more strongly distorted,
photorealistic image ($\sigma_{\varepsilon}=0.25$, PSNR of noisy image
$12.02\,$dB, MMSE PSNR $20.08\,$dB).

##  6 Conclusion and Open Problems

We have analyzed a recently proposed algorithm for sampling from log-concave,
non-smooth potentials which makes use of alternating primal-dual steps as is
well-known in convex non-smooth optimization. We showed that the sampling
algorithm corresponds to a stochastic differential equation and formulated the
mean field limit Fokker-Planck equation. Using a simple coupling argument,
contraction of the continuous-time limit to a stationary solution was proved.
Unlike for purely primal overdamped Langevin diffusion, the primal-dual
Fokker-Planck equation does not have the target distribution as its stationary
solution. A proposed correction of the Fokker-Planck equation forcing the
target to be its stationary solution is possible, but eventually not useful
for numerical implementation. However, the target distribution can be related
to the stationary solution of the uncorrected scheme in the limit of the ratio
of primal and dual step sizes growing to infinity. We made use of coupling
techniques to prove a corresponding bound in Wasserstein distance to the
target for a finite step size ratio. With similar coupling arguments in
discrete time, estimates in the typical spirit of convergence analysis of the
optimization algorithm can be transferred to the sampling case, showing
stability of the sampling method in Wasserstein distance.

We are currently working on a final coupling result, filling the gap between
the stationary solution in discrete time and the stationary solution in
continuous time. This, together with the stability in discrete time and the
bound between continuous-time stationary state and the target would allow a
convergence analysis between the algorithm and the target. Numerically, the
results suggest that if there is a way to compute the subgradient of the
dualized potential term, the additional bias due to the finite step size ratio
makes the primal-dual samples consistently over-dispersed compared to the
method of [3]. We currently see the primal-dual method’s use cases in two
settings. If the subgradient of the dualized term is not available, but the
proximal mapping of its convex conjugate is, the method is up to our knowledge
the first applicable algorithm able to draw samples from an approximation of
the posterior. Secondly, if the additional bias due to a larger finite step
size ratio does not matter (because a correction step is added or because
downstream tasks can be carried out with biased samples), the method can be
applied with potentially larger step sizes than a purely primal method,
allowing for faster mixing time. To find out whether this makes the primal-
dual method more efficient in applications, further numerical studies, in
particular in combination with a Metropolis-Hastings correction step, will be
needed.

## Acknowledgements

We thank Jan Modersitzki (Lübeck) for fruitful discussions and predicting non-
concentration of the dual variable.

## References

  * [1] Ziruo Cai, Jinglai Li and Xiaoqun Zhang “Approximate Primal-Dual Fixed-Point based Langevin Algorithms for Non-smooth Convex Potentials” arXiv, 2023 DOI: 10.48550/ARXIV.2304.04544
  * [2] Matthias J. Ehrhardt, Lorenz Kuger and Carola-Bibiane Schönlieb “Proximal Langevin Sampling With Inexact Proximal Mapping”, 2023 arXiv:2306.17737 [stat.CO]
  * [3] Andreas Habring, Martin Holler and Thomas Pock “Subgradient Langevin Methods for Sampling from Non-smooth Potentials” arXiv, 2023 DOI: 10.48550/ARXIV.2308.01417
  * [4] 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” arXiv, 2023 DOI: 10.48550/ARXIV.2308.09460
  * [5] Tim Tsz-Kit Lau, Han Liu and Thomas Pock “Non-Log-Concave and Nonsmooth Sampling via Langevin Monte Carlo Algorithms” arXiv, 2023 DOI: 10.48550/ARXIV.2305.15988
  * [6] Savvas Melidonis, Paul Dobson, Yoann Altmann, Marcelo Pereyra and Konstantinos Zygalakis “Efficient Bayesian Computation for Low-Photon Imaging Problems” In _SIAM Journal on Imaging Sciences_ 16.3 Society for Industrial & Applied Mathematics (SIAM), 2023, pp. 1195–1234 DOI: 10.1137/22m1502240
  * [7] Alain Durmus, Éric Moulines and Marcelo Pereyra “A Proximal Markov Chain Monte Carlo Method for Bayesian Inference in Imaging Inverse Problems: When Langevin Meets Moreau” In _SIAM Review_ 64.4 Society for Industrial & Applied Mathematics (SIAM), 2022, pp. 991–1028 DOI: 10.1137/22m1522917
  * [8] Tim Tsz-Kit Lau and Han Liu “Bregman Proximal Langevin Monte Carlo via Bregman-Moreau Envelopes” In _Proceedings of the 39th International Conference on Machine Learning_ 162, Proceedings of Machine Learning Research PMLR, 2022, pp. 12049–12077 URL: https://proceedings.mlr.press/v162/lau22a.html
  * [9] Ruilin Li, Molei Tao, Santosh S. Vempala and Andre Wibisono “The Mirror Langevin Algorithm Converges with Vanishing Bias” In _Proceedings of The 33rd International Conference on Algorithmic Learning Theory_ 167, Proceedings of Machine Learning Research PMLR, 2022, pp. 718–742 URL: https://proceedings.mlr.press/v167/li22b.html
  * [10] Dominik Narnhofer, Andreas Habring, Martin Holler and Thomas Pock “Posterior-Variance-Based Error Quantification for Inverse Problems in Imaging” arXiv, 2022 DOI: 10.48550/ARXIV.2212.12499
  * [11] Kwangjun Ahn and Sinho Chewi “Efficient constrained sampling via the mirror-Langevin algorithm” In _Advances in Neural Information Processing Systems_ , 2021 URL: https://openreview.net/forum?id=vh7qBSDZW3G
  * [12] Remi Flamary et al. “POT: Python Optimal Transport” In _Journal of Machine Learning Research_ 22.78, 2021, pp. 1–8 URL: http://jmlr.org/papers/v22/20-451.html
  * [13] Nicolas Fournier and Benoît Perthame “Transport distances for PDEs: the coupling method” In _EMS Surveys in Mathematical Sciences_ 7.1 European Mathematical Society - EMS - Publishing House GmbH, 2021, pp. 1–31 DOI: 10.4171/emss/35
  * [14] Adil Salim and Peter Richtarik “Primal Dual Interpretation of the Proximal Stochastic Gradient Langevin Algorithm” In _Advances in Neural Information Processing Systems_ 33 Curran Associates, Inc., 2020, pp. 3786–3796 URL: https://proceedings.neurips.cc/paper/2020/file/2779fda014fbadb761f67dd708c1325e-Paper.pdf
  * [15] Arnak S. Dalalyan and Avetik Karagulyan “User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient” In _Stochastic Processes and their Applications_ 129.12, 2019, pp. 5278–5311 DOI: 10.1016/j.spa.2019.02.016
  * [16] Alain Durmus, Szymon Majewski and Błażej Miasojedow “Analysis of Langevin Monte Carlo via Convex Optimization” In _J. Mach. Learn. Res._ 20.1, 2019, pp. 2666–2711
  * [17] Alain Durmus and Éric Moulines “High-dimensional Bayesian inference via the unadjusted Langevin algorithm” In _Bernoulli_ 25.4A Bernoulli Society for Mathematical StatisticsProbability, 2019, pp. 2854–2882 DOI: 10.3150/18-BEJ1073
  * [18] Espen Bernton “Langevin Monte Carlo and JKO splitting” In _Proceedings of the 31st Conference On Learning Theory_ 75, Proceedings of Machine Learning Research PMLR, 2018, pp. 1777–1798 URL: https://proceedings.mlr.press/v75/bernton18a.html
  * [19] Sébastien Bubeck, Ronen Eldan and Joseph Lehec “Sampling from a Log-Concave Distribution with Projected Langevin Monte Carlo” In _Discrete and Computational Geometry_ 59.4 Springer ScienceBusiness Media LLC, 2018, pp. 757–783 DOI: 10.1007/s00454-018-9992-1
  * [20] Alain Durmus, Éric Moulines and Marcelo Pereyra “Efficient Bayesian Computation by Proximal Markov Chain Monte Carlo: When Langevin Meets Moreau” In _SIAM Journal on Imaging Sciences_ 11.1 Society for IndustrialApplied Mathematics, 2018, pp. 473–506 DOI: 10.1137/16M1108340
  * [21] Andre Wibisono “Sampling as optimization in the space of measures: The Langevin dynamics as a composite optimization problem” In _Proceedings of the 31st Conference On Learning Theory_ PMLR, 2018, pp. 2093–3027
  * [22] Nicolas Brosse, Alain Durmus, Éric Moulines and Marcelo Pereyra “Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo” In _Proceedings of the 2017 Conference on Learning Theory_ 65, Proceedings of Machine Learning Research PMLR, 2017, pp. 319–342 URL: https://proceedings.mlr.press/v65/brosse17a.html
  * [23] Alain Durmus and Éric Moulines “Nonasymptotic convergence analysis for the unadjusted Langevin algorithm” In _The Annals of Applied Probability_ 27.3 Institute of Mathematical Statistics, 2017 DOI: 10.1214/16-aap1238
  * [24] Antonin Chambolle and Thomas Pock “An introduction to continuous optimization for imaging” In _Acta Numerica_ 25 Cambridge University Press (CUP), 2016, pp. 161–319 DOI: 10.1017/s096249291600009x
  * [25] Marcelo Pereyra “Proximal Markov chain Monte Carlo algorithms” In _Statistics and Computing_ 26.4, 2016, pp. 745–760 DOI: 10.1007/s11222-015-9567-4
  * [26] Kristian Bredies and Martin Holler “A TGV-Based Framework for Variational Image Decompression, Zooming, and Reconstruction. Part II: Numerics” In _SIAM Journal on Imaging Sciences_ 8.4 Society for Industrial & Applied Mathematics (SIAM), 2015, pp. 2851–2886 DOI: 10.1137/15m1023877
  * [27] Grigorios A. Pavliotis “Stochastic Processes and Applications” Springer New York, 2014 DOI: 10.1007/978-1-4939-1323-7
  * [28] Heinz H. Bauschke and Patrick L. Combettes “Convex Analysis and Monotone Operator Theory in Hilbert Spaces” Springer New York, 2011 DOI: 10.1007/978-1-4419-9467-7
  * [29] Antonin Chambolle and Thomas Pock “A First-Order Primal-Dual Algorithm for Convex Problems with Applications to Imaging” In _Journal of Mathematical Imaging and Vision_ 40.1, 2011, pp. 120–145
  * [30] Mark Girolami and Ben Calderhead “Riemann Manifold Langevin and Hamiltonian Monte Carlo Methods” In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ 73.2 Oxford University Press (OUP), 2011, pp. 123–214 DOI: 10.1111/j.1467-9868.2010.00765.x
  * [31] E.J. Candes, J. Romberg and T. Tao “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information” In _IEEE Transactions on Information Theory_ 52.2 Institute of ElectricalElectronics Engineers (IEEE), 2006, pp. 489–509 DOI: 10.1109/tit.2005.862083
  * [32] Richard Jordan, David Kinderlehrer and Felix Otto “The Variational Formulation of the Fokker–Planck Equation” In _SIAM Journal on Mathematical Analysis_ 29.1 Society for Industrial & Applied Mathematics (SIAM), 1998, pp. 1–17 DOI: 10.1137/s0036141096303359
  * [33] Gareth O. Roberts and Richard L. Tweedie “Exponential Convergence of Langevin Distributions and Their Discrete Approximations” In _Bernoulli_ 2.4 JSTOR, 1996, pp. 341 DOI: 10.2307/3318418
  * [34] Lars Hörmander “Hypoelliptic second order differential equations” In _Acta Mathematica_ 119.0 International Press of Boston, 1967, pp. 147–171 DOI: 10.1007/bf02392081