# The Mirror Langevin Algorithm Converges with Vanishing Bias

Ruilin Li Georgia Institute of Technology & Hudson River Trading. Email:
<EMAIL_ADDRESS>Molei Tao Georgia Institute of Technology, School of
Mathematics. Email<EMAIL_ADDRESS>Santosh S. Vempala Georgia Institute of
Technology, College of Computing. Email<EMAIL_ADDRESS>Andre Wibisono
Yale University, Department of Computer Science. Email:
<EMAIL_ADDRESS>
###### Abstract

The technique of modifying the geometry of a problem from Euclidean to Hessian
metric has proved to be quite effective in optimization, and has been the
subject of study for sampling. The Mirror Langevin Diffusion (MLD) is a
sampling analogue of mirror flow in continuous time, and it has nice
convergence properties under log-Sobolev or Poincare inequalities relative to
the Hessian metric, as shown by Chewi et al. (2020). In discrete time, a
simple discretization of MLD is the Mirror Langevin Algorithm (MLA) studied by
Zhang et al. (2020), who showed a biased convergence bound with a non-
vanishing bias term (does not go to zero as step size goes to zero). This
raised the question of whether we need a better analysis or a better
discretization to achieve a vanishing bias. Here we study the Mirror Langevin
Algorithm and show it indeed has a vanishing bias. We apply mean-square
analysis based on Li et al. (2019) and Li et al. (2021) to show the mixing
time bound for MLA under the modified self-concordance condition introduced by
Zhang et al. (2020).

##  1 Introduction

Suppose we wish to sample from a probability distribution $\nu(x)\propto
e^{-f(x)}$ supported on a convex set $\mathcal{X}\subseteq\mathbb{R}^{d}$
where $f\colon\mathcal{X}\to\mathbb{R}$ is differentiable. A popular algorithm
is the Unadjusted Langevin Algorithm (ULA), which is a basic discretization of
the Langevin Dynamics in continuous time:

$\displaystyle dX_{t}=-\nabla f(X_{t})\,dt+\sqrt{2}\,dW_{t}.$

Langevin Dynamics has an optimization interpretation as the gradient flow for
minimizing relative entropy (KL divergence) with respect to $\nu$ using the
Wasserstein metric $W_{2}$ in the space of probability distributions on
$\mathbb{R}^{d}$, starting from the seminal work of Jordan et al. (1998); see
also Wibisono (2018). In continuous time, Langevin Dynamics has convergence
guarantees in various distances, including $W_{2}$ distance, KL divergence, or
$\chi^{2}$-divergence, under various conditions, such as strong log-concavity,
or functional inequalities such as the Log-Sobolev Inequality (LSI) or
Poincaré inequality. In discrete time, ULA is a biased discretization of the
Langevin Dynamics, and it has an asymptotic bias which scales with the step
size. In particular, by setting a small enough step size, we can obtain a
mixing time bound of ULA which has inverse polynomial dependence on the error
threshold; see for example Dalalyan (2017); Durmus and Moulines (2017, 2019);
Dalalyan and Karagulyan (2019); Vempala and Wibisono (2019); Li et al. (2019);
Li et al. (2021).

In many settings, the problem of interest is non-smooth or constrained (e.g.,
the $L_{1}$ ball or a general polytope), and the basic Langevin algorithm does
not apply. In optimization, this is handled effectively (and elegantly) by
interior-point methods, which use a convex “barrier” function to define a non-
Euclidean metric. The resulting metric is given locally by the Hessian of the
barrier function. This method results is a convergence bound that scales as
$\sqrt{d}$ for linear and convex optimization.

It is natural to wonder whether such a modification of the geometry could be
useful for sampling. Early evidence of this is the Dikin walk, which replaces
the ball walk (a discrete-time implementation of constrained Brownian motion)
by using an ellipsoid at each step, defined by the Hessian of the logarithmic
barrier function. This walk was shown to converge in $\tilde{O}(md)$ steps for
uniformly sampling a polytope with $m$ facets in $d$ dimension (Kannan and
Narayanan, 2012). It was recently refined using a weighted barrier function to
improve the convergence time to $\tilde{O}(d^{2})$ (Laddha et al., 2020).

A related approach that also originated in optimization is mirror descent,
which uses a mirror map (the gradient of the barrier function) to change the
geometry favorably, and in the context of sampling, can be seen as a
generalization of the Langevin algorithm by changing the metric. For
constrained sampling, the Mirror Langevin Dynamics was introduced by Zhang et
al. (2020) using a mirror map to constrain the domain; see also Hsieh et al.
(2018) for a related earlier approach. Mirror Langevin Dynamics is the
Langevin dynamics for sampling from $\nu$ using the Hessian metric generated
by the mirror map. In continuous time, Mirror Langevin Dynamics has nice
convergence guarantees under an analogous notion of mirror Poincare inequality
relative to the Hessian metric, as shown by Chewi et al. (2020); see also
Appendix A for a review. In discrete time, the Mirror Langevin Algorithm (MLA)
is a simple discretization of MLD proposed by Zhang et al. (2020), who showed
a biased convergence analysis but with a non-vanishing bias (does not go to
$0$ with step size, but remains a constant). This raised the question of
whether we need a better analysis of MLA or a better discretization of MLD to
achieve a vanishing bias. This also led to an alternative discretization
proposed by Ahn and Chewi (2020) which achieves a vanishing bias, but requires
an exact simulation of the Brownian motion with changing covariance.

In this paper, we study the Mirror Langevin Algorithm and show that it indeed
has a vanishing bias. The tool we will use is the mean-square analysis
framework, proposed by Li et al. (2019) and then refined by Li et al. (2021);
the latter version will be used. It will help establish a biased convergence
analysis of MLA under relative smoothness, strong convexity and modified self-
concordance; these are a subset of the conditions assumed by Zhang et al.
(2020). We show that the bias of MLA with step size $h$ scales as $\sqrt{h}$;
this leads to a $\tilde{O}(d/\epsilon^{2})$ mixing time bound for MLA (see
Theorem 3.1 and Corollary 3.2).

##  2 Algorithm and Problem Set-Up

###  2.1 Problem set-up

Suppose we want to sample from a probability distribution $\nu$ supported on a
convex set $\mathcal{X}\subseteq\mathbb{R}^{d}$. We assume $\nu$ is absolutely
continuous with respect to the Lebesgue measure on $\mathbb{R}^{d}$ and has
density $\nu(x)\propto e^{-f(x)}$ for some differentiable
$f\colon\mathcal{X}\to\mathbb{R}$.

Let $\phi\colon\mathcal{X}\to\mathbb{R}$ be a twice-differentiable strictly
convex function which is of Legendre type (Rockafellar, 1970). This implies
$\nabla\phi(\mathcal{X})=\mathbb{R}^{d}$, and in particular the gradient map
$\nabla\phi\colon\mathcal{X}\to\mathbb{R}^{d}$ is bijective. We also have
$\nabla^{2}\phi(x)\succ 0$ for all $x\in\mathcal{X}$. Moreover, we require
that $\|\nabla\phi(x)\|\to\infty$ and $\nabla^{2}\phi(x)\to\infty$ as $x$
approaches the boundary of $\mathcal{X}$. Using the Hessian metric
$\nabla^{2}\phi$ on $\mathcal{X}$ will prevent the iterates from leaving the
domain $\mathcal{X}$. We call $\nabla\phi\colon\mathcal{X}\to\mathbb{R}^{d}$
the mirror map and $\mathcal{Y}=\nabla\phi(\mathcal{X})=\mathbb{R}^{d}$ the
dual space.

Let $\phi^{\ast}\colon\mathbb{R}^{d}\to\mathbb{R}$ be the dual function of
$\phi$, defined by $\phi^{\ast}(y)=\sup_{x\in\mathcal{X}}\langle
x,y\rangle-\phi(x)$. Recall
$\nabla\phi^{\ast}(y)=\arg\max_{x\in\mathcal{X}}\langle x,y\rangle-\phi(x)$,
and we have $\nabla\phi^{\ast}=(\nabla\phi)^{-1}$, so
$\nabla\phi(\nabla\phi^{*}(y))=y$ for all $y\in\mathbb{R}^{d}$. Furthermore,
$\nabla^{2}\phi(x)=\nabla^{2}\phi^{*}(\nabla\phi(x))^{-1}$ for all
$x\in\mathcal{X}$.

For a vector $v\in\mathbb{R}^{d}$, let $\|v\|=\sqrt{\langle v,v\rangle}$ be
the $\ell_{2}$-norm. For a matrix $A\in\mathbb{R}^{d\times d}$, let
$\|A\|_{\mathrm{HS}}=\sqrt{\operatorname{Tr}(AA^{\top})}$ be the Hilbert-
Schmidt norm.

###  2.2 Mirror Langevin Algorithm

In this paper we study the Mirror Langevin Algorithm:

$\displaystyle x_{k+1}=\nabla\phi^{*}\left(\nabla\phi(x_{k})-h\nabla
f(x_{k})+\sqrt{2h}\sqrt{\nabla^{2}\phi(x_{k})}\,z_{k}\right)$ (1)

where $h>0$ is step size and $z_{k}\sim\mathcal{N}(0,I)$ is an independent
Gaussian random variable in $\mathbb{R}^{d}$. Here $\sqrt{\nabla^{2}\phi(x)}$
is a square-root of $\nabla^{2}\phi(x)$, i.e. any matrix
$C(x)\in\mathbb{R}^{d\times d}$ satisfying
$C(x)C(x)^{\top}=\nabla^{2}\phi(x)$. This algorithm can be seen as a sampling
version of the mirror descent algorithm from optimization, since we can write
the update of MLA in the following form which resembles mirror descent:

$\displaystyle x_{k+1}=\arg\min_{x\in\mathcal{X}}\left\\{\langle h\nabla
f(x_{k})-\sqrt{2h}\sqrt{\nabla^{2}\phi(x_{k})}\,z_{k},\,x-x_{k}\rangle+D_{\phi}(x,x_{k})\right\\}$

where
$D_{\phi}(x,x^{\prime})=\phi(x)-\phi(x^{\prime})-\langle\nabla\phi(x^{\prime}),x-x^{\prime}\rangle$
is the Bregman divergence of $\phi$. In particular, in the Euclidean case,
i.e. when $\mathcal{X}=\mathbb{R}^{d}$ and $\phi(x)=\frac{1}{2}\|x\|^{2}$, MLA
recovers the usual Unadjusted Langevin Algorithm (ULA).

MLA can be seen as a coordinate-transformed Euler-Maruyama discretization of
the Mirror Langevin Dynamics in continuous time, given by

$\begin{cases}Y_{t}&=\nabla\phi(X_{t})\\\ dY_{t}&=-\nabla
f(X_{t})dt+\sqrt{2}\sqrt{\nabla^{2}\phi(X_{t})}dW_{t}\end{cases}.$

See Section 2.3 for a reformulation purely in the dual space, and Appendix A
for more details on the continuous-time dynamics. Zhang et al. (2020) studied
MLA as a simple discretization of the Mirror Langevin Dynamics, and showed
that under certain assumptions, the iterates of MLA converge to a Wasserstein
ball around the target with some radius which depends on the modified self-
concordance parameter of $\phi$ (see Section 3.1 for more detail). In the
Euclidean case (when $\phi$ is quadratic) this radius is $0$, so MLA converges
to $\nu$ with sufficiently small step size, recovering the typical bound for
ULA. However, for general $\phi$, the radius is positive. Therefore, the
result of Zhang et al. (2020) only guarantees MLA enters a ball around the
target, but it may not converge to the target even when the step size goes to
$0$. They further conjectured the bias is unavoidable. This raises an
interesting question of whether the non-vanishing bias of MLA is indeed
unavoidable because we are discretizing a diffusion process with changing
covariance, or whether there is a better analysis of MLA with vanishing bias.
Here we show that indeed MLA has a vanishing bias, by applying the mean-square
analysis framework of Li et al. (2019) and Li et al. (2021).

###  2.3 Mirror Langevin Algorithm in the Dual Space

Let us work in the dual space
$\mathcal{Y}=\nabla\phi(\mathcal{X})=\mathbb{R}^{d}$ via the mirror map
$\nabla\phi\colon\mathcal{X}\to\mathbb{R}^{d}$. Given $x\in\mathcal{X}$, we
define the dual variable

$\displaystyle y=\nabla\phi(x)\in\mathbb{R}^{d}$ (2)

and its inverse is given by $x=\nabla\phi^{*}(y)$. The target distribution
$\tilde{\nu}$ on the dual space is the pushforward of the original target
$\nu\propto e^{-f}$ under the mirror map: $\tilde{\nu}=(\nabla\phi)_{\\#}\nu$.
If we write the density as $\tilde{\nu}(y)\propto e^{-\tilde{f}(y)}$, then we
have $\tilde{f}(y)=f(\nabla\phi^{*}(y))-\log\det\nabla^{2}\phi^{*}(y)$.
Moreover, the Hessian metric $\nabla^{2}\phi(x)$ on $\mathcal{X}$ corresponds
to the Hessian metric $\nabla^{2}\phi^{*}(y)$ on $\mathbb{R}^{d}$ generated by
the dual function $\phi^{*}$; that is, $\nabla^{2}\phi^{*}$ on
$\mathbb{R}^{d}$ is the pullback metric of $\nabla^{2}\phi$ on $\mathcal{X}$
under the inverse mirror map
$\nabla\phi^{*}\colon\mathbb{R}^{d}\to\mathcal{X}$. Therefore, the metric
space $(\mathcal{X},\nabla^{2}\phi)$ is isometric to
$(\mathbb{R}^{d},\nabla^{2}\phi^{*})$.

If $x_{k}\in\mathcal{X}$ follows the Mirror Langevin Algorithm (1), then
$y_{k}=\nabla\phi(x_{k})\in\mathbb{R}^{d}$ follows the Mirror Langevin
Algorithm in the dual space:

$\displaystyle y_{k+1}=y_{k}-h\nabla
f(\nabla\phi^{*}(y_{k}))+\sqrt{2h}\sqrt{\nabla^{2}\phi^{*}(y_{k})^{-1}}\,z_{k}.$
(3)

MLA in the dual space (3) can be seen as a discretization of the mirror
Langevin dynamics to sample from $\tilde{\nu}\propto e^{-\tilde{f}}$ with the
Hessian metric $\nabla^{2}\phi^{*}$ on $\mathbb{R}^{d}$.

Let us define $g\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ and
$A\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$ by

$\displaystyle g(y)$ $\displaystyle=\nabla f(\nabla\phi^{*}(y))$ (4)
$\displaystyle A(y)$ $\displaystyle=\sqrt{\nabla^{2}\phi^{*}(y)^{-1}}.$ (5)

Note here $A(y)$ is any square-root of $\nabla^{2}\phi^{*}(y)^{-1}$. Then we
can write MLA in the dual space as

$\displaystyle y_{k+1}=y_{k}-hg(y_{k})+\sqrt{2h}A(y_{k})\,z_{k}.$ (6)

As $h\to 0$, MLA converges to the Mirror Langevin Dynamics, which is a
continuous-time stochastic process $Y_{t}\in\mathbb{R}^{d}$ following the
stochastic differential equation:

$\displaystyle dY_{t}=-g(Y_{t})dt+\sqrt{2}A(Y_{t})dW_{t}$

where $W_{t}$ is the standard Brownian motion in $\mathbb{R}^{d}$; see Section
4.2.1 for more properties.

###  2.4 Wasserstein distance in dual space

Along MLA in the dual space (3), let $\tilde{\rho}_{k}$ denote the
distribution of the random variable $y_{k}\in\mathbb{R}^{d}$. We will show a
convergence analysis of MLA in the dual space in terms of the Euclidean
Wasserstein distance $W_{2}$ between $\tilde{\rho}_{k}$ and $\tilde{\nu}$ on
$\mathbb{R}^{d}$:

$\displaystyle
W_{2}(\tilde{\rho},\tilde{\nu})^{2}=\inf_{y\sim\tilde{\rho},y^{*}\sim\tilde{\nu}}\mathbb{E}[\|y-y^{*}\|^{2}].$
(7)

Note that this distance does not use the Hessian metric $\nabla^{2}\phi^{*}$
on $\mathbb{R}^{d}$. In the original space $\mathcal{X}$, this gives a
modified $W_{2}$ distance under the mirror map:

$\displaystyle
W_{2,\phi}(\rho,\nu)^{2}=\inf_{x\sim\rho,x^{\prime}\sim\nu}\mathbb{E}[\|\nabla\phi(x)-\nabla\phi(x^{\prime})\|^{2}].$
(8)

That is, if $\tilde{\rho}=(\nabla\phi)_{\\#}\rho$ and
$\tilde{\nu}=(\nabla\phi)_{\\#}\nu$, then
$W_{2,\phi}(\rho,\nu)=W_{2}(\tilde{\rho},\tilde{\nu})$. This is the same
modified Wasserstein distance that is used in Zhang et al. (2020). This
corresponds to using the squared Hessian metric $(\nabla^{2}\phi(x))^{2}$ on
$\mathcal{X}$, which is isometric to the Euclidean metric $I$ on
$\mathbb{R}^{d}$ (rather than the Hessian metric $\nabla^{2}\phi(x)$ on
$\mathcal{X}$, which is isometric to the Hessian metric
$\nabla^{2}\phi^{*}(y)$ on $\mathbb{R}^{d}$, and which is used in the
continuous-time analysis in Chewi et al. (2020)).

##  3 Main Result: Mixing Time Bound for MLA

We present our main result on the mixing time bound of MLA. We need the
following assumptions.

  1. (A1)

$\phi$ satisfies the modified self-concordance property with parameter
$\alpha>0$, which means:

$\displaystyle\|\sqrt{\nabla^{2}\phi(x^{\prime})}-\sqrt{\nabla^{2}\phi(x)}\|_{\mathrm{HS}}\leq\sqrt{\alpha}\|\nabla\phi(x^{\prime})-\nabla\phi(x)\|_{2}~{}~{}~{}~{}~{}\forall\,x^{\prime},x\in\mathcal{X}.$
(9)

Equivalently, $A(y)=\sqrt{\nabla^{2}\phi^{*}(y)^{-1}}$ is
$\sqrt{\alpha}$-Lipschitz in the Hilbert-Schmidt norm:

$\displaystyle\|A(y^{\prime})-A(y)\|_{\mathrm{HS}}\leq\sqrt{\alpha}\|y^{\prime}-y\|_{2}~{}~{}~{}~{}~{}\forall\,y^{\prime},y\in\mathbb{R}^{d}.$
(10)

  2. (A2)

$f$ is $M$-smooth with respect to $\phi$ for some $0<M<\infty$, which means:

$\displaystyle\|\nabla f(x^{\prime})-\nabla f(x)\|_{2}\leq
M\|\nabla\phi(x^{\prime})-\nabla\phi(x)\|_{2}~{}~{}~{}~{}~{}~{}\forall\,x^{\prime},x\in\mathcal{X}.$
(11)

Equivalently, $g(y)=\nabla f(\nabla\phi^{*}(y))$ is $M$-Lipschitz:

$\displaystyle\|g(y^{\prime})-g(y)\|_{2}\leq
M\|y^{\prime}-y\|_{2}~{}~{}~{}~{}~{}~{}\forall\,y^{\prime},y\in\mathcal{Y}.$
(12)

  3. (A3)

$f$ is $m$-strongly convex with respect to $\phi$ for some $0<m\leq M$, which
means:

$\displaystyle\langle\nabla f(x^{\prime})-\nabla
f(x),\nabla\phi(x^{\prime})-\nabla\phi(x)\rangle\geq
m\|\nabla\phi(x^{\prime})-\nabla\phi(x)\|^{2}_{2}~{}~{}~{}~{}~{}\forall\,x^{\prime},x\in\mathcal{X}.$
(13)

Equivalently, $g(y)=\nabla f(\nabla\phi^{*}(y))$ is $m$-monotone:

$\displaystyle\langle g(y^{\prime})-g(y),y^{\prime}-y\rangle\geq
m\|y^{\prime}-y\|^{2}_{2}~{}~{}~{}~{}~{}\forall\,y^{\prime},y\in\mathbb{R}^{d}.$
(14)

These are a subset of the assumptions in Zhang et al. (2020). In particular,
we do not assume a bound on the commutator of $\nabla^{2}f$ and
$\nabla^{2}\phi$. Our main result is the following.

######  Theorem 3.1.

Assume (A1), (A2), (A3), and assume $\alpha<m$. There is a maximum step size
$h_{\max}=\mathcal{O}\left(\frac{(m-\alpha)^{2}}{M^{2}(1+4\alpha)^{2}}\right)$
and constant
$C_{\mathrm{MLA}}=\mathcal{O}\left(\frac{M(1+4\alpha)\sqrt{d}}{m-\alpha}\right)$,
such that if we run MLA (1) with $0<h\leq h_{\max}$ from any
$x_{0}\sim\rho_{0}$, then the iterates $x_{k}\sim\rho_{k}$ satisfy:

$\displaystyle
W_{2,\phi}(\rho_{k},\nu)\leq\sqrt{2}e^{-(m-\alpha)hk}W_{2,\phi}(\rho_{0},\nu)+C_{\mathrm{MLA}}\sqrt{2h}.$
(15)

Equivalently, if we run MLA in the dual space (3) with $0<h\leq h_{1}$ from
any $y_{0}\sim\tilde{\rho}_{0}$, then the iterates $y_{k}\sim\tilde{\rho}_{k}$
satisfy:

$\displaystyle
W_{2}(\tilde{\rho}_{k},\tilde{\nu})\leq\sqrt{2}e^{-(m-\alpha)hk}W_{2}(\tilde{\rho}_{0},\tilde{\nu})+C_{\mathrm{MLA}}\sqrt{2h}.$
(16)

See Section 4.3 for the proof of Theorem 3.1 and explicit forms of the
constant $C_{\mathrm{MLA}}$ and maximum step size $h_{\max}$. This result
shows MLA has a bias that is vanishing with step size, and thus we can reach
an arbitrary accuracy by using a small enough step size. In particular, this
improves on the analysis in Zhang et al. (2020), which has a non-vanishing
bias and under stronger assumptions. By choosing a small step size, we obtain
the following mixing time bound for MLA.

######  Corollary 3.2.

For any (small) error threshold $\epsilon>0$, to reach
$W_{2}(\tilde{\rho}_{k},\tilde{\nu})\leq\epsilon$, it suffices to run MLA in
the dual space (3) with step size
$h=\frac{\epsilon^{2}}{4C_{\mathrm{MLA}}^{2}}$ for $k=\tau_{W_{2}}(\epsilon)$
iterations where

$\displaystyle\tau_{W_{2}}(\epsilon)$
$\displaystyle\leq\frac{1}{(m-\alpha)h}\log\frac{2\sqrt{2}W_{2}(\tilde{\rho}_{0},\tilde{\nu})}{\epsilon}=\tilde{O}\left(\frac{C_{\mathrm{MLA}}^{2}}{(m-\alpha)\epsilon^{2}}\right)=\tilde{O}\left(\frac{M^{2}(1+4\alpha)^{2}d}{(m-\alpha)^{3}\epsilon^{2}}\right).$

###  3.1 Discussion of result

Theorem 3.1 shows that MLA has a biased convergence guarantee where the bias
scales as $O(\sqrt{dh})$ where $d$ is dimension and $h$ is step size (assuming
$m,M,\alpha$ are independent of $d$ for now). This leads to a mixing time
bound of $\tilde{O}(d/\epsilon^{2})$ for MLA.

Let us compare MLA with ULA (i.e., MLA in the Euclidean case with
$\phi(x)=\frac{1}{2}\|x\|^{2}$). Recall for ULA, the mean-square analysis by
Li et al. (2021) yields a biased convergence guarantee where the bias scales
as $O(\sqrt{d}h)$ under an additional 3rd-order regularity condition on $f$.
This leads to a mixing time bound of $\tilde{O}(\sqrt{d}/\epsilon)$ for ULA.
We see the bias of MLA has a worse dependence on $h$ than the bias of ULA.
This is because the continuous-time Mirror Langevin Dynamics (26) of MLA has a
changing covariance, while the usual continuous-time Langevin Dynamics of ULA
has a constant covariance; therefore, MLA incurs an additional stochastic
error from the Brownian motion part, which is not incurred by ULA. Formally,
this is reflected in the orders of error of the two algorithms: We show below
that MLA has local weak and strong errors of orders $p_{1}=\frac{3}{2}$ at
least and $p_{2}=1$ (note the local weak order of MLA is actually $p_{1}=2$,
because it is the Euler-Maruyama discretization of an SDE; the multiplicative
noise causes the strong error to lose half an order, but not the weak error
(see e.g., (Milstein and Tretyakov, 2013, page 14)); however, we will see that
as long as $p_{1}\geq p_{2}+\frac{1}{2}$, the order of the final sampling
error is determined by $p_{2}$ but not $p_{1}$, and even though our
$p_{1}=\frac{3}{2}$ bound is not tight in order, its constants can be made
very explicit and hence helpful to later analysis). On the other hand, it is
well known that ULA has local weak and strong error of orders $p_{1}=2$ and
$p_{2}=\frac{3}{2}$ because it is the Euler-Maruyama discretization of an SDE
with additive noise (see Milstein and Tretyakov (2013) for the general theory
and Li et al. (2021) for details of worked out constants). It would be
interesting to understand whether we can improve the local errors and the bias
of MLA, perhaps using more sophisticated discretization of MLD to improve the
stochastic error.

Our result improves on the analysis of Zhang et al. (2020), who assume
stronger assumptions (our assumptions (A1), (A2), (A3), along with two
assumptions on the moment of $\nabla^{2}\phi$ and a bound on the commutator of
$\nabla^{2}f$ and $\nabla^{2}\phi$), and prove a biased convergence analysis
where the bias scales as $O(\sqrt{dh}+r_{0})$, where $r_{0}=O(\sqrt{\alpha
d})$ does not depend on $h$. Note in the Euclidean case (when
$\phi(x)=\frac{1}{2}\|x\|^{2}$), the modified self-concordance parameter is
$\alpha=0$, and thus $r_{0}=0$; but for general $\phi$, the asymptotic radius
is positive: $r_{0}>0$, so the result of Zhang et al. (2020) does not
guarantee convergence to the target. With our mean-square analysis, we have
shown that in fact there is no dependence on this radius $r_{0}$, and the bias
indeed scales as $O(\sqrt{dh})$.

We note our result uses the modified self-concordance property, as also used
in Zhang et al. (2020). In one-dimension ($d=1$), modified self-concordance is
equivalent to the classical self-concordance property: Both are equivalent to
the condition that $x\mapsto 1/\sqrt{\phi^{\prime\prime}(x)}$ is a Lipschitz
function. However, in higher dimension, they are different. In particular,
modified self-concordance is not an affine-invariant property (in contrast to
the classical self-concordance), and the parameter $\alpha$ can be arbitrarily
large; see example in Appendix D. This is problematic since our convergence
bound only holds when $\alpha$ is less than $m$ (the strong convexity
parameter). It would be desirable to have an analysis of MLA with the more
natural self-concordance property.

Our result in Theorem 3.1 shows that to obtain a consistent algorithm (with a
vanishing bias) from MLD, it suffices to apply a simple discretization such as
MLA. This shows we do not need to use an exact simulator of the Brownian
motion with changing covariance, as proposed by Ahn and Chewi (2020), which
allows a nice analysis under self-concordance property. It would be
interesting to bridge the analysis technique to MLA.

The relative smoothness (A2) and relative strong convexity (A3) conditions
imply that the Hessian of $f$ are bounded by the Hessian of $\phi$:

$\displaystyle m\nabla^{2}\phi(x)\preceq\nabla^{2}f(x)\preceq
M\nabla^{2}\phi(x)~{}~{}~{}~{}~{}~{}\forall~{}x\in\mathcal{X}.$

See (Zhang et al., 2020, Appendix B) for more details. Since we assume $\phi$
is a Legendre function, $\nabla^{2}\phi(x)\to\infty$ as
$x\to\partial\mathcal{X}$; then for our result to hold, we need
$\nabla^{2}f\to\infty$ as $x\to\partial\mathcal{X}$. This restricts the
applicability of the result; for example, it does not apply when $\nu$ is a
uniform ($f=0$) or Gaussian distribution ($f$ is quadratic) restricted on a
polytope with $\phi$ being the log-barrier function. It is desirable to have a
more general convergence analysis of MLA under weaker conditions on $f$ and
$\phi$.

##  4 Proof of main result

The proof of Theorem 3.1 uses the mean-square analysis framework described in
Li et al. (2021). We review the mean-square analysis framework in Section 4.1.
We verify the conditions hold for MLA in Section 4.2, and apply the mean-
square analysis to prove Theorem 3.1 in Section 4.3.

###  4.1 A review of the mean-square analysis framework

Mean-square analysis was a classical tool for analyzing the integration error
of SDEs (e.g., Milstein and Tretyakov (2013)). Li et al. (2019) extended it to
obtain non-asymptotic sampling error bound of an algorithm which is a
discretization of a decaying stochastic differential equation (SDE). While Li
et al. (2019) required the local errors to satisfy uniform bounds, Li et al.
(2021) relaxes this requirement and only needs non-uniform bounds. We will
establish non-uniform local error bounds for MLA, and thus use the version of
mean-square analysis in Li et al. (2021). The results will be reviewed in a
simplified setting; see (Li et al., 2021, Section 3) for details.

##### Contractive SDE.

Consider a continuous-time process $Y_{t}\in\mathbb{R}^{d}$ which evolves
following the SDE:

$\displaystyle dY_{t}=-g(Y_{t})\,dt+\sqrt{2}A(Y_{t})\,dW_{t}$ (17)

for some vector field $g\colon\mathbb{R}^{d}\to\mathbb{R}^{d}$ and matrix
$A\colon\mathbb{R}^{d}\to\mathbb{R}^{d\times d}$. We assume $g$ and $A$ are
Lipschitz continuous. Here $W_{t}$ is the standard Brownian motion in
$\mathbb{R}^{d}$.

We say the SDE (17) is contractive with rate $\beta>0$ if there exists
$t_{0}>0$ such that any two solutions $Y_{t},Y_{t}^{\prime}$ with synchronous
coupling (i.e. driven by the same Brownian motion) satisfy:

$\displaystyle\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}]\leq e^{-2\beta
t}\mathbb{E}[\|Y_{0}-Y_{0}^{\prime}\|^{2}]~{}~{}~{}~{}~{}~{}\forall~{}t\in(0,t_{0}).$
(18)

If the SDE (17) is contractive, then it has a stationary distribution
$\tilde{\nu}$.

##### Short-time deviation.

Since $g$ and $A$ are Lipschitz continuous, one can show (Milstein and
Tretyakov, 2013, Lemma 1.3) that there exist a maximum time $t_{0}>0$ and a
constant $C_{0}>0$ such that for any solutions $Y_{t},Y_{t}^{\prime}$ with
synchronous coupling:

$\displaystyle\mathbb{E}[\|(Y_{t}^{\prime}-Y_{0}^{\prime})-(Y_{t}-Y_{0})\|^{2}_{2}]\leq
C_{0}\,\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]\,t~{}~{}~{}~{}~{}~{}\forall~{}0<t\leq
t_{0}.$ (19)

##### Algorithm and local error.

Suppose we have an algorithm $\mathsf{Alg}_{h}$ depending on a step size $h>0$
that simulates the solution $Y_{t}$ of the SDE (17) at time $t=h$.

For any $Y_{0}\in\mathbb{R}^{d}$, let $Y_{h}$ denote the solution of the SDE
(17) at time $t=h$, and let $\bar{Y}_{1}=\mathsf{Alg}_{h}(Y_{0})$ denote the
output of the algorithm from $Y_{0}$. We say that the algorithm has (non-
uniform) local weak error of order $p_{1}$ if there exist a maximum step size
$h_{1}>0$ and constants $C_{1},D_{1}\geq 0$ such that

$\displaystyle\|\mathbb{E}[Y_{h}-\bar{Y}_{1}]\|\leq\left(C_{1}+D_{1}\sqrt{\mathbb{E}[\|Y_{0}\|^{2}}\right)h^{p_{1}}~{}~{}~{}~{}~{}~{}\forall~{}0<h\leq
h_{1}.$ (20)

We say the algorithm has (non-uniform) local strong error of order $p_{2}$ if
there exist a maximum step size $h_{2}>0$ and constants $C_{2},D_{2}\geq 0$
such that

$\displaystyle\mathbb{E}[\|Y_{h}-\bar{Y}_{1}\|^{2}]\leq\left(C_{2}^{2}+D_{2}^{2}\mathbb{E}[\|Y_{0}\|^{2}\right)h^{2p_{2}}~{}~{}~{}~{}~{}~{}\forall~{}0<h\leq
h_{2}.$ (21)

Here $Y_{h}$ and $\bar{Y}_{1}$ are coupled by sharing the same filtration
(i.e. the algorithm $\mathsf{Alg}_{h}$ has access to the realization of the
Wiener process that generates $Y_{h}$).

When $D_{1}=D_{2}=0$, the bounds are termed as uniform bounds in Li et al.
(2019).

##### Bound on global error.

With the set-up above, the mean-square analysis framework produces the
following bound on the global (long-term) error.

######  Theorem 4.1 ((Li et al., 2021, Theorem 3.3, 3.4)).

Assume the SDE (17) is contractive with rate $\beta>0$. Assume the algorithm
$\mathsf{Alg}_{h}$ has local weak error of order $p_{1}$ and local strong
error of order $p_{2}$ with $\frac{1}{2}<p_{2}\leq p_{1}-\frac{1}{2}$. Let us
define a maximum step size $h_{\max}>0$ by

$\displaystyle
h_{\max}=\min\left\\{t_{0},h_{1},h_{2},\frac{1}{4\beta},\left(\frac{\sqrt{\beta}}{4\sqrt{2}D_{2}}\right)^{\frac{1}{p_{2}-\frac{1}{2}}},\,\left(\frac{\beta}{8\sqrt{2}(D_{1}+C_{0}D_{2})}\right)^{\frac{1}{p_{2}-\frac{1}{2}}}\right\\}$
(22)

and constants
$U=\sqrt{4\mathbb{E}[\|Y_{0}\|^{2}]+6\mathbb{E}_{\tilde{\nu}}[\|Y\|^{2}]}$ and
$C>0$ by

$\displaystyle
C=\frac{2}{\sqrt{\beta}}\left(\frac{C_{1}+C_{0}C_{2}+\sqrt{2}U(D_{1}+C_{0}D_{2})}{\sqrt{\beta}}+C_{2}+\sqrt{2}D_{2}U\right).$
(23)

Starting from any $Y_{0}=\bar{Y}_{0}\sim\tilde{\rho}_{0}$, suppose we run the
algorithm $\mathsf{Alg}_{h}$ with step size $0<h\leq h_{\max}$ to produce
iterates $\bar{Y}_{k}=\mathsf{Alg}_{h}(\bar{Y}_{k-1})\sim\tilde{\rho}_{k}$.
Let $Y_{hk}$ denote the solution to the SDE (17) at time $t=hk$. Then
$\bar{Y}_{k}$ is close to $Y_{hk}$ at all time:

$\displaystyle\sqrt{\mathbb{E}[\|Y_{hk}-\bar{Y}_{k}\|^{2}]}\leq
Ch^{p_{2}-\frac{1}{2}}~{}~{}~{}~{}~{}~{}\forall~{}k\geq 0.$ (24)

Furthermore, the distribution of $\bar{Y}_{k}\sim\tilde{\rho}_{k}$ has the
following biased convergence guarantee:

$\displaystyle W_{2}(\tilde{\rho}_{k},\tilde{\nu})\leq\sqrt{2}e^{-\beta
kh}W_{2}(\tilde{\rho}_{0},\tilde{\nu})+\sqrt{2}Ch^{p_{2}-\frac{1}{2}}~{}~{}~{}~{}~{}~{}\forall~{}k\geq
0.$ (25)

###  4.2 Application to MLA

For our sampling problem, we wish to apply the mean-square analysis framework
to the Mirror Langevin Algorithm in the dual space (3). The continuous-time
SDE (17) of MLA is the Mirror Langevin Dynamics, which we review in the next
section. We establish the local error orders of MLA in the following section.

####  4.2.1 Mirror Langevin Dynamics

Consider the Mirror Langevin Dynamics (MLD), which is a stochastic process
$Y_{t}\in\mathbb{R}^{d}$ following the SDE:

$\displaystyle dY_{t}=-g(Y_{t})dt+\sqrt{2}A(Y_{t})dW_{t}$ (26)

where as defined in (4) and (5), $g(y)=\nabla f(\nabla\phi^{\ast}(y))$ and
$A(y)=\sqrt{\nabla^{2}\phi^{\ast}(y)^{-1}}$. The stationary distribution of
MLD (26) is the target distribution in the dual space:
$\tilde{\nu}=(\nabla\phi)_{\\#}\nu$.

By assumptions (A1) and (A2), $g$ and $A$ are Lipschitz continuous. Let us
establish the contractivity and deviation bound on MLD. The proofs are
provided in Appendix B.

######  Lemma 4.2.

Assume (A1) and (A2) with $\alpha<m$. Then MLD (26) is contractive with rate
$\beta=m-\alpha$.

######  Lemma 4.3.

Assume (A1), (A2), and (A3) with $\alpha<m$. Then any two solutions
$Y_{t},Y_{t}^{\prime}$ of MLD (26) with synchronous coupling satisfy

$\displaystyle\mathbb{E}[\|(Y_{t}^{\prime}-Y_{0}^{\prime})-(Y_{t}-Y_{0})\|^{2}]\leq
4M\,\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}]\,t~{}~{}~{}~{}~{}~{}\forall~{}t\geq
0.$ (27)

We also need the following bound on MLD. Let
$x^{\ast}=\arg\min_{x\in\mathcal{X}}f(x)$ and
$y^{\ast}=\nabla\phi(x^{*})\in\mathbb{R}^{d}$.

######  Lemma 4.4.

Assume (A1), (A2), and (A3). Along MLD (26), for
$0<t\leq\frac{1}{M^{2}+4\alpha}$,

$\displaystyle\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]\leq\gamma\,t$ (28)

where
$\gamma=8(1+4\alpha)\mathbb{E}[\|Y_{0}\|^{2}]+8(1+4\alpha)\|y^{\ast}\|^{2}+16\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{4}{M^{2}}\|g(y^{\ast})\|^{2}$.

######  Remark 1.

In Lemma 4.2 we show MLD is contracting if $\alpha<m$. In general, a bound on
$\alpha$ (the Lipschitz constant of the covariance) is necessary for an SDE
with multiplicative noise to contract; see the example of the geometric
Brownian motion in Appendix C.

####  4.2.2 Local Errors of the Mirror Langevin Algorithm

Let us now consider the algorithm $\mathsf{Alg}_{h}$ to be the Mirror Langevin
Algorithm in the dual space (3). We can show MLA has the following local
errors. The proofs are provided in Appendix B.

######  Lemma 4.5.

Assume (A1), (A2), and (A3). Then MLA (3) has local weak error at least of
order $p_{1}=\frac{3}{2}$, with maximum step size
$h_{1}=\frac{1}{M^{2}+4\alpha}$ and constants

$\displaystyle C_{1}$
$\displaystyle=3M\sqrt{1+4\alpha}\,\left(\|y^{\ast}\|+\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{1}{M}\|g(y^{\ast})\|\right)$
$\displaystyle D_{1}$ $\displaystyle=2M\sqrt{1+4\alpha}.$

######  Lemma 4.6.

Assume (A1), (A2), and (A3). Then MLA (3) has local strong error at least of
order $p_{2}=1$, with maximum step size $h_{2}=\frac{1}{M^{2}+4\alpha}$ and
constants

$\displaystyle C_{2}$
$\displaystyle=7(1+4\alpha)\left(\|y^{\ast}\|+\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{1}{M}\|g(y^{\ast})\|\right)$
$\displaystyle D_{2}$ $\displaystyle=5(1+4\alpha).$

###  4.3 Proof of Theorem 3.1: Convergence Rate of MLA

###### Proof of Theorem 3.1.

Assume (A1), (A2), and (A3) with $\alpha<m$. We have verified that MLA
satisfies the conditions in the mean-square analysis framework: In Lemma 4.2
we show MLD is contractive with rate $\beta=m-\alpha$. We derive the deviation
bound in Lemma 4.3 with $C_{0}=4M$. In Lemmas 4.5 and 4.6 we show MLA has
local weak error of order $p_{1}=\frac{3}{2}$ and local strong error of order
$p_{2}=1$, and indeed $p_{2}\leq p_{1}-\frac{1}{2}$.

Then by Theorem 4.1, we can compute the maximum step size:

$\displaystyle h_{\max}$
$\displaystyle=\min\left\\{\frac{1}{M^{2}+4\alpha},\frac{1}{4\beta},\left(\frac{\sqrt{\beta}}{4\sqrt{2}D_{2}}\right)^{\frac{1}{p_{2}-\frac{1}{2}}},\,\left(\frac{\beta}{8\sqrt{2}(D_{1}+C_{0}D_{2})}\right)^{\frac{1}{p_{2}-\frac{1}{2}}}\right\\}$
$\displaystyle=\min\Bigg{\\{}\frac{1}{M^{2}+4\alpha},\frac{1}{4(m-\alpha)},\frac{m-\alpha}{800(1+4\alpha)^{2}},\frac{(m-\alpha)^{2}}{128\left(2M\sqrt{(1+4\alpha)}+20M(1+8\alpha)\right)^{2}}\Bigg{\\}}$
$\displaystyle=\mathcal{O}\left(\frac{(m-\alpha)^{2}}{M^{2}(1+4\alpha)^{2}}\right).$

Recall $\tilde{\nu}=(\nabla\phi)_{\\#}\nu$ is the target distribution of MLD
(26). We can compute the constant

$\displaystyle
U=\sqrt{4\mathbb{E}[\|Y_{0}\|^{2}_{2}]+6\mathbb{E}_{\tilde{\nu}}[\|Y\|^{2}_{2}]}=O(\sqrt{d}).$

Note that
$\|A(y^{*})\|_{\mathrm{HS}}=\sqrt{\operatorname{Tr}(A(y^{*})A(y^{*})^{\top})}=\sqrt{\operatorname{Tr}(\nabla^{2}\phi^{*}(y^{*})^{-1})}=\sqrt{\operatorname{Tr}(\nabla^{2}\phi(x^{*}))}=O(\sqrt{d})$.
Let us define
$V:=\|y^{\ast}\|+\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{1}{M}\|g(y^{\ast})\|=O(\sqrt{d})$.
Then the resulting constant is

$\displaystyle C_{\mathrm{MLA}}$
$\displaystyle=\frac{2}{\beta}\left(C_{1}+C_{0}C_{2}+\sqrt{2}U(D_{1}+C_{0}D_{2})\right)+\frac{2}{\sqrt{\beta}}\left(C_{2}+\sqrt{2}D_{2}U\right)$
$\displaystyle=\frac{2}{m-\alpha}\Bigg{(}3M\sqrt{(1+4\alpha)}V+28M(1+4\alpha)V+\sqrt{2}U\left(2M\sqrt{(1+4\alpha)}+20M(1+4\alpha)\right)\Bigg{)}$
$\displaystyle~{}~{}~{}~{}+\frac{2}{\sqrt{m-\alpha}}\left(7(1+4\alpha)V+5\sqrt{2}(1+4\alpha)U\right)$
$\displaystyle=\mathcal{O}\left(\frac{M(1+4\alpha)\sqrt{d}}{m-\alpha}\right).$

The conclusion of Theorem 3.1 follows from Theorem 4.1. ∎

##  5 Discussion

Our result leaves open many questions, including the following. It would be
interesting to consider a more sophisticated discretization of MLD such that
the mean-square analysis framework will show improved local errors and smaller
bias.

It would be interesting to have a better analysis of MLA under more natural
conditions on $\phi$, such as self-concordance (rather than modified self-
concordance), and under relaxed requirements on $f$ and $\phi$ (e.g. that
allows us to sample from a uniform or Gaussian distribution on a polytope). It
would be desirable to have a convergence analysis of MLA in the Wasserstein
distance generated by the Hessian metric $\nabla^{2}\phi$ rather than the
Euclidean metric, or in other measures such as KL or $\chi^{2}$-divergence.

It would be interesting to understand whether we can discretize the Newton
Langevin Dynamics (which is the case when $\phi=f$ as described in Appendix
A.4 and which is affine-invariant in continuous time) and obtain a discrete-
time algorithm with a convergence guarantee which is also affine-invariant.

It would also be interesting to understand whether we can derive a more
general discrete-time analysis framework that works under a relaxed condition,
e.g. without requiring contraction in continuous time, but only exponential
convergence in function value (which is known for ULA under the log-Sobolev
inequality, see for example Vempala and Wibisono (2019)).

## References

  * Ahn and Chewi (2020) Kwangjun Ahn and Sinho Chewi.  Efficient constrained sampling via the mirror-Langevin algorithm.  _arXiv preprint arXiv:2010.16212_ , 2020.
  * Brascamp and Lieb (1976) Herm Jan Brascamp and Elliott H Lieb.  On extensions of the Brunn-Minkowski and Prékopa-Leindler theorems, including inequalities for log concave functions, and with an application to the diffusion equation.  _Journal of Functional Analysis_ , 22(4):366–389, 1976.  ISSN 0022-1236.
  * Chewi et al. (2020) Sinho Chewi, Thibaut Le Gouic, Chen Lu, Tyler Maunu, Philippe Rigollet, and Austin Stromme.  Exponential ergodicity of mirror-Langevin diffusions.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 19573–19585. Curran Associates, Inc., 2020.
  * Dalalyan (2017) Arnak Dalalyan.  Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent.  In _Proceedings of the 2017 Conference on Learning Theory_ , volume 65 of _Proceedings of Machine Learning Research_ , pages 678–689. PMLR, 07–10 Jul 2017.  URL http://proceedings.mlr.press/v65/dalalyan17a.html.
  * Dalalyan and Karagulyan (2019) Arnak S Dalalyan and Avetik Karagulyan.  User-friendly guarantees for the Langevin Monte Carlo with inaccurate gradient.  _Stochastic Processes and their Applications_ , 2019.
  * Durmus and Moulines (2017) Alain Durmus and Eric Moulines.  Nonasymptotic convergence analysis for the unadjusted Langevin algorithm.  _Annals of Applied Probability_ , 27(3):1551–1587, 2017.
  * Durmus and Moulines (2019) Alain Durmus and Eric Moulines.  High-dimensional bayesian inference via the unadjusted Langevin algorithm.  _Bernoulli_ , 25(4A):2854–2882, 2019.
  * Fathi (2019) Max Fathi.  _Quelques applications du transport optimal en analyse et en probabilités_.  Habilitation à diriger des recherches, Université Paul Sabatier (Toulouse 3), April 2019.
  * Hsieh et al. (2018) Ya-Ping Hsieh, Ali Kavis, Paul Rolland, and Volkan Cevher.  Mirrored Langevin Dynamics.  In _Advances in Neural Information Processing Systems 31: NeurIPS 2018, Montréal, Canada_ , pages 2883–2892, 2018.
  * Jordan et al. (1998) Richard Jordan, David Kinderlehrer, and Felix Otto.  The variational formulation of the Fokker–Planck equation.  _SIAM Journal on Mathematical Analysis_ , 29(1):1–17, January 1998.
  * Kannan and Narayanan (2012) Ravindran Kannan and Hariharan Narayanan.  Random walks on polytopes and an affine interior point method for linear programming.  _Mathematics of Operations Research_ , 37(1):1–20, 2012.  ISSN 0364765X, 15265471.
  * Laddha et al. (2020) Aditi Laddha, Yin Tat Lee, and Santosh Vempala.  Strong self-concordance and sampling.  In _Proceedings of the 52nd Annual ACM SIGACT Symposium on Theory of Computing_ , pages 1212–1222, 2020.
  * Li et al. (2021) Ruilin Li, Hongyuan Zha, and Molei Tao.  Sqrt($d$) dimension dependence of Langevin Monte Carlo.  _arXiv preprint arXiv:2109.03839_ , 2021.
  * Li et al. (2019) Xuechen Li, Yi Wu, Lester Mackey, and Murat A Erdogdu.  Stochastic Runge-Kutta accelerates Langevin Monte Carlo and beyond.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Milstein and Tretyakov (2013) Grigori Noah Milstein and Michael V Tretyakov.  _Stochastic numerics for mathematical physics_.  Springer Science & Business Media, 2013.
  * Rockafellar (1970) R. Tyrrell Rockafellar.  _Convex Analysis_.  Princeton Landmarks in Mathematics and Physics. Princeton University Press, 1970.  ISBN 978-1-4008-7317-3.
  * Vempala and Wibisono (2019) Santosh Vempala and Andre Wibisono.  Rapid convergence of the Unadjusted Langevin Algorithm: Isoperimetry suffices.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc., 2019.
  * Wibisono (2018) Andre Wibisono.  Sampling as optimization in the space of measures: The Langevin dynamics as a composite optimization problem.  In _Conference On Learning Theory, COLT 2018, Stockholm, Sweden, 6-9 July 2018_ , pages 2093–3027, 2018.  URL http://proceedings.mlr.press/v75/wibisono18a.html.
  * Zhang et al. (2020) Kelvin Shuangjian Zhang, Gabriel Peyré, Jalal Fadili, and Marcelo Pereyra.  Wasserstein control of Mirror Langevin Monte Carlo.  In _Conference on Learning Theory, COLT 2020, 9-12 July 2020, Virtual Event [Graz, Austria]_ , volume 125 of _Proceedings of Machine Learning Research_ , pages 3814–3841. PMLR, 2020.  URL http://proceedings.mlr.press/v125/zhang20a.html.

##  Appendix A Riemannian and Mirror Langevin Dynamics in Continuous Time

Consider the problem of sampling from $\nu\propto e^{-f}$ on
$\mathcal{X}\subseteq\mathbb{R}^{d}$ as described in Section 2.

Suppose we endow $\mathcal{X}$ with a Riemannian metric $\mathsf{g}$, which we
write as a positive definite matrix: $\mathsf{g}(x)\succ 0$ for all
$x\in\mathcal{X}$. This means at each point $x\in\mathcal{X}$ we measure local
norm using the metric $\mathsf{g}(x)$:

$\langle u,v\rangle_{x}:=u^{\top}\mathsf{g}(x)v$

for all $u,v$ in the tangent space. We assume $x\mapsto\mathsf{g}(x)$ is
differentiable. Let $M(x)=\mathsf{g}(x)^{-1}$ be the inverse matrix, and let
$\sqrt{M(x)}$ be a square-root of $M(x)$. Let $\nabla\cdot
M(x)\in\mathbb{R}^{d}$ be the divergence of $M$, which is a vector-valued
function whose entries are the divergences of the columns of $M$. We assume
$\mathsf{g}(x)\to\infty$ (equivalently, $M(x)\to 0$) as $x$ approaches the
boundary of $\mathcal{X}$.

###  A.1 Review for optimization

Recall in optimization, the Riemannian gradient flow (RGF) (or natural
gradient flow) for minimizing $f$ using the metric $\mathsf{g}$ is the
solution $X_{t}$ to the differential equation:

$\dot{X}_{t}=\frac{d}{dt}X_{t}=-M(X_{t})\,\nabla f(X_{t}).$

Here we use the inverse metric $M(x)=\mathsf{g}(x)^{-1}$ to turn the
$\ell_{2}$-gradient $\nabla f(x)=(\frac{\partial f(x)}{\partial
x_{1}},\dots,\frac{\partial f(x)}{\partial x_{d}})$ into a gradient tangent
vector $\mathrm{grad}\,f(x)=M(x)\nabla f(x)$ under the Riemannian metric
$\mathsf{g}(x)$. RGF has nice properties when the objective function $f$
satisfies some properties. For example, if $f$ is geodesically strongly convex
(which means $f$ is strongly convex along geodesics generated by the
Riemannian metric $\mathsf{g}$), then RGF is exponentially contracting.
Moreover, if $f$ is gradient dominated with respect to $\mathsf{g}$, then the
function value $f(X_{t})$ converges exponentially fast along RGF.

Consider when the metric $\mathsf{g}(x)$ is given by the Hessian of a convex
Legendre function $\phi$: $\mathsf{g}(x)=\nabla^{2}\phi(x)\succ 0$. Then the
RGF becomes:

$\dot{X}_{t}=-\nabla^{2}\phi(X_{t})^{-1}\,\nabla f(X_{t}).$

In terms of the dual variable $Y_{t}=\nabla\phi(X_{t})$, this becomes the
mirror flow:

$\dot{Y}_{t}=-\nabla f(X_{t})=-\nabla f(\nabla\phi^{*}(Y_{t})).$

Recall by the mirror map $\nabla\phi$, the metric $\nabla^{2}\phi$ on
$\mathcal{X}$ becomes the Hessian metric $\nabla^{2}\phi^{*}$ on
$\mathcal{Y}=\nabla\phi(X)=\mathbb{R}^{d}$. The mirror flow is also the
Riemannian gradient flow for minimizing the pushforward function
$\tilde{f}(y)=f(\nabla\phi^{*}(y))$ under the Hessian metric
$\nabla^{2}\phi^{*}(y)$ (because
$\mathrm{grad}\,\tilde{f}(y)=\nabla^{2}\phi^{*}(y)^{-1}\nabla\tilde{f}(y)=\nabla
f(\nabla\phi^{*}(y))$). Discretizing the mirror flow gives the mirror descent
algorithm in optimization.

###  A.2 Riemannian Langevin Dynamics

The Riemannian Langevin Dynamics (RLD) for sampling from $\nu\propto e^{-f}$
on $\mathcal{X}$ using the metric $\mathsf{g}(x)$ is the solution $X_{t}$ to
the stochastic differential equation:

$\displaystyle dX_{t}=\left(\nabla\cdot M(X_{t})-M(X_{t})\,\nabla
f(X_{t})\right)dt+\sqrt{2}\sqrt{M(X_{t})}\,dW_{t}.$ (29)

Here $W_{t}$ is the standard Brownian motion in $\mathbb{R}^{d}$. Since
$M(x)\to 0$ as $x\to\partial\mathcal{X}$, the process does not leave
$\mathcal{X}$: If $X_{0}\in\mathcal{X}$, then $X_{t}\in\mathcal{X}$ for all
$t>0$.

The additional drift term $\nabla\cdot M(X_{t})$ accounts for the covariance
$M(X_{t})$ in the Brownian motion. The stationary distribution for RLD is
$\nu(x)\propto e^{-f(x)}$ (the density is with respect to the Lebesgue measure
$dx$ on $\mathbb{R}^{d}$). This can be seen, for example, from the following
Fokker-Planck equation.

If $X_{t}\in\mathcal{X}$ follows RLD (29), then its density
$\rho_{t}\colon\mathcal{X}\to\mathbb{R}$ evolves following the partial
differential equation (PDE):

$\displaystyle\frac{\partial\rho_{t}}{\partial
t}=\nabla\cdot\left(\rho_{t}M\nabla\log\frac{\rho_{t}}{\nu}\right).$ (30)

Clearly if $\rho_{t}=\nu$ then the dynamics is stationary. Furthermore, the
PDE above can be interpreted as the gradient flow for minimizing relative
entropy with respect to the Wasserstein metric on the metric space
$(\mathcal{X},\mathsf{g})$.

From the Fokker-Planck equation (30), we can derive how fast the dynamics RLD
approaches the target distribution $\nu$ in various measures.

For example, recall the $\chi^{2}$-divergence of a probability distribution
$\rho$ with respect to $\nu$ is

$\displaystyle\chi^{2}_{\nu}(\rho)=\mathrm{Var}_{\nu}\left(\frac{\rho}{\nu}\right)=\int_{\mathcal{X}}\nu(x)\left(\frac{\rho(x)}{\nu(x)}-1\right)^{2}dx=\int_{\mathcal{X}}\frac{\rho(x)^{2}}{\nu(x)}dx-1.$

Then a standard calculation reveals that the $\chi^{2}$-divergence is
decreasing along RLD (30):

$\displaystyle\frac{d}{dt}\chi^{2}_{\nu}(\rho_{t})=-2G_{\nu}(\rho_{t})$

where

$\displaystyle
G_{\nu}(\rho)=\mathbb{E}_{\nu}\left[\left\|\nabla\left(\frac{\rho}{\nu}\right)\right\|^{2}_{M}\right]=\int_{\mathcal{X}}\nu(x)\left\langle\nabla\left(\frac{\rho(x)}{\nu(x)}\right),M(x)\nabla\left(\frac{\rho(x)}{\nu(x)}\right)\right\rangle\,dx.$

Therefore, if $\nu$ satisfies a Poincaré inequality with respect to
$\mathsf{g}$, which means for any differentiable function
$h\colon\mathcal{X}\to\mathbb{R}$ we have

$\displaystyle\mathrm{Var}_{\nu}(h)\leq
C_{\mathrm{P}}\,\mathbb{E}_{\nu}[\|\nabla h\|^{2}_{M}],$ (31)

then we can conclude RLD converges exponentially fast in
$\chi^{2}$-divergence: $\chi^{2}_{\nu}(\rho_{t})\leq
e^{-\frac{2}{C_{\mathrm{P}}}t}\chi^{2}_{\nu}(\rho_{0})$.

Similarly, recall the relative entropy (or KL divergence) of $\rho$ with
respect to $\nu$ is

$\displaystyle
H_{\nu}(\rho)=\mathbb{E}_{\nu}\left[\frac{\rho}{\nu}\log\frac{\rho}{\nu}\right]=\int_{\mathcal{X}}\rho(x)\log\frac{\rho(x)}{\nu(x)}\,dx.$

Then along RLD (30), KL divergence is decreasing:

$\displaystyle\frac{d}{dt}H_{\nu}(\rho_{t})=-J_{\nu}(\rho_{t})$

where $J_{\nu}(\rho)$ is the relative Fisher information of $\rho$ with
respect to $\nu$ under the metric $\mathsf{g}$:

$\displaystyle
G_{\nu}(\rho)=\mathbb{E}_{\rho}\left[\left\|\nabla\log\frac{\rho}{\nu}\right\|^{2}_{M}\right].$

Therefore, if $\nu$ satisfies a log-Sobolev inequality with respect to
$\mathsf{g}$, which means for any $\rho$ we have

$H_{\nu}(\rho)\leq C_{\mathrm{LSI}}\,J_{\nu}(\rho),$

then we can conclude RLD converges exponentially fast in KL divergence:
$H_{\nu}(\rho_{t})\leq e^{-\frac{t}{C_{\mathrm{LSI}}}}H_{\nu}(\rho_{0})$.

###  A.3 Mirror Langevin Dynamics

Suppose now the metric $\mathsf{g}(x)$ is given by the Hessian of a convex
Legendre function $\phi$: $\mathsf{g}(x)=\nabla^{2}\phi(x)\succ 0$. The
Riemannian Langevin dynamics (29) becomes the following SDE, which is also
studied by Zhang et al. (2020) and Chewi et al. (2020):

$\displaystyle
dX_{t}=\left(\nabla\cdot(\nabla^{2}\phi(X_{t})^{-1})-\nabla^{2}\phi(X_{t})^{-1}\,\nabla
f(X_{t})\right)dt+\sqrt{2\nabla^{2}\phi(X_{t})^{-1}}\,dW_{t}.$ (32)

If $\nu$ satisfies log-Sobolev or Poincaré inequality (which is called mirror
Poincaré inequality in Chewi et al. (2020)), then we can conclude exponential
convergence rate in KL or $\chi^{2}$ divergence along (32).

The SDE (32) requires $\nabla\cdot(\nabla^{2}\phi(x)^{-1})$, which may be
complicated. Consider the dual variable $Y_{t}=\nabla\phi(X_{t})$. By Itô’s
lemma, $Y_{t}$ evolves following the Mirror Langevin Dynamics:

$\displaystyle dY_{t}=-\nabla
f(\nabla\phi^{\ast}(Y_{t}))\,dt+\sqrt{2\nabla^{2}\phi^{\ast}(Y_{t})^{-1}}\,dW_{t}.$

In particular, the drift term simplifies and there is no third derivative
involved. The mirror Langevin dynamics is also the Riemannian Langevin
dynamics (29) for sampling from the pushforward distribution
$\tilde{\nu}=(\nabla\phi)_{\\#}\nu$ using the Hessian metric
$\nabla^{2}\phi^{*}$. Furthermore, the $\chi^{2}$-divergence and KL divergence
are invariant under the mirror map. Therefore, $\nu$ satisfies LSI or Poincaré
inequality with respect to $\nabla^{2}\phi$ if and only if $\tilde{\nu}$ also
satisfies LSI or Poincaré inequality with respect to $\nabla^{2}\phi^{*}$.
Therefore, we get the same convergence guarantee in both primal and dual
spaces.

###  A.4 Newton Langevin Dynamics

A particularly nice choice of $\phi$ is when $\phi=f$. This gives the Newton
Langevin Dynamics, which in the primal space takes the form:

$\displaystyle
dX_{t}=\left(\nabla\cdot(\nabla^{2}f(X_{t})^{-1})-\nabla^{2}f(X_{t})^{-1}\,\nabla
f(X_{t})\right)dt+\sqrt{2}\sqrt{\nabla^{2}f(X_{t})^{-1}}\,dW_{t}.$ (33)

A remarkable property of NLD, as pointed out by Chewi et al. (2020), is that
the Poincaré inequality of $\nu\propto e^{-f}$ with respect to its Hessian
metric $\nabla^{2}f$ is always true with a uniform constant $C_{\mathrm{P}}=1$
for any strictly log-concave distribution $\nu$, by the virtue of the
Brascamp-Lieb inequality (Brascamp and Lieb, 1976). This gives a uniform
exponential convergence rate along NLD in $\chi^{2}$-divergence as well as the
Wasserstein distance with respect to the metric $\nabla^{2}f$; see detailed
exposition and additional consequences in Chewi et al. (2020).

In the dual space, Newton Langevin Dynamics has a simple drift:

$\displaystyle
dY_{t}=-Y_{t}\,dt+\sqrt{2\nabla^{2}f^{\ast}(Y_{t})^{-1}}\,dW_{t}$ (34)

since $\nabla f(\nabla f^{*}(y))=y$. The target distribution of NLD in the
dual space is the pushforward distribution $\tilde{\nu}=(\nabla f)_{\\#}\nu$
where $\nu\propto e^{-f}$. The SDE (34) for sampling from $\tilde{\nu}$ was
also pointed out by Fathi (2019) from the study of Stein’s kernel.

##  Appendix B Proofs of Lemmas

###  B.1 Proof of Lemma 4.2: Contraction of MLD

###### Proof of Lemma 4.2.

Assume (A1) and (A2). We will show MLD (26) is contractive if
$\alpha<\frac{m}{2}$.

Suppose we have two solutions $Y_{t}^{\prime},Y_{t}$ of MLD (26) with the same
Brownian motion:

$\displaystyle dY_{t}^{\prime}$
$\displaystyle=-g(Y_{t}^{\prime})dt+\sqrt{2}A(Y_{t}^{\prime})dW_{t}$
$\displaystyle dY_{t}$ $\displaystyle=-g(Y_{t})dt+\sqrt{2}A(Y_{t})dW_{t}.$

Then the difference satisfies the SDE

$\displaystyle d(Y_{t}^{\prime}-Y_{t})$
$\displaystyle=-(g(Y_{t}^{\prime})-g(Y_{t}))dt+\sqrt{2}(A(Y_{t}^{\prime})-A(Y_{t}))dW_{t}.$
(35)

Recall in general that if $V_{t}\in\mathbb{R}^{d}$ follows a general SDE
$dV_{t}=b(V_{t})dt+G(V_{t})dW_{t}$, then

$\displaystyle\frac{d}{dt}\mathbb{E}[\|V_{t}\|^{2}]=\mathbb{E}[2\langle
b(V_{t}),V_{t}\rangle+\|G(V_{t})\|^{2}_{\mathrm{HS}}].$

Then for the SDE (35) of the difference $V_{t}=Y_{t}^{\prime}-Y_{t}$, and by
applying assumptions (A1) and (A3), we have

$\displaystyle\frac{d}{dt}\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}]$
$\displaystyle=-2\mathbb{E}[\langle
g(Y_{t}^{\prime})-g(Y_{t}),Y_{t}^{\prime}-Y_{t}\rangle]+2\mathbb{E}[\|A(Y_{t}^{\prime})-A(Y_{t})\|^{2}_{\mathrm{HS}}]$
$\displaystyle\leq-2m\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}_{2}]+2\alpha\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}_{2}]$
$\displaystyle=-2(m-\alpha)\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}_{2}].$

We see that we have an exponential contraction if $\alpha<m$:

$\displaystyle\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}]\leq\exp\left(-2(m-\alpha)t\right)\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}]~{}~{}~{}~{}~{}~{}\forall~{}t\geq
0.$ (36)

This shows that MLD (26) is contractive with rate $\beta=m-\alpha$. ∎

###  B.2 Proof of Lemma 4.3: Deviation bound of MLD

###### Proof of Lemma 4.3.

Assume (A1), (A2), and (A3) with $\alpha<m$.

Suppose we have two solutions $Y_{t}^{\prime},Y_{t}$ of MLD (26) with the same
Brownian motion:

$\displaystyle dY_{t}^{\prime}$
$\displaystyle=-g(Y_{t}^{\prime})dt+\sqrt{2}A(Y_{t}^{\prime})dW_{t}$
$\displaystyle dY_{t}$ $\displaystyle=-g(Y_{t})dt+\sqrt{2}A(Y_{t})dW_{t}.$

Consider the shifted variables
$\tilde{Y}_{t}^{\prime}=Y_{t}^{\prime}-Y_{0}^{\prime}$ and
$\tilde{Y}_{t}=Y_{t}-Y_{0}$, which satisfy:

$\displaystyle d\tilde{Y}_{t}^{\prime}$
$\displaystyle=-g(\tilde{Y}_{t}^{\prime}+Y_{0}^{\prime})dt+\sqrt{2}A(\tilde{Y}_{t}^{\prime}+Y_{0}^{\prime})dW_{t}$
$\displaystyle d\tilde{Y}_{t}$
$\displaystyle=-g(\tilde{Y}_{t}+Y_{0})dt+\sqrt{2}A(\tilde{Y}_{t}+Y_{0})dW_{t}.$

Then the difference
$\tilde{Y}_{t}^{\prime}-\tilde{Y}_{t}=(Y_{t}^{\prime}-Y_{0}^{\prime})-(Y_{t}-Y_{0})$
satisfies:

$\displaystyle d(\tilde{Y}_{t}^{\prime}-\tilde{Y}_{t})$
$\displaystyle=-(g(\tilde{Y}_{t}^{\prime}+Y_{0}^{\prime})-g(\tilde{Y}_{t}+Y_{0}))dt+\sqrt{2}(A(\tilde{Y}_{t}^{\prime}+Y_{0}^{\prime})-A(\tilde{Y}_{t}+Y_{0}))dW_{t}.$

By Lemma 4.2, we have the contraction result (36), which implies
$\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}_{2}]\leq\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]$
for all $t\geq 0$. Then by applying (A1) and (A2) and using $\alpha<m\leq M$,
we get

$\displaystyle\frac{d}{dt}$
$\displaystyle\mathbb{E}[\|(Y_{t}^{\prime}-Y_{0}^{\prime})-(Y_{t}-Y_{0})\|^{2}_{2}]$
$\displaystyle=\frac{d}{dt}\mathbb{E}[\|\tilde{Y}_{t}^{\prime}-\tilde{Y}_{t}\|^{2}_{2}]$
$\displaystyle=-2\mathbb{E}[\langle
g(\tilde{Y}_{t}^{\prime}+Y_{0}^{\prime})-g(\tilde{Y}_{t}+Y_{0}),\tilde{Y}_{t}^{\prime}-\tilde{Y}_{t}\rangle]+2\mathbb{E}[\|A(\tilde{Y}_{t}^{\prime}+Y_{0}^{\prime})-A(\tilde{Y}_{t}+Y_{0})\|^{2}_{\mathrm{HS}}]$
$\displaystyle=-2\mathbb{E}[\langle
g(Y_{t}^{\prime})-g(Y_{t}),Y_{t}^{\prime}-Y_{t}-(Y_{0}^{\prime}-Y_{0})\rangle]+2\mathbb{E}[\|A(Y_{t}^{\prime})-A(Y_{t})\|^{2}_{\mathrm{HS}}]$
$\displaystyle\leq 2\mathbb{E}[\langle
g(Y_{t}^{\prime})-g(Y_{t}),Y_{0}^{\prime}-Y_{0}\rangle]+2\alpha\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}_{2}]$
$\displaystyle\leq
2\mathbb{E}[\|g(Y_{t}^{\prime})-g(Y_{t})\|^{2}_{2}]^{1/2}\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]^{\frac{1}{2}}+2\alpha\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]$
$\displaystyle\leq
2M\mathbb{E}[\|Y_{t}^{\prime}-Y_{t}\|^{2}_{2}]^{\frac{1}{2}}\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]^{\frac{1}{2}}+2\alpha\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]$
$\displaystyle\leq(2M+2\alpha)\,\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]$
$\displaystyle\leq 4M\,\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}].$

Integrating, we conclude that for all $t\geq 0$,

$\displaystyle\mathbb{E}[\|(Y_{t}^{\prime}-Y_{0}^{\prime})-(Y_{t}-Y_{0})\|^{2}_{2}]\leq
4M\,\mathbb{E}[\|Y_{0}^{\prime}-Y_{0}\|^{2}_{2}]\,t.$

∎

###  B.3 Proof of Lemma 4.4: Growth bound of MLD

###### Proof of Lemma 4.4.

Assume (A1), (A2), and (A3). Consider the solution $Y_{t}$ of MLD (26)
starting from $Y_{0}$. The centered variable $\tilde{Y}_{t}=Y_{t}-Y_{0}$
follows the SDE

$\displaystyle
d\tilde{Y}_{t}=-g(\tilde{Y}_{t}+Y_{0})dt+\sqrt{2}A(\tilde{Y}_{t}+Y_{0})dW_{t}.$

Then

$\displaystyle\frac{d}{dt}\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}_{2}]=\frac{d}{dt}\mathbb{E}[\|\tilde{Y}_{t}\|^{2}_{2}]$
$\displaystyle=-2\mathbb{E}[\langle
g(\tilde{Y}_{t}+Y_{0}),\tilde{Y}_{t}\rangle]+2\mathbb{E}[\|A(\tilde{Y}_{t}+Y_{0})\|^{2}_{\mathrm{HS}}]$
$\displaystyle=\underbrace{-2\mathbb{E}[\langle
g(Y_{t}),Y_{t}-Y_{0}\rangle]}_{=~{}I}+\underbrace{2\mathbb{E}[\|A(Y_{t})\|^{2}_{\mathrm{HS}}]}_{=~{}II}.$

Let us bound the two terms above. Let $x^{*}=\arg\min_{x\in\mathcal{X}}f(x)$
and $y^{*}=\nabla\phi(x^{*})$.

##### First term:

By (A2) and (A3),

$\displaystyle I$ $\displaystyle=-2\mathbb{E}[\langle
g(Y_{t}),Y_{t}-Y_{0}\rangle]$ $\displaystyle=-2\mathbb{E}[\langle
g(Y_{t})-g(Y_{0}),Y_{t}-Y_{0}\rangle]-2\mathbb{E}[\langle
g(Y_{0}),Y_{t}-Y_{0}\rangle]$ $\displaystyle\leq-2\mathbb{E}[\langle
g(Y_{0}),Y_{t}-Y_{0}\rangle]$ $\displaystyle\leq
2\mathbb{E}[\|g(Y_{0})\|\cdot\|Y_{t}-Y_{0}\|]$
$\displaystyle\leq\frac{1}{M^{2}}\mathbb{E}[\|g(Y_{0})\|^{2}]+M^{2}\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]$
$\displaystyle\leq
2\mathbb{E}[\|Y_{0}-y^{\ast}\|^{2}]+M^{2}\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]+\frac{2}{M^{2}}\|g(y^{\ast})\|^{2}.$

In the last step we have used $\|g(y)\|^{2}_{2}\leq
2\|g(y)-g(y^{\ast})\|^{2}+2\|g(y^{\ast})\|^{2}\leq
2M^{2}\|y-y^{\ast}\|^{2}+2\|g(y^{\ast})\|^{2}$.

##### Second term:

By triangle inequality and (A1),

$\displaystyle\|A(Y_{t})\|^{2}_{\mathrm{HS}}$ $\displaystyle\leq
2\|A(Y_{t})-A(Y_{0})\|^{2}_{\mathrm{HS}}+2\|A(Y_{0})\|^{2}_{\mathrm{HS}}$
$\displaystyle\leq
2\|A(Y_{t})-A(Y_{0})\|^{2}_{\mathrm{HS}}+4\|A(Y_{0})-A(y^{\ast})\|^{2}_{\mathrm{HS}}+4\|A(y^{\ast})\|^{2}_{\mathrm{HS}}$
$\displaystyle\leq
2\alpha\|Y_{t}-Y_{0}\|^{2}_{2}+4\alpha\|Y_{0}-y^{\ast}\|^{2}+4\|A(y^{\ast})\|^{2}_{\mathrm{HS}}.$

Therefore,

$\displaystyle II$ $\displaystyle=2\mathbb{E}[\|A(Y_{t})\|^{2}_{\mathrm{HS}}]$
$\displaystyle\leq
4\alpha\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}_{2}]+8\alpha\mathbb{E}[\|Y_{0}-y^{\ast}\|^{2}]+8\|A(y^{\ast})\|^{2}_{\mathrm{HS}}.$

Combining the two terms above, we get that along MLD (26):

$\displaystyle\frac{d}{dt}\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}_{2}]$
$\displaystyle\leq(M^{2}+4\alpha)\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}_{2}]+D$ (37)

where

$\displaystyle D$
$\displaystyle=(2+8\alpha)\mathbb{E}[\|Y_{0}-y^{\ast}\|^{2}]+8\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{2}{M^{2}}\|g(y^{\ast})\|^{2}$
$\displaystyle\leq
4(1+4\alpha)\mathbb{E}[\|Y_{0}\|^{2}]+4(1+4\alpha)\|y^{\ast}\|^{2}+8\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{2}{M^{2}}\|g(y^{\ast})\|^{2}.$

Recall in general if $V_{t}\geq 0$ satisfies $\frac{d}{dt}V_{t}\leq CV_{t}+D$
for some $C,D>0$, then $V_{t}\leq e^{Ct}V_{0}+\frac{D}{C}(e^{Ct}-1)$.
Furthermore, if $V_{0}=0$ and $0<t\leq\frac{1}{C}$, then
$V_{t}\leq\frac{D}{C}2Ct=2Dt$. Applying this to
$V_{t}=\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]$ which satisfies (37) and $V_{0}=0$, we
conclude that if $0<t\leq\frac{1}{M^{2}+4\alpha}$, then

$\displaystyle\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}_{2}]\leq\gamma t$

where $\gamma=2D\leq
8(1+4\alpha)\mathbb{E}[\|Y_{0}\|^{2}]+8(1+4\alpha)\|y^{\ast}\|^{2}+16\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{4}{M^{2}}\|g(y^{\ast})\|^{2}$.
∎

###  B.4 Proof of Lemma 4.5: Local weak error of MLA

###### Proof of Lemma 4.5.

Assume (A1), (A2), and (A3). Starting from $Y_{0}\in\mathbb{R}^{d}$, let
$Y_{t}$ be the solution to the MLD (26), and let $Y_{t}^{\prime}$ be the
solution to the modified SDE with constant drift, driven by the same Brownian
motion:

$\displaystyle dY_{t}$ $\displaystyle=-g(Y_{t})dt+\sqrt{2}A(Y_{t})dW_{t}$
$\displaystyle dY_{t}^{\prime}$
$\displaystyle=-g(Y_{0})dt+\sqrt{2}A(Y_{0})dW_{t}.$

The value $Y_{h}^{\prime}$ at time $t=h$ is the output $\bar{Y}_{1}$ of MLA
(3) from $Y_{0}$. We wish to bound $\|\mathbb{E}[Y_{h}-Y_{h}^{\prime}]\|$.

Since $Y_{t},Y_{t}^{\prime}$ are coupled using the same Brownian motion, the
difference $Y_{t}-Y_{t}^{\prime}$ satisfies

$\displaystyle d(Y_{t}-Y_{t}^{\prime})$
$\displaystyle=-(g(Y_{t})-g(Y_{0}))dt+\sqrt{2}(A(Y_{t})-A(Y_{0}))dW_{t}.$

Integrating, and since $Y_{0}=Y_{0}^{\prime}$, this means

$\displaystyle
Y_{h}-Y_{h}^{\prime}=-\int_{0}^{h}(g(Y_{t})-g(Y_{0}))dt+\sqrt{2}\int_{0}^{h}(A(Y_{t})-A(Y_{0}))dW_{t}.$

Taking expectation gives

$\displaystyle\mathbb{E}[Y_{h}-Y_{h}^{\prime}]=-\int_{0}^{h}\mathbb{E}[g(Y_{t})-g(Y_{0})]dt.$
(38)

By (A2) and Lemma 4.4, for $0<t\leq\frac{1}{M^{2}+4\alpha}$ we have

$\displaystyle\mathbb{E}[\|g(Y_{t})-g(Y_{0})\|]\leq
M\mathbb{E}[\|Y_{t}-Y_{0}\|]\leq M\sqrt{\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]}\leq
M\sqrt{\gamma t}.$

Therefore, by triangle inequality on (38), for
$0<h\leq\frac{1}{M^{2}+4\alpha}$ we have

$\displaystyle\|\mathbb{E}[Y_{h}-Y_{h}^{\prime}]\|$
$\displaystyle\leq\int_{0}^{h}\mathbb{E}[\|g(Y_{t})-g(Y_{0})\|]dt$
$\displaystyle\leq M\sqrt{\gamma}\int_{0}^{h}\sqrt{t}\,dt$
$\displaystyle=\frac{2}{3}M\sqrt{\gamma}\,h^{\frac{3}{2}}$
$\displaystyle=\frac{2}{3}M\left(8(1+4\alpha)\mathbb{E}[\|Y_{0}\|^{2}]+8(1+4\alpha)\|y^{\ast}\|^{2}+16\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{4}{M^{2}}\|g(y^{\ast})\|^{2}\right)^{\frac{1}{2}}h^{\frac{3}{2}}$
$\displaystyle\leq\frac{2}{3}M\left(\sqrt{8(1+4\alpha)}\sqrt{\mathbb{E}[\|Y_{0}\|^{2}]}+\sqrt{8(1+4\alpha)}\|y^{\ast}\|+4\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{2}{M}\|g(y^{\ast})\|\right)h^{\frac{3}{2}}$
$\displaystyle=\left(C_{1}+D_{1}\sqrt{\mathbb{E}[\|Y_{0}\|^{2}]}\right)h^{3/2}.$

This shows the local weak error order is at least $p_{1}=\frac{3}{2}$, with
maximum step size $h_{1}=\frac{1}{M^{2}+8\alpha}$ and constants

$\displaystyle C_{1}$
$\displaystyle=\frac{2}{3}M\left(\sqrt{8(1+4\alpha)}\|y^{\ast}\|+4\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{2}{M}\|g(y^{\ast})\|\right)$
$\displaystyle\leq
3M\sqrt{(1+4\alpha)}\,\left(\|y^{\ast}\|+\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{1}{M}\|g(y^{\ast})\|\right)$
$\displaystyle D_{1}$ $\displaystyle=\frac{2}{3}M\sqrt{8(1+4\alpha)}$
$\displaystyle\leq 2M\sqrt{1+4\alpha}.$

∎

###  B.5 Proof of Lemma 4.6: Local strong error of MLA

###### Proof of Lemma 4.6.

Assume (A1), (A2), and (A3). As in the proof of Lemma 4.5, consider two
dynamics $Y_{t}^{\prime},Y_{t}$ starting from $Y_{0}^{\prime}=Y_{0}$ following
the SDEs coupled with the same Brownian motion:

$\displaystyle dY_{t}$ $\displaystyle=-g(Y_{t})dt+\sqrt{2}A(Y_{t})dW_{t}$
$\displaystyle dY_{t}^{\prime}$
$\displaystyle=-g(Y_{0})dt+\sqrt{2}A(Y_{0})dW_{t}.$

We wish to bound $\mathbb{E}[\|Y_{h}-Y_{h}^{\prime}\|^{2}]$. The difference
$Y_{t}-Y_{t}^{\prime}$ satisfies

$\displaystyle d(Y_{t}-Y_{t}^{\prime})$
$\displaystyle=-(g(Y_{t})-g(Y_{0}))dt+\sqrt{2}(A(Y_{t})-A(Y_{0}))dW_{t}.$

By (A1), (A2), and Lemma 4.4, for $0<t\leq\frac{1}{M^{2}+4\alpha}$ we have

$\displaystyle\frac{d}{dt}\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}_{2}]$
$\displaystyle=-2\mathbb{E}[\langle
g(Y_{t})-g(Y_{0}),Y_{t}-Y_{t}^{\prime}\rangle]+2\mathbb{E}[\|A(Y_{t})-A(Y_{0})\|^{2}_{\mathrm{HS}}]$
$\displaystyle\leq
2\mathbb{E}[\|g(Y_{t})-g(Y_{0})\|^{2}]^{\frac{1}{2}}\,\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}]^{\frac{1}{2}}+2\alpha\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]$
$\displaystyle\leq
2M\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]^{\frac{1}{2}}\,\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}]^{\frac{1}{2}}+2\alpha\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]$
$\displaystyle\leq
M^{2}\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}]+(1+2\alpha)\mathbb{E}[\|Y_{t}-Y_{0}\|^{2}]$
$\displaystyle\leq
M^{2}\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}]+(1+2\alpha)\gamma t.$

Equivalently,
$\frac{d}{dt}(e^{-M^{2}t}\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}_{2}])\leq
e^{-M^{2}t}(1+2\alpha)\gamma t\leq(1+2\alpha)\gamma t$, so

$\displaystyle\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}_{2}]\leq
e^{M^{2}t}\frac{(1+2\alpha)}{2}\gamma t^{2}.$

Furthermore, since $t\leq\frac{1}{M^{2}+4\alpha}\leq\frac{1}{M^{2}}$, we have
$e^{M^{2}t}\leq e<3$, so

$\displaystyle\mathbb{E}[\|Y_{t}-Y_{t}^{\prime}\|^{2}_{2}]$
$\displaystyle\leq\frac{3}{2}(1+2\alpha)\gamma t^{2}$
$\displaystyle=3(1+2\alpha)\left(8(1+4\alpha)\mathbb{E}[\|Y_{0}\|^{2}]+8(1+4\alpha)\|y^{\ast}\|^{2}+16\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{4}{M^{2}}\|g(y^{\ast})\|^{2}\right)t^{2}$
$\displaystyle=(C_{2}^{2}+D_{2}^{2}\,\mathbb{E}[\|Y_{0}\|^{2}])\,t^{2}.$

This shows the local strong error order is at least $p_{2}=1$ with maximum
step size $h_{2}=\frac{1}{M^{2}+4\alpha}$ and constants

$\displaystyle C_{2}$
$\displaystyle=\left(24(1+2\alpha)(1+4\alpha)\|y^{\ast}\|^{2}+48(1+2\alpha)\|A(y^{\ast})\|^{2}_{\mathrm{HS}}+\frac{12(1+2\alpha)}{M^{2}}\|g(y^{\ast})\|^{2}\right)^{\frac{1}{2}}$
$\displaystyle\leq
5(1+4\alpha)\|y^{\ast}\|+7\sqrt{1+2\alpha}\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{4\sqrt{1+2\alpha}}{M}\|g(y^{\ast})\|$
$\displaystyle\leq
7(1+4\alpha)\left(\|y^{\ast}\|+\|A(y^{\ast})\|_{\mathrm{HS}}+\frac{1}{M}\|g(y^{\ast})\|\right)$
$\displaystyle D_{2}$ $\displaystyle=\sqrt{24(1+2\alpha)(1+4\alpha)}$
$\displaystyle\leq 5(1+4\alpha).$

∎

##  Appendix C An Analogy: Geometric Brownian Motion

We wondered if our requirement on the modified self-concordance parameter
$\alpha$ being upper-bounded is an artifact of our proof technique. Thus we
did some simple calculations on Geometric Brownian Motion (GBM) which is an
SDE with multiplicative noise and yet admitting close-form solution. It is not
an exact example of MLD but only an analogy; nevertheless, GBM does need
$\alpha$ to be bounded in order to converge.

More precisely, consider GBM on $\mathbb{R}_{+}=(0,\infty)$ which follows the
stochastic differential equation:

$\displaystyle dY_{t}=-Y_{t}\,dt+\sqrt{2\alpha}\,Y_{t}\,dW_{t}$ (39)

where $dW_{t}$ is the standard Brownian motion on $\mathbb{R}$. This has exact
solution

$\displaystyle
Y_{t}=Y_{0}\exp\left(-(1+\alpha)t+\sqrt{2\alpha}\,W_{t}\right).$

By a standard calculation, we see there is a threshold $\alpha<1$ for the
convergence of $Y_{t}$ as $t\to\infty$. Recall since
$W_{t}\sim\mathcal{N}(0,t)$, $\mathbb{E}[\exp(\sigma
W_{t})]=e^{\sigma^{2}t/2}$ for all $\sigma>0$. Then

$\displaystyle\mathbb{E}[Y_{t}^{2}]$
$\displaystyle=\mathbb{E}[Y_{0}^{2}]\,e^{-2(1+\alpha)t}\,\mathbb{E}[\exp(2\sqrt{2\alpha}W_{t})]=\mathbb{E}[Y_{0}^{2}]\,e^{-2(1-\alpha)t}.$

Therefore,

$\displaystyle\lim_{t\to\infty}\mathbb{E}[Y_{t}^{2}]=\begin{cases}0~{}~{}~{}&\text{
if }\alpha<1\\\ \mathbb{E}[Y_{0}^{2}]&\text{ if }\alpha=1\\\ \infty&\text{ if
}\alpha>1.\end{cases}$

Now consider a synchronous coupling $Y_{t},\tilde{Y}_{t}$ following GBM (39)
with the same Brownian motion:

$\displaystyle Y_{t}$
$\displaystyle=Y_{0}\exp\left(-(1+\alpha)t+\sqrt{2\alpha}\,W_{t}\right)$
$\displaystyle\tilde{Y}_{t}$
$\displaystyle=\tilde{Y}_{0}\exp\left(-(1+\alpha)t+\sqrt{2\alpha}\,W_{t}\right).$

Then

$\displaystyle\mathbb{E}[(Y_{t}-\tilde{Y}_{t})^{2}]$
$\displaystyle=\mathbb{E}[(Y_{0}-\tilde{Y}_{0})^{2}]\,e^{-2(1+\alpha)t}\,\mathbb{E}[\exp(2\sqrt{2\alpha}W_{t})]=\mathbb{E}[(Y_{0}-\tilde{Y}_{0})^{2}]\,e^{-2(1-\alpha)t}.$

Thus, we see that GBM is a contraction if and only if $\alpha<1$. In
particular, we also have

$\displaystyle\lim_{t\to\infty}\mathbb{E}[(Y_{t}-\tilde{Y}_{t})^{2}]=\begin{cases}0~{}~{}~{}&\text{
if }\alpha<1\\\ \mathbb{E}[(Y_{0}-\tilde{Y}_{0})^{2}]&\text{ if }\alpha=1\\\
\infty&\text{ if }\alpha>1.\end{cases}$

GBM (39) is an instance of MLD (26) (and in fact NLD (33)) with $\phi=f$ where

$\displaystyle\frac{1}{\sqrt{(\phi^{\ast})^{\prime\prime}(y)}}=\sqrt{\alpha}y$
(40)

which is $\sqrt{\alpha}$-Lipschitz, so it satisfies modified self-concordance
(A1) with parameter $\alpha$. Since $\phi=f$, it satisfes relative smoothness
(A2) and relative strong convexity (A3) with $M=m=1$. Note our assumption in
Theorem 3.1 is $\alpha<m=1$, which c is tight for GBM to contract, as well as
to determine if there is a $t\to\infty$ limit.

##  Appendix D Example: Log-Barrier on a Polytope

Let $\mathcal{X}$ be the polytope (not necessarily bounded)

$\mathcal{X}=\\{x\in\mathbb{R}^{d}\colon a_{i}^{\top}x\geq
b_{i}~{}~{}\forall\,i=1,\dots,m\\}$

for some $a_{1},\dots,a_{m}\in\mathbb{R}^{d}$ and
$b_{1},\dots,b_{m}\in\mathbb{R}$. Consider the log-barrier function defined in
the interior of $\mathcal{X}$:

$\displaystyle\phi(x)=-\sum_{i=1}^{m}\log(a_{i}^{\top}x-b_{i}).$ (41)

Recall that $\phi$ satisfies the classical self-concordance condition with a
constant parameter $2$. Let $\alpha$ be the modified self-concordance
parameter of $\phi$. For some polytopes, such as the positive orthant,
$\alpha$ is also a constant (because the Hessian is diagonal and the
dimensions are independent). For general polytopes, however, $\alpha$ can be
arbitrarily large. Here we show $\alpha$ can be as large as the square inverse
of the smallest singular value of the constraint matrix; we also construct an
explicit example in two dimension.

Without loss of generality we may assume $\|a_{i}\|=1$ for $i=1,\dots,m$. Let
$A=(a_{1},\cdots,a_{m})\in\mathbb{R}^{d\times m}$ be the constraint matrix, so
the polytope is described by $A^{\top}x\geq b$. Let the singular values of $A$
be $\sigma_{1}\geq\dots\geq\sigma_{d}\geq 0$ (assuming $d\leq m$). Then
$\sum_{i=1}^{d}\sigma_{i}^{2}=\operatorname{Tr}(AA^{\top})=\operatorname{Tr}(A^{\top}A)=\sum_{i=1}^{m}\|a_{i}\|^{2}=m$;
but $\sigma_{d}=\min_{i}\sigma_{i}$ can be small or $0$. For
$x\in\mathcal{X}$, let $S_{x}\in\mathbb{R}^{m\times m}$ be the diagonal matrix
with entries $a_{i}^{\top}x-b_{i}$.

The gradient of $\phi$ is

$\displaystyle\nabla\phi(x)=-\sum_{i=1}^{m}\frac{a_{i}}{a_{i}^{\top}x-b_{i}}.$

Then for $x,x^{\prime}\in\mathcal{X}$, we have

$\displaystyle\nabla\phi(x^{\prime})-\nabla\phi(x)$
$\displaystyle=\sum_{i=1}^{m}\left(\frac{1}{a_{i}^{\top}x-b_{i}}-\frac{1}{a_{i}^{\top}x^{\prime}-b_{i}}\right)a_{i}$
$\displaystyle=\sum_{i=1}^{m}\frac{a_{i}^{\top}(x^{\prime}-x)}{(a_{i}^{\top}x-b_{i})(a_{i}^{\top}x^{\prime}-b_{i})}a_{i}$
$\displaystyle=\sum_{i=1}^{m}\frac{a_{i}a_{i}^{\top}}{(a_{i}^{\top}x-b_{i})(a_{i}^{\top}x^{\prime}-b_{i})}(x^{\prime}-x)$
$\displaystyle=AS_{x}^{-1}S_{x^{\prime}}^{-1}A^{\top}(x^{\prime}-x).$

Therefore,

$\displaystyle\|\nabla\phi(x^{\prime})-\nabla\phi(x)\|^{2}$
$\displaystyle=\|AS_{x}^{-1}S_{x^{\prime}}^{-1}A^{\top}(x^{\prime}-x)\|^{2}$
$\displaystyle=(x^{\prime}-x)^{\top}AS_{x^{\prime}}^{-1}S_{x}^{-1}A^{\top}AS_{x}^{-1}S_{x^{\prime}}^{-1}A^{\top}(x^{\prime}-x)$
$\displaystyle=v(x,x^{\prime})^{\top}A^{\top}Av(x,x^{\prime})$

where

$v(x,x^{\prime})=S_{x}^{-1}S_{x^{\prime}}^{-1}A^{\top}(x^{\prime}-x)\in\mathbb{R}^{m}.$

The Hessian is

$\displaystyle\nabla^{2}\phi(x)=\sum_{i=1}^{m}\frac{a_{i}a_{i}^{\top}}{(a_{i}^{\top}x-b_{i})^{2}}=AS_{x}^{-2}A^{\top}.$

As a square-root, we can choose:

$\displaystyle\sqrt{\nabla^{2}\phi(x)}=AS_{x}^{-1}=\begin{pmatrix}\dfrac{a_{1}}{a_{1}^{\top}x-b_{1}}&\cdots&\dfrac{a_{m}}{a_{m}^{\top}x-b_{m}}\end{pmatrix}$

since indeed
$\sqrt{\nabla^{2}\phi(x)}\sqrt{\nabla^{2}\phi(x)}^{\top}=AS_{x}^{-1}S_{x}^{-1}A^{\top}=\nabla^{2}\phi(x)$.

For $x,x^{\prime}\in\mathcal{X}$, we have that

$\displaystyle\sqrt{\nabla^{2}\phi(x^{\prime})}-\sqrt{\nabla^{2}\phi(x)}$
$\displaystyle=A(S_{x^{\prime}}-S_{x})$
$\displaystyle=-\begin{pmatrix}\dfrac{a_{1}a_{1}^{\top}(x^{\prime}-x)}{(a_{1}^{\top}x^{\prime}-b_{1})(a_{1}^{\top}x-b_{1})}&\cdots&\dfrac{a_{m}a_{m}^{\top}(x^{\prime}-x)}{(a_{m}^{\top}x^{\prime}-b_{m})(a_{m}^{\top}x-b_{m})}\end{pmatrix}.$

Therefore,

$\displaystyle\|\sqrt{\nabla^{2}\phi(x^{\prime})}-\sqrt{\nabla^{2}\phi(x)}\|^{2}_{\mathrm{HS}}$
$\displaystyle=\sum_{i=1}^{m}\left\|\frac{a_{1}a_{1}^{\top}(x^{\prime}-x)}{(a_{1}^{\top}x^{\prime}-b_{1})(a_{1}^{\top}x-b_{1})}\right\|^{2}$
$\displaystyle=\sum_{i=1}^{m}(x^{\prime}-x)^{\top}\frac{a_{1}a_{1}^{\top}a_{1}a_{1}^{\top}}{(a_{1}^{\top}x^{\prime}-b_{1})^{2}(a_{1}^{\top}x-b_{1})^{2}}(x^{\prime}-x)$
$\displaystyle=(x^{\prime}-x)^{\top}\left(\sum_{i=1}^{m}\frac{a_{1}a_{1}^{\top}}{(a_{1}^{\top}x^{\prime}-b_{1})^{2}(a_{1}^{\top}x-b_{1})^{2}}\right)(x^{\prime}-x)$
$\displaystyle=(x^{\prime}-x)^{\top}AS_{x^{\prime}}^{-2}S_{x}^{-2}A^{\top}(x^{\prime}-x)$
$\displaystyle=\|v(x,x^{\prime})\|^{2}.$

##### Modified self-concordance.

The modified self-concordance parameter is

$\displaystyle\alpha$
$\displaystyle=\sup_{x,x^{\prime}\in\mathcal{X}}\frac{\|\sqrt{\nabla^{2}\phi(x^{\prime})}-\sqrt{\nabla^{2}\phi(x)}\|^{2}_{\mathrm{HS}}}{\|\nabla\phi(x^{\prime})-\nabla\phi(x)\|^{2}_{2}}$
$\displaystyle=\sup_{x,x^{\prime}\in\mathcal{X}}\frac{\|v(x,x^{\prime})\|^{2}}{v(x,x^{\prime})^{\top}A^{\top}Av(x,x^{\prime})}$
$\displaystyle\leq\sup_{v\in\mathbb{R}^{d}}\frac{\|v\|^{2}}{v^{\top}A^{\top}Av}$
$\displaystyle=\max_{i=1,\dots,d}\frac{1}{\sigma_{i}^{2}}$
$\displaystyle=\frac{1}{\sigma_{d}^{2}}.$

This shows the modified self-concordance parameter can be as large as
$\frac{1}{\sigma_{d}^{2}}$, by choosing appropriate $x,x^{\prime}$. For some
polyhedra $\sigma_{d}\approx 0$, so $\alpha\approx 1/\sigma_{d}^{2}$ can be
arbitrarily large.

##### Example in two dimension.

Let $d=2$, and consider

$a_{1}=\begin{pmatrix}1\\\
0\end{pmatrix},~{}~{}~{}~{}~{}~{}a_{2}=\begin{pmatrix}\sqrt{1-\epsilon^{2}}\\\
\epsilon\end{pmatrix}$

for some small $\epsilon>0$, and $b_{1}=b_{2}=0$. This defines the
intersection of two halfspaces:

$\mathcal{X}=\\{x=(x_{1},x_{2})\colon x_{1}\geq
0,~{}\sqrt{1-\epsilon^{2}}\,x_{1}+\epsilon x_{2}\geq 0\\}.$

The constraint matrix is $A=\begin{pmatrix}1&\sqrt{1-\epsilon^{2}}\\\
0&\epsilon\end{pmatrix}.$ We have
$A^{\top}A=\begin{pmatrix}1&\sqrt{1-\epsilon^{2}}\\\
\sqrt{1-\epsilon^{2}}&1\end{pmatrix}$ which has eigenvalues
$\sigma_{1}^{2}=1+\sqrt{1-\epsilon^{2}}$ and
$\sigma_{2}^{2}=1-\sqrt{1-\epsilon^{2}}$. Note that if $\epsilon$ is small,
$\sigma_{1}^{2}\approx 2$ and $\sigma_{2}^{2}\approx\epsilon^{2}/2$. The
corresponding eigenvectors are $v_{1}=\begin{pmatrix}1\\\ 1\end{pmatrix}$ and
$v_{2}=\begin{pmatrix}1\\\ -1\end{pmatrix}$.

Let us choose

$x=\begin{pmatrix}1\\\
0\end{pmatrix},~{}~{}~{}~{}~{}x^{\prime}=\begin{pmatrix}a\\\ b\end{pmatrix}$

for some constant $a,b\in\mathbb{R}$. For simplicity let
$s=\sqrt{1-\epsilon^{2}}$. We require $x^{\prime}\in\mathcal{X}$, so $a\geq 0$
and $b\geq-\frac{s}{\epsilon}a$. We have

$A^{\top}x=\begin{pmatrix}1\\\
s\end{pmatrix},~{}~{}~{}~{}~{}~{}A^{\top}x^{\prime}=\begin{pmatrix}a\\\
sa+\epsilon b\end{pmatrix}$

and

$A^{\top}(x^{\prime}-x)=\begin{pmatrix}a-1\\\ s(a-1)+\epsilon b\end{pmatrix}.$

We also have

$S_{x}=\begin{pmatrix}1&0\\\
0&s\end{pmatrix},~{}~{}~{}~{}~{}S_{x}=\begin{pmatrix}a&0\\\ 0&sa+\epsilon
b\end{pmatrix}.$

Then

$\displaystyle v(x,x^{\prime})$
$\displaystyle=S_{x}^{-1}S_{x^{\prime}}^{-1}A^{\top}(x^{\prime}-x)$
$\displaystyle=\begin{pmatrix}\frac{1}{a}&0\\\ 0&\frac{1}{s(sa+\epsilon
b)}\end{pmatrix}\begin{pmatrix}a-1\\\ s(a-1)+\epsilon b\end{pmatrix}$
$\displaystyle=\begin{pmatrix}\frac{a-1}{a}\\\ \frac{s(a-1)+\epsilon
b}{s(sa+\epsilon b)}\end{pmatrix}$

We want this to be proportional to $v_{2}=\begin{pmatrix}1\\\
-1\end{pmatrix}$, so we want

$\displaystyle\frac{a-1}{a}+\frac{s(a-1)+\epsilon b}{s(sa+\epsilon b)}=0.$
(42)

We can solve for $b$ in terms of $a$:

$b=-\frac{a(a-1)s(s+1)}{\epsilon((a-1)s+a)}.$

We can see that for all $a\geq 0$, this choice of $b$ satisfies the constraint
$b\geq-\frac{s}{\epsilon}a$, so $x^{\prime}\in\mathcal{X}$.

Explicitly, suppose we choose

$\displaystyle a$ $\displaystyle=2$ $\displaystyle b$
$\displaystyle=-\frac{2s(s+1)}{\epsilon(s+2)}$

which satisfies the condition $b\geq-\frac{2s}{\epsilon}$. We can verify
directly that the condition (42) holds:

$\displaystyle\frac{a-1}{a}+\frac{s(a-1)+\epsilon b}{s(sa+\epsilon b)}$
$\displaystyle=\frac{1}{2}+\frac{s-\frac{2s(s+1)}{(s+2)}}{s(2s-\frac{2s(s+1)}{(s+2)})}=\frac{1}{2}+\frac{-\frac{s^{2}}{s+2}}{s(\frac{2s}{(s+2)})}=\frac{1}{2}-\frac{1}{2}=0.$

Then with this choice

$x=\begin{pmatrix}1\\\
0\end{pmatrix},~{}~{}~{}~{}~{}x^{\prime}=\begin{pmatrix}2\\\
-\frac{2s(s+1)}{\epsilon(s+2)}\end{pmatrix}$

we have that $v(x,x^{\prime})=\frac{1}{2}v_{2}$, i.e. proportional to the
eigenvector of $A^{\top}A$ with small eigenvalue $\sigma_{2}^{2}$. Then
$A^{\top}Av(x,x^{\prime})=\sigma_{2}^{2}v(x,x^{\prime})$, and this gives the
bound for the modified self-concordance parameter:

$\displaystyle\alpha\geq\frac{\|v(x,x^{\prime})\|^{2}}{v(x,x^{\prime})^{\top}A^{\top}Av(x,x^{\prime})}=\frac{\|v(x,x^{\prime})\|^{2}}{\sigma_{2}^{2}\|v(x,x^{\prime})\|^{2}}=\frac{1}{\sigma_{2}^{2}}=\frac{1}{1-\sqrt{1-\epsilon^{2}}}\approx\frac{2}{\epsilon^{2}}.$

Thus, by setting $\epsilon\to 0$ we can make $\alpha$ as large as we want.
However, note that the case $\epsilon=0$ is nice and we have $\alpha=1$,
because the domain is a half-space and the problem reduces to one dimension.
This example shows the definition of modified self-concordance is not stable.