# Convergence Error Analysis of Reflected Gradient Langevin Dynamics for
Globally Optimizing Non-Convex Constrained Problems

Kanji Sato  Graduate School of Information Science and Technology, The
University of Tokyo  Akiko Takeda  Graduate School of Information Science and
Technology, The University of Tokyo  Center for Advanced Intelligence Project,
RIKEN  Reiichiro Kawai  Graduate School of Arts and Sciences, The University
of Tokyo  Taiji Suzuki  Graduate School of Information Science and Technology,
The University of Tokyo  Center for Advanced Intelligence Project, RIKEN

###### Abstract

Non-convex optimization problems have various important applications, whereas
many algorithms have been proven only to converge to stationary points.
Meanwhile, gradient Langevin dynamics (GLD) and its variants have attracted
increasing attention as a framework to provide theoretical convergence
guarantees for a global solution in non-convex settings. The studies on GLD
initially treated unconstrained convex problems and very recently expanded to
convex constrained non-convex problems by Lamperski, (2021). In this work, we
can deal with non-convex problems with some kind of non-convex feasible
region. This work analyzes reflected gradient Langevin dynamics (RGLD), a
global optimization algorithm for smoothly constrained problems, including
non-convex constrained ones, and derives a convergence rate to a solution with
$\epsilon$-sampling error. The convergence rate is faster than the one given
by Lamperski, (2021) for convex constrained cases. Our proofs exploit the
Poisson equation to effectively utilize the reflection for the faster
convergence rate.

##  1 Introduction

We are interested in the following (generally) non-convex problems over a
domain $K\subset\mathbb{R}^{d}$ defined by

$\displaystyle\begin{split}&\min_{x\in\mathbb{R}^{d}}f(x)\quad\mathrm{s.t.}\quad
x\in K,\end{split}$ (1)

where $f:{\mathbb{R}}^{d}\rightarrow\mathbb{R}$ is a non-convex function and
$K$ is a (possibly non-convex) bounded region with smooth boundary. Let us
also assume that $K$ are thick enough to contain a Euclidean ball. The
assumptions on $K$ and $f$ will be specified in Section 3.1 in detail.

We consider the following reflected gradient Langevin dynamics (RGLD) as an
optimization algorithm to minimize the non-convex objective on a possibly non-
convex feasible region $K$. The RGLD algorithm is basically the gradient
Langevin dynamics (or unadjusted Langevin algorithm) but includes a reflection
step so that the intermediate solution does not go outside the feasible
region:

$\displaystyle\begin{cases}X^{\prime}_{k+1}=X_{k}-\eta\gradient
f(X_{k})+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\\\
X_{k+1}=\mathcal{R}_{K}\left(X^{\prime}_{k+1}\right),\end{cases}$

where $\mathcal{R}_{K}$ is the reflection operator defined later in (2),
$\xi_{k+1}$ is an i.i.d. Rademacher random vector in $\mathbb{R}^{d}$ (i.e.,
each element of which is i.i.d. a Rademacher random variable), $\eta(>0)$ is
the step size, and $\beta(>0)$ is the inverse temperature parameter. Although
it was proved that the continuous-time limit of $\\{X_{k}\\}$ converges to a
unique stationary distribution (Bubeck et al.,, 2018), i.e., the Gibbs
distribution $\pi(x)\sim\exp(-\beta f(x))$ supported on $K$, which
concentrates around the global minima of $f(x)$ on $K$, the error analysis of
the discrete-time algorithm has not been extensively explored yet. The main
purpose of this work is to give the error analysis including the convergence
rate and the discretization error of the discrete-time RGLD algorithm to
optimize the problem (1).

Non-convex optimization has been essential in many machine learning
applications such as tensor factorization (Signoretto et al.,, 2013; Suzuki et
al.,, 2016), robust classification using non-convex losses (Masnadi-shirazi
and Vasconcelos,, 2009), deep learning (Kingma and Ba,, 2014). While
unconstrained problems have been widely studied in various contexts,
constrained problems also emerge in many important applications such as
portfolio optimization (Barkhagen et al.,, 2019) and sparse estimation with
non-convex regularization (Wen et al.,, 2018; Zhang et al.,, 2018). However,
the existence of constraints sometimes complicates the problems, which leads
to deteriorated computational efficiency. While various algorithms have been
proposed to solve constrained problems (Li and Lin, (2015); Xie and Wright,
(2021); Hinder and Ye, (2018); Lin et al., (2019)), many of them are proposed
to achieve first-order optimality and thus their theoretical analyses are
limited to proving convergence to stationary points since non-convex
optimization is known to be difficult to globally optimize in general.

However, a line of works on Langevin dynamics, originating from stochastic
diffusions, has focused on global optimization of non-convex problems in
limited situations, where the objective function is required to be smooth to
some extent (Geman and Hwang,, 1986; Chiang et al.,, 1987; Gelfand and
Mitter,, 1991). The algorithms appearing in these works are often called
gradient Langevin dynamics (GLD). Nevertheless, constrained non-convex
problems have not been analyzed until very recently because of the theoretical
difficulties of the underlying stochastic differential equations (SDE); even
proving the existence of the solution of the underlying SDE is non-trivial
(see Tanaka, (1979); Lions and Sznitman, (1984)). Very recently, Lamperski,
(2021) analyzed projected gradient Langevin dynamics (PGLD), proposed by
Bubeck et al., (2018), for convex constrained non-convex problems and gave its
convergence rate. However, his work does not cover non-convex constraints and
leaves a gap between the obtained convergence rates in the unconstrained
setting and that in the constrained one (see Table 1).

| objective function
---|---
| convex | non-convex
unconstrained |  Dalalyan, (2017) $\tilde{O}(\epsilon^{-2})$ |  Zou et al., (2021) $\tilde{O}(d^{4}\epsilon^{-2})$
constrained | convex |  Hsieh et al., (2018) $\tilde{O}(d\epsilon^{-2})$ |  Lamperski, (2021) $\tilde{O}(d^{4}\lambda_{*}^{-1}\epsilon^{-4})$
non-convex |  This work $\tilde{O}(d^{3}\lambda_{*}^{-3}\epsilon^{-3})$
Table 1: Iteration complexity of gradient Langevin dynamics to converge to the
Gibbs distribution $\pi$. $\tilde{O}$ is the order ignoring polylogarithmic
factors. $\lambda_{*}$ is the spectral gap.

##### Our Contributions

This paper focuses on enriching the theoretical analysis of constrained
optimization problems, especially those having (generally non-convex) smooth
constraints, by considering the RGLD algorithm.

While Lamperski, (2021) adopted a coupling argument for the 1-Wasserstein
metric with the convexity assumption of the feasible region $K$ and showed
geometric ergodicity, its time discretization error was $O((\eta\log
k)^{1/4})$ which is slower than the unconstrained counterpart $O(\eta^{1/2})$
in terms of the convergence with respect to the 1-Wasserstein metric. This is
mainly because they employed the projected gradient Langevin dynamics (PGLD)
and gave the rough estimate of the discretization error by the convergence in
mean. Although PGLD has already been equipped with a convergence guarantee for
the convex-constrained problems (e.g., Bubeck et al., (2018)), it induces a
suboptimal rate.

To overcome this difficulty, we replace the projection step with a so-called
“reflection” step to keep the trajectory inside the feasible region. Indeed,
this strategy enables us to give a better convergence guarantee by canceling
out a certain error term. Although the reflection step incurs further
theoretical difficulties, we exploit the Poisson equation of the corresponding
continuous-time dynamics as in Leimkuhler et al., (2020), which yields a
simpler proof by adopting some known results on the dynamics. Moreover, we
explicitly specify the parameter dependence of the error bound to derive the
convergence rate, which was not covered by Leimkuhler et al., (2020, Theorem
4.2). This is important because the optimization requires employing a large
inverse temperature $\beta$, and thus the dependence on $\beta$ and $d$ is of
great interest in optimization problems. As shown in Table 1, our proven
convergence rate is sharper than that of Lamperski, (2021) thanks to the
better discretization error while there remains a gap between the best known
rate under unconstrained formulations and that under constrained ones.

Moreover, this is the first work to address non-convex constrained situations,
requiring a certain smoothness assumption on the feasible region, which has
not been imposed so far (e.g., by Lamperski, (2021)).

##### Related Works on GLD

GLD, including its stochastic variant (SGLD), has been studied for a long
time, starting with Geman and Hwang, (1986). It has also been studied in the
context of Bayesian learning, for instance, by Welling and Teh, (2011). Chiang
et al., (1987); Gelfand and Mitter, (1991), and their following works examined
the convergence of GLD to the global minimum of the objective functions,
whereas Zhang et al., (2017) obtained the convergence to local optima.

Taking advantage of the accumulated probabilistic analyses of GLD and SGLD,
Raginsky et al., (2017); Xu et al., (2018); Vempala and Wibisono, (2019)
proved a non-asymptotic convergence rate to global optima by analyzing GLD for
unconstrained non-convex problems. Zou et al., (2021) is the best known work
that provides the fastest convergence rate. A line of studies has extended the
conventional GLD to its variants. Chen et al., (2021) combined SGLD with the
stochastic variance-reduced gradient (SVRG), while Durmus et al., (2018)
worked on a non-smooth objective function by making use of the Moreau
envelope.

Among such varied researches, Bubeck et al., (2018) is considered to be the
first known work that analyzed projected gradient Langevin dynamics (PGLD) for
constrained sampling. They aimed to sample from a log-concave distribution,
which translates to convex optimization in our context. Brosse et al., (2017)
also focused on sampling from a constrained log-concave distribution by
utilizing a proximal operation. Proximal-type algorithms are also considered
in Pereyra, (2016); Salim et al., (2019); Salim and Richtarik, (2020). Hsieh
et al., (2018); Zhang et al., 2020b  proposed mirrored Langevin dynamics as a
way to handle constraints. They used the mirror map to transform the target
distribution into an unconstrained one, to which many existing methods can be
applied. Finding an appropriate mirror map is a non-trivial task in practice,
though. Patterson and Teh, (2013); Wang et al., (2020) can be located in the
context of constrained optimization when the feasible region forms a
Riemannian manifold such as the Stiefel manifold or the Grassmann manifold.

Parpas et al., (2006) firstly proposed a GLD algorithm that converges to the
global optimum for linearly constrained non-convex problems. Very recently,
Lamperski, (2021) proved a non-asymptotic convergence rate of PGLD for convex
constrained non-convex problems. The key components of the theoretical
analysis are reflected stochastic differential equations (RSDE) and their
discretization, which are also essential for our analysis.

Tanaka, (1979); Lions and Sznitman, (1984) discussed the existence and
uniqueness of solutions of reflected diffusions, named the “Skorokhod
problem”. Dupuis and Ramanan, (1999); Andres, (2011) studied the regularity of
the solutions, and Banerjee and Budhiraja, (2020); Freidlin, (1985) present
the existence of a unique invariant measure of RSDEs. RSDEs have been analyzed
from the viewpoint of partial differential equations (PDE), for instance, by
Miranda, (1970); Cerrai, (1998); Barbu and Prato, (2005). The discretization
error of RSDEs has been a key component to connect the analysis on SDEs with
the convergence rate to optimization algorithms; Leimkuhler et al., (2020);
Cattiaux et al., (2017); Ding and Zhang, (2008) studied this aspect.

##  2 Our Proposed Algorithm: RGLD

Here, we describe the details of the reflected gradient Langevin dynamics
(RGLD) method. To properly describe the algorithm, we define some notations.

##### Notation

We call a random variable $\xi$ a Rademacher random variable if it satisfies
$\mathbb{P}(\xi=+1)=\mathbb{P}(\xi=-1)=1/2$, and call a vector-valued random
variable a Rademacher random vector if each element of the vector is an
independent Rademacher random variable. For a compact set
$K\subset\mathbb{R}^{d}$, we define the projection $\mathcal{P}_{K}$ and
reflection $\mathcal{R}_{K}$ respectively as

$\displaystyle\mathcal{P}_{K}(x):=\operatorname*{\arg\min}_{y\in
K}\norm{y-x},\quad\mathcal{R}_{K}(x):=2\mathcal{P}_{K}(x)-x.$ (2)

Let $\partial K,K^{\mathrm{o}}$, $K^{c}$ and $\bar{K}$ denote the boundary,
interior, complement and closure of $K$ respectively. Denote by
$K_{-r^{\prime}}$ the set of the points outside $K$ whose distance from $K$
are at most $r^{\prime}>0$. $\tilde{\mathcal{O}}$ represents the order symbol
ignoring the poly-log arithmetric factors, that is, for sequences $(f_{n})$
and $(g_{n})$, $f_{n}=\tilde{\mathcal{O}}(g_{n})$ means that there exists
constants $C,c>0$ such that $|f_{n}|\leq C|g_{n}|\log^{c}(|g_{n}|)$ for all
$n$. Similarly, we write $f\preceq g$ for positive functions $f$ and $g$ if
there exists a constant $C>0$ such that $f\leq Cg$, and write $f\sim g$ if
$f\preceq g\land f\succeq g.$ Finally, we let $[n]:=\\{0,1,\ldots,n-1\\}$.

The most common approach to a constrained optimization problem would be the
projected gradient descent (PG), the procedure of which is given by

$\displaystyle X_{k+1}=\mathcal{P}_{K}\left(X_{k}-\eta\gradient
f(X_{k})\right).$

However, PG may well be stuck on a saddle point or a suboptimal local optimal
solution. To escape from such a globally suboptimal critical point, it is
helpful to add a stochastic perturbation. The gradient Langevin dynamics (GLD)
is one of the most common approaches in that direction (Welling and Teh,,
2011; Raginsky et al.,, 2017). To adapt GLD to a constrained optimization
problem, GLD is combined with the projected operation as in the projected GLD
(PGLD) algorithm and its theoretical property has been analyzed especially
with a convex constraint (Bubeck et al.,, 2018; Lamperski,, 2021).
Conventionally, Gaussian noise is used to approximate the Brownian motion term
appearing in the continuous-time counterpart of GLD but it sometimes creates
so huge a noise that the solution moves far away from the feasible region.
Therefore, we employ the Rademacher random vector instead of the Gaussian
noise perturbation to keep the solution close to the domain since a Rademacher
random variable is bounded.

While the projection step has been commonly used to keep the feasibility
condition, we employ the reflection operation instead because we find that it
can cancel out the second-order term in the discretization error analysis and
achieves higher-order approximation error, which will be fully discussed in
Section C.1. This is analogous to the fact that the reflection step in the
Douglas-Rachford algorithm speeds up convergence in the feasibility problem as
shown in Artacho et al., (2020). Keeping these notions in our mind, we arrive
at the following RGLD method which is analyzed in this paper.

######  Algorithm 2.1 (Reflected gradient Langevin dynamics (RGLD)).

$\displaystyle\begin{cases}X^{\prime}_{k+1}=X_{k}-\eta\gradient
f(X_{k})+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\\\
X_{k+1}=\mathcal{R}_{K}\left(X^{\prime}_{k+1}\right),\end{cases}$

where $\eta>0$ and $\beta>0$ are the step size and the inverse temperature
parameter respectively, and $\xi_{k}$ is an i.i.d. Rademacher random vector in
$\mathbb{R}^{d}$.

Note that the boundedness of the Rademacher noise assures that the reflection
is uniquely defined for every iteration even if the domain is non-convex when
the boundary of the feasible region is sufficiently smooth. One has to reject
a huge increment by the Gaussian noise as in Bossy et al., (2004), which could
degrade the convergence.

RGLD algorithm (2.1) can be viewed as a time-discretized version of the
continuous-time reflected gradient Langevin diffusion described by

$\displaystyle\differential X_{t}=-\gradient f(X_{t})\differential
t+\sqrt{\frac{2}{\beta}}\differential W_{t}-v(X_{t})L(\differential
t),\quad(t\geq 0),$ (3)

where $(W_{t})$ is the standard $d$-dimensional Brownian motion, $v(X_{t})$ is
an outer unit normal vector at $X_{t}$, and $L$ is the local time supported on
the set $\\{t\in\mathbb{R}_{+}\mid X_{t}\in\partial K\\}$. Here, we say that
$v(x)$ is an outer unit normal vector at $x\in\partial K$ if $\norm{v(x)}=1$
and there exists $\epsilon>0$ such that

$\displaystyle\text{$x^{\prime}\in K$ satisfies
$\norm{x-x^{\prime}}<\epsilon$}\quad\Rightarrow\quad\langle
x-x^{\prime},v(x)\rangle\geq 0.$

We note that $v(x)=\bm{0}$ is the outer unit normal vector for $x\in
K^{\mathrm{o}}$, and thus RGLD is reduced to the standard gradient Langevin
dynamics on $K^{\mathrm{o}}$. We notice that there exists the reflection term
$-v(X_{t})L(\differential t)$ which is required to cancel out the velocity
directed to outside the region, as thoroughly discussed in Tanaka, (1979);
Lions and Sznitman, (1984); Dupuis and Ramanan, (1999). This term imposes a
non-trivial treatment when we discretize the dynamics. Indeed, there are at
least two choices, projection or reflection, to implement in the discrete-time
setting.

According to Bubeck et al., (2018), it is known that under suitable
conditions, the continuous-time RGLD (3) weakly converges to a unique
stationary distribution, the Gibbs distribution with the density function
$\pi$ defined by

$\displaystyle\pi(x)=\frac{\exp(-\beta f(x))}{\int_{K}\exp(-\beta
f(y))\differential y}\quad(x\in K),$ (4)

and $\pi(x)=0$ for $x\not\in K$. This explicit representation of the Gibbs
distribution enables us to bound the discrepancy between the expectation of
the objective with respect to the stationary distribution and the optimal
value as usual (Raginsky et al.,, 2017). Our main strategy is to show the weak
convergence of the discrete-time dynamics to the stationary distribution by
bounding the discretization error through the Poisson equation.

##  3 Convergence Analysis

In this section, we state the main theorem on the convergence rate of RGLD
under the assumptions given in Section 3.1 and give a sketch of its proof.

###  3.1 Assumptions

First, we prepare the assumptions on $K$ and $f$ for our convergence analysis.

######  Assumption 3.1 (Conditions on $K$).

  1. 1.

$\mathcal{P}_{K}$ is efficiently computable by some oracle as long as it is
unique.

  2. 2.

$K\subset\mathbb{R}^{d}$ is the closure of a possibly non-convex open and
connected domain such that $0\in K$ and, as a consequence, $K$ contains a
Euclidean ball of radius $r$.

  3. 3.

$K$ has a smooth boundary such that $\partial K\in C^{4}$.

  4. 4.

$K$ is bounded; that is, $K$ is included in a sphere of radius $R$ centered at
the origin.

  5. 5.

The problem (1) admits at least one optimal solution $x^{*}\in K$.

Assumption 3.1.1 is required to ensure that the reflection step in our
algorithm can be easily computed. Assumption 3.1.2 is necessary to guarantee
the uniqueness of the invariant law of the reflected diffusion process (3)
(Miranda,, 1970; Benchérif-Madani and Pardoux,, 2009), which is a key
component of the global optimization of RGLD. Moreover, it also allows a
reflection operation lest the operation moves the solution outside the domain
again. Assumption 3.1.3 is a relatively strict condition because many
constrained problems have multiple constraints and, thus, may have several
indifferentiable extreme points. An example of a domain with a smooth boundary
is a thick-walled sphere (Chepurnenko et al.,, 2021; Zhang et al., 2020a, ). A
Riemannian manifold such as the Stiefel manifold or the Grassmann manifold,
also satisfies the assumptions, covered by some of the existing works
(Patterson and Teh,, 2013; Wang et al.,, 2020). The $l_{q}$ norm ($q<1$) is
one of the most popular non-convex constraints for sparse estimation but it
violates Assumption 3.1.3 at points with zero elements. However, by applying
some smoothing operation, we may also make it have a smooth boundary.

Next, we describe the assumption on $f$.

######  Assumption 3.2 (Conditions on $f$).

  1. 1.

$f$ is a $C^{4}$ function.

  2. 2.

$\gradient f$ is accessible by some oracle.

From these assumptions, we can conclude that $f$ and the norm of $\gradient f$
are bounded by a constant $G$ because the feasible region $K$ is compact.
Moreover, $M$-smoothness of the objective function is easily established.

######  Proposition 3.3.

Under Assumptions 3.1.4 and 3.2.1, there exists a constant $M>0$ such that $f$
is $M$-smooth, that is, for any $x,y\in K$, it holds that $\norm{\gradient
f(x)-\gradient f(y)}\leq M\norm{x-y}$.

Under these assumptions, we can confirm that the reflection is unique for
sufficiently small $\eta$. Indeed, we can show
$\norm{\mathcal{P}_{K}(X^{\prime}_{k})-X^{\prime}_{k}}=O(\eta+\sqrt{d\eta/\beta})$,
which is sufficient to ensure the uniqueness under Assumption 3.1.3 (See Bossy
et al., (2004); Leimkuhler et al., (2020) for a more thorough statement).
Then, under Assumption 3.1.2, we can guarantee
$\mathcal{R}_{K}(X^{\prime}_{k})\in K$ for sufficiently small $\eta$, which
confirms the well-definedness of RGLD.

###  3.2 Main Result: Convergence Rate of RGLD

Here, we give our main result about the convergence rate of RGLD. In our
analysis, we utilize the geometric ergodicity of the continuous-time RGLD. The
rate of its weak convergence to the unique invariant measure is governed by
the spectral gap $\lambda_{*}$ defined by

$\displaystyle\lambda_{*}:=\inf_{\phi\in W^{1,2}(K^{o}),\phi\neq
0}\frac{\int_{K}\norm{\gradient\phi}^{2}\differential\pi}{\beta\int_{K}\phi^{2}\differential\pi},$
(5)

where $W^{1,2}(K^{o})$ is the Sobolev space such that the function itself and
its weak derivative have finite $L^{2}$ norms on $K^{o}$. Note that
$\lambda_{*}$ is the smallest, or the principal, eigenvalue of the
infinitesimal generator of $(X_{t})$ which is studied in, for example, Chen
and Lou, (2012); Peng et al., (2019).

Under this preparation, we obtain the following rate of convergence.

######  Theorem 3.4 (Main Theorem).

Suppose that Assumptions 3.1 and 3.2 hold. Then, for any $\beta>0$, the
following bound holds

$\displaystyle\mathbb{E}\left[\min_{k\in[N]}f(X_{k})\right]-\min f\leq
C\left(\frac{1}{\lambda_{*}N\eta}+\frac{\sqrt{\eta}\left(\beta\sqrt{\eta}+(\beta\eta+d)^{3/2}\right)}{\lambda_{*}}+\frac{d}{\beta}\log\beta\right),$

where $C>0$ is a constant independent of $\lambda_{*},N,\eta,\beta,d$. In
particular, for any $0<\epsilon\ll 1$, by setting each parameter as

$\displaystyle\beta>0,\;\eta\preceq\min\left\\{\frac{\lambda_{*}^{2}\epsilon^{2}}{d^{3}},\frac{\lambda_{*}\epsilon}{\beta}\right\\},\;N\succeq\frac{1}{\lambda_{*}\epsilon\eta}\succeq\max\left\\{\frac{d^{3}}{\lambda_{*}^{3}\epsilon^{3}},\frac{\beta}{\lambda_{*}^{2}\epsilon^{2}}\right\\},$

we have

$\displaystyle\mathbb{E}\left[\min_{k\in[N]}f(X_{k})\right]-\min
f\leq\epsilon+\frac{d\log\beta}{\beta}.$

Informally speaking, the main theorem states that, by selecting a sufficiently
small $\eta$, an $\epsilon$-sampling error is achieved after
$\tilde{O}\left(\max\left\\{{d^{3}}{\lambda_{*}^{-3}\epsilon^{-3}},{\beta}{\lambda_{*}^{-2}\epsilon^{-2}}\right\\}\right)$
iterations of RGLD. Now, remember that the iteration complexity for the PGLD
obtained by the previous work (Lamperski,, 2021) was evaluated as
$N=\tilde{O}(d^{4}\lambda_{*}^{-1}\epsilon^{-4})$ (see Table 1). This is worse
than our evaluation $N=\tilde{O}(d^{3}\lambda_{*}^{-3}\epsilon^{-3})$ when
$\beta$ is sufficiently large as $\beta=\tilde{\Omega}(d)$ (which is the
situation of the interest).

This advantage is mainly due to the employment of the reflection operator
instead of the projection. The discretization error can be evaluated by
approximating the time differentiation for the infinitesimal generator by a
discrete-time differentiation and in that approximation, we can cancel out the
second-order derivatives thanks to the symmetric relation between $X_{k}$ and
$X^{\prime}_{k}$. Hence, we obtain a sharper bound than that of PGLD. On the
other hand, Lamperski, (2021, Lemmas 5, 6) showed that the discretization
error between $f(X_{k})$ and $f(X_{k\eta})$ is $\tilde{O}(\eta^{1/4})$ which
is worse than our analysis. Although we might be able to improve the rate even
for PGLD, it would require ensuring the hitting time on the boundary is
negligibly small so that the difference of projection and reflection could be
ignored. Another reason of the suboptimality is the crude estimate of
1-Wasserstein metric by the pathwise distance
$\mathbb{E}\left[\|X_{k}-X_{k\eta}\|\right]$, which should be less tight than
the direct evaluation of the difference in a weak convergence sense. We cope
with this issue by exploiting the solution of the corresponding Poisson
equation.

######  Remark 3.5.

The dependence of $\lambda_{*}$ on $\beta$ is an open problem and has been
studied in Chen and Lou, (2012); Peng et al., (2019). $\lambda_{*}$ is
supposed to be proportional to the exponential of $\beta$ since the
exponential dependence is proved for unconstrained GLD under the similar
assumptions in Raginsky et al., (2017); Gayrard et al., (2005); Bardet et al.,
(2018). If that is true, the task of finding an $\epsilon$-optimal solution
suffers from the curse of dimensionality, where $\epsilon$-optimal solution
$x$ is defined as the solution which satisfies
$\mathbb{E}\left[f(x)\right]-\min_{K}f\leq\epsilon$. Indeed, it is necessary
to set $\beta\succeq d\epsilon^{-1}\log\epsilon^{-1}$ to assert
$\mathbb{E}\left[\min_{k\in[N]}f(X_{k})\right]-\min_{K}f\leq\epsilon$. A
simple calculation shows that this implies
$N\succeq\exp(c_{1}d\epsilon^{-1})$. $N$ depends on both the exponential of
the dimensionality $d$ and that of the inverse of the admitted error
$\epsilon$.

###  3.3 Proof Sketch of Theorem 3.4

In this subsection, we give a proof sketch of Theorem 3.4. We will denote any
constant independent of $\lambda_{*},N,\eta,\beta,d$ by the same notation $C$
for simplicity. Theorem 3.4 immediately follows from the following lemma by
noticing $\min_{k\in[N]}f(X_{k})\leq\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})$.

######  Lemma 3.6.

Under Assumptions 3.1 and 3.2, if $\beta>0$, then it holds that

$\displaystyle\mathbb{E}\left[\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})\right]-\min
f\leq
C\left(\frac{1}{\lambda_{*}N\eta}+\frac{\sqrt{\eta}\left(\beta\sqrt{\eta}+(\beta\eta+d)^{3/2}\right)}{\lambda_{*}}+\frac{d}{\beta}\log(C\beta)\right).$

Given this evaluation, we can easily obtain the parameter configuration in
Theorem 3.4 to obtain a solution with $\epsilon$-sampling error by controlling
the first two terms on the right-hand side to make them smaller than
$\epsilon$. For completeness, we attach the proof of Theorem 3.4 from Lemma
3.6 in Appendix A.

Hereafter, we give a proof of Lemma 3.6. Recall that the algorithm
$\\{X_{k}\\}$ is a time-discretized version of the continuous-time version
$(X_{t})$ following the SDE (3). One finds that $\left(X_{t}\right)$
geometrically converges in distribution to the Gibbs distribution $\pi$ given
in (4) (see Lemma C.4). It can be shown that the distribution $\pi$
concentrates around $\operatorname*{\arg\min}_{x\in K}f(x)$. Thereby, we just
need to prove that $\\{X_{k}\\}$ well approximates $\pi$ for sufficiently
large $k$ and small $\eta$. Along with this idea, we decompose the expected
risk as follows: $\mathbb{E}[N^{-1}\sum_{k=0}^{N-1}f(X_{k})]-\min
f=(\mathbb{E}[N^{-1}\sum_{k=0}^{N-1}f(X_{k})]-\mathbb{E}_{\pi}f)+(\mathbb{E}_{\pi}f-\min
f)$. Then, each term on the right hand side can be bounded in the following
two lemmas.

######  Lemma 3.7.

Under Assumptions 3.1 and 3.2, if $\beta>0$, then it holds that

$\displaystyle\mathbb{E}\left[\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})\right]-\mathbb{E}_{\pi}f\leq
C\left(\frac{1}{\lambda_{*}N\eta}+\frac{\sqrt{\eta}\left(\beta\sqrt{\eta}+(\beta\eta+d)^{3/2}\right)}{\lambda_{*}}\right).$
(6)

######  Lemma 3.8 (Lemma 15 from Lamperski, (2021)).

Under Assumptions 3.1.4 and 3.2.1, we have

$\displaystyle\mathbb{E}_{\pi}f-\min_{x\in
K}f(x)\leq\frac{d}{\beta}\log(2R\max\left\\{\frac{2}{r},\frac{L\beta\left(r+\sqrt{r^{2}+R^{2}}\right)}{r\log{2}}\right\\}),$

where $L$ is the Lipschitz constant of $f$.

We can easily see that Lemma 3.6 follows after the above two lemmas. Lemma 3.8
evaluates how well the expectation of the objective with respect to $\pi$
approximates the minimum value of the given function. The reader may refer to
Lamperski, (2021) for its proof. A similar formula of this bound was also
shown by Raginsky et al., (2017) where an unconstrained problem was analyzed.

Lemma 3.7 is the key lemma for our purpose to show that the sequence generated
by RGLD approximates sampling from the invariant law $\pi$ given a
sufficiently small step size $\eta$ and a sufficiently large number of
iterations $N$.

###  3.4 Strategy for Proving Lemma 3.7

Remember that the continuous-time RGLD follows the SDE given in (3). It is
very involved to prove that the process $(X_{t})$ stays on $K$ almost surely
due to the reflection term. The existence and uniqueness of reflected SDEs
have been studied under the name of the Skorokhod problem; see Tanaka, (1979);
Lions and Sznitman, (1984).

In a similar way to the unconstrained setting, we can ensure the existence and
uniqueness of the invariant measure of the SDE (3) (see e.g., Bubeck et al.,
(2018)). A concise proof is given in Appendix B for completeness.

######  Lemma 3.9 (Existence of Unique Invariant Measure).

Under Assumptions 3.1.3 and 3.2.1, $(X_{t})$ has a unique invariant measure
$\pi$, where $\pi$ is the Gibbs measure with a density defined by (4).

Then, our purpose is reduced to showing the weak convergence of the time
discretized version and bound the error for finite $k$ and positive $\eta$. To
do so, we consider the following Neumann problem associated with the Poisson
equation:

$\displaystyle\begin{cases}\mathcal{A}u(x)=f(x)-\mathbb{E}_{\pi}f,&x\in K,\\\
\gradient u(x)\cdot v(x)=0,&x\in\partial K,\end{cases}$ (7)

where $\mathcal{A}$ is the infinitesimal generator of the dynamics defined as
$\mathcal{A}u(x)=\beta^{-1}\Delta u(x)-\langle\gradient f(x),\gradient
u(x)\rangle$ for $x\in K$ where $\Delta$ is the Laplace operator. As is known
in Miranda, (1970); Freidlin, (1985), the solvability of this problem requires
the compatibility condition,

$\displaystyle\int_{K}\left(f(x)-\mathbb{E}_{\pi}f\right)\differential\pi(x)=0,$

which holds by definition. From Assumptions 3.1.3 and 3.2.1, the Neumann
problem has a unique solution (up to an additive constant) $u\in C^{4}(K)$
(see Miranda, (1970), Gilbarg and Trudinger, (1983, Chapter 6), and Lieberman,
(1986, Theorem 3)). Indeed, following Benchérif-Madani and Pardoux, (2009,
Theorem 4), we have that the solution can be written as
$u(x)=\int_{0}^{\infty}\mathbb{E}\left[f(X_{t}^{x})-\mathbb{E}_{\pi}f\right]\differential
t$ where $X_{t}^{x}$ denotes the solution of (3) at time $t$ with the initial
value $x$. Then, by the Ito-formula and (7), we have that

$\frac{\partial\mathbb{E}[u(X_{t}^{x})]}{\partial
t}\Big{|}_{t=0}=\mathcal{A}u(x)=f(x)-\mathbb{E}_{\pi}f.$

By substituting $x=X_{k}$, we obtain an estimate of
$f(X_{k})-\mathbb{E}_{\pi}f$ through the solution of the Poisson equation.
Thus, by taking the average over $k=0,\dots,N-1$, we have
$\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})-\mathbb{E}_{\pi}f=\frac{1}{N}\sum_{k=0}^{N-1}{\partial_{t}\mathbb{E}[u(X_{t}^{X_{k}})]}\big{|}_{t=0}$,
which is what we want to evaluate. Then, intuitively speaking, we approximate
each term ${\partial_{t}\mathbb{E}[u(X_{t}^{X_{k}})]}\big{|}_{t=0}$ on the
right hand side as
${\partial_{t}\mathbb{E}[u(X_{t}^{X_{k}})]}\big{|}_{t=0}\approx{\eta}^{-1}\left(\mathbb{E}[u(X_{k+1})]-\mathbb{E}[u(X_{k})]\right)$,
and arrive at
$\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})-\mathbb{E}_{\pi}f\approx\frac{1}{N\eta}\left(\mathbb{E}[u(X_{N})]-\mathbb{E}[u(X_{0})]\right)$.
The right hand side of this equation can be bounded by
$O(1/(\lambda_{*}N\eta))$ from Lemma C.1, which corresponds to the first term
of the bound (6) in Lemma 3.7. The second term of the bound (6) stems from the
discretization error which can be evaluated as follows.

To evaluate the discretization error, we employ a strategy that is parallel to
that of Leimkuhler et al., (2020, Theorem 4.2), but we additionally require
the parameter dependence of their bounds on $\beta$ and $d$. For that purpose,
we adapt the Taylor expansion of $u(X_{k+1})-u(X_{k})$ up to the third order
with respect to $X_{k+1}-X_{k}$. Then, we bound the derivatives, which is non-
trivial, by employing the Bismut-Elworthy-Li formula borrowed from Andres,
(2011); Cerrai, (1998) for our Neumann problem. In that evaluation, we need to
show that the derivatives of
$\mathbb{E}\left[f(X_{t}^{x})\right]-\mathbb{E}_{\pi}f$ with respect to $x$
decays exponentially as $t$ goes to infinity. We prove it by extending the
argument made in the proof of Bréhier, (2014, Lemma 5.4). The detailed proof
of Lemma 3.7 is deferred to Appendix C. Fortunately, the reflection operation
enables us to cancel the second-order derivative term and we can obtain a
tighter evaluation of the discretization error.

##  4 Numerical Experiments

We conducted the numerical experiments to see how well RGLD performs on
constrained non-convex problems in comparison with a classical algorithm, PG.

We quantified the optimization error by the gap of objective function values
at the global optimal solution and the solution of the algorithms. Our focus
was on how efficiently RGLD, compared with PG, finds the global optimum of a
constrained non-convex problem. We also compared the performances of RGLD
under different hyperparameters $\eta$, $\beta$, and dimensionalities $d$.

The latter half of the experiments is described in Appendix D.

###  4.1 Two Dimensional Problem

We start by looking at a two dimensional problem for ease of visualization.
The test function is a Gaussian mixture density given by

$\displaystyle
f(x)=-\sum_{i=1}^{M}w_{i}\exp(-\frac{1}{2}(x-m_{i})^{\top}(x-m_{i})),$ (8)

where $w_{i}$ is randomly generated and $m_{i}$ is a point in the meshgrid
$\\{-2,-1,0,1,2\\}^{2}$.

Figure 1: The objective function $f$ in (8). Figure 2: Convergence of the
iterates by PG and RGLD. $\eta=0.05,\>\beta=1.0$.
---

The feasible region is set to the region between spheres of radius 0.9 and 4
centered at the origin. Such a region is called a thick-walled sphere and
appears in Chepurnenko et al., (2021); Zhang et al., 2020a . The region is
non-convex and satisfies our assumptions because the boundary is smooth.
Figure 2 shows a 3D plot of the objective function, which has 25 local minima.
The global minimizer is at $(0,-2)$.

####  4.1.1 Comparison with PG

We found that RGLD can explore better solutions while PG becomes trapped at a
local minimum. Figure 2 shows the plot of the optimization error of PG and
RGLD for $\eta=0.05$, $\beta=1.0$ and $X_{0}=(0.5,0.5)$.

We also visualized the trajectory of each algorithm, as shown in Figure 3. We
found that RGLD explores a wide range of the feasible region but seems not to
converge to a global minimum at all. However, as shown in Figure 4, a larger
$\beta$ stabilizes the optimization, which causes RGLD to become stuck in each
mode of the objective function.

PG

RGLD

Figure 3: Trajectories of the iterates by PG and RGLD.

$\beta=1.0$.

$\beta=8.0$.

Figure 4: Trajectory of RGLD with different inverse temperature parameters
$\beta$.

####  4.1.2 Sensitivity to Hyperparameters

Our main theorem provides an insight on the appropriate $\eta$ and $\beta$.
This motivated us to examine how the values of each hyperparameter could
affect the optimization. Here, we plot the cumulative minimum of the
optimization error generated by the iterates of RGLD for ease of viewing.
Figures 7 and 7 show the convergence of RGLD for different values of $\beta$
or $\eta$. We can see that a larger $\beta$ limits the solution to a narrower
range of the region, resulting in slow convergence. As for $\eta$, there is a
certain trade-off; that is, we need to set $\eta$ large enough to search a
broad range of the region within a limited number of iterates and, at the same
time, the range has to be small enough to find narrow loss valleys.

####  4.1.3 Comparison of Reflection and Projection

We also compare the effects of the reflection operation and the projection
operation. Figure 7 shows the convergence of the iterates generated by RGLD
and PGLD. The reflection is preferable in our theoretical analysis but the
result shows that these two operations make little difference in practice.

Figure 5: RGLD under
different $\beta$. $\eta=.05$. Figure 6: RGLD under
different $\eta$. $\beta=1.0$. Figure 7: Comparison between RGLD and PGLD.
---

##  5 Conclusion and Future Work

We proved the sub-linear convergence rate of RGLD in solving smoothly
constrained problems. Additionally, the rate faster than that of Lamperski,
(2021) is obtained by employing the reflection and the direct estimate of the
discretization error by the Poisson equation. However, the order of
$\lambda_{*}$ is still unclear because the smallest eigenvalue problems with
Neumann boundary conditions have not been studied well. Moreover, the
smoothness of $\partial K$ is essential in the analysis but not satisfied in
many applications, though a smoothing operation can make a non-smooth domain
have a smooth boundary.

The future works will alleviate the currently imposed assumptions towards more
general unbounded non-smooth constraints. One common approach is to assume a
dissipativity condition instead of the boundedness of $K$, as imposed in
Raginsky et al., (2017); Xu et al., (2018). In this case, Barbu and Prato,
(2005) could be a key component with which to establish the necessary
theoretical background. It is also important to extend our result to the
stochastic gradient oracle version.

## Acknowledgments

I would like to express my very great appreciation to our former research
fellow, Michael Metel for his valuable and constructive suggestions and
cooperations during the planning and development of this research. TS was
partially supported by JSPS KAKENHI (18H03201), Japan Digital Design and JST
CREST.

## References

  * Andres, (2011) Andres, S. (2011).  Pathwise differentiability for SDEs in a smooth domain with reflection.  Electronic Journal of Probability, 16:845–879.
  * Artacho et al., (2020) Artacho, F. J. A., Campoy, R., and Tam, M. K. (2020).  The Douglas–Rachford algorithm for convex and nonconvex feasibility problems.  Mathematical Methods of Operations Research, 91(2):201–240.
  * Banerjee and Budhiraja, (2020) Banerjee, S. and Budhiraja, A. (2020).  Parameter and dimension dependence of convergence rates to stationarity for reflecting Brownian motions.  The Annals of Applied Probability, 30(5):2005–2029.
  * Barbu and Prato, (2005) Barbu, V. and Prato, G. D. (2005).  The Neumann problem on unbounded domains of $\mathbb{R}^{d}$ and stochastic variational inequalities.  Communications in Partial Differential Equations, 30(8):1217–1248.
  * Bardet et al., (2018) Bardet, J.-B., Gozlan, N., Malrieu, F., and Zitt, P.-A. (2018).  Functional inequalities for Gaussian convolutions of compactly supported measures: Explicit bounds and dimension dependence.  Bernoulli, 24(1):333–353.
  * Barkhagen et al., (2019) Barkhagen, M., Fleming, B., Quiles, S. G., Gondzio, J., Kalcsics, J., Kroeske, J., Sabanis, S., and Staal, A. (2019).  Optimising portfolio diversification and dimensionality.  arXiv preprint arXiv:1906.00920.
  * Benchérif-Madani and Pardoux, (2009) Benchérif-Madani, A. and Pardoux, É. (2009).  A probabilistic formula for a Poisson equation with Neumann boundary condition.  Stochastic Analysis and Applications, 27(4):739–746.
  * Bossy et al., (2004) Bossy, M., Gobet, E., and Talay, D. (2004).  A symmetrized Euler scheme for an efficient approximation of reflected diffusions.  Journal of Applied Probability, 41(3):877–889.
  * Bréhier, (2014) Bréhier, C.-E. (2014).  Approximation of the invariant measure with an Euler scheme for stochastic PDEs driven by space-time white noise.  Potential Analysis, 40(1):1–40.
  * Brosse et al., (2017) Brosse, N., Durmus, A., Moulines, É., and Pereyra, M. (2017).  Sampling from a log-concave distribution with compact support with proximal Langevin Monte Carlo.  In Kale, S. and Shamir, O., editors, Proceedings of the 2017 Conference on Learning Theory, volume 65 of Proceedings of Machine Learning Research, pages 319–342. PMLR.
  * Bubeck et al., (2018) Bubeck, S., Eldan, R., and Lehec, J. (2018).  Sampling from a log-concave distribution with projected Langevin Monte Carlo.  Discrete & Computational Geometry, 59(4):757–783.
  * Cattiaux et al., (2017) Cattiaux, P., León, J. R., and Prieur, C. (2017).  Invariant density estimation for a reflected diffusion using an Euler scheme.  Monte Carlo Methods and Applications, 23(2):71–88.
  * Cerrai, (1998) Cerrai, S. (1998).  Some results for second order elliptic operators having unbounded coefficients.  Differential and Integral Equations, 11(4):561–588.
  * Chen et al., (2021) Chen, P., Lu, J., and Xu, L. (2021).  Approximation to stochastic variance reduced gradient Langevin dynamics by stochastic delay differential equations.  arXiv preprint arXiv:2106.04357.
  * Chen and Lou, (2012) Chen, X. and Lou, Y. (2012).  Effects of diffusion and advection on the smallest eigenvalue of an elliptic operator and their applications.  Indiana University Mathematics Journal, 61(1):45–80.
  * Chepurnenko et al., (2021) Chepurnenko, A., Litvinov, S., Meskhi, B., and Beskopylny, A. (2021).  Optimization of thick-walled viscoelastic hollow polymer cylinders by artificial heterogeneity creation: Theoretical aspects.  Polymers, 13(15).
  * Chiang et al., (1987) Chiang, T.-S., Hwang, C.-R., and Sheu, S. J. (1987).  Diffusion for global optimization in $\mathbb{R}^{n}$.  SIAM Journal on Control and Optimization, 25(3):737–753.
  * Dalalyan, (2017) Dalalyan, A. (2017).  Further and stronger analogy between sampling and optimization: Langevin Monte Carlo and gradient descent.  In Kale, S. and Shamir, O., editors, Proceedings of the 2017 Conference on Learning Theory, volume 65 of Proceedings of Machine Learning Research, pages 678–689. PMLR.
  * Ding and Zhang, (2008) Ding, D. and Zhang, Y.-y. (2008).  Numerical solutions for reflected stochastic differential equations.
  * Dupuis and Ramanan, (1999) Dupuis, P. and Ramanan, K. (1999).  Convex duality and the Skorokhod problem. I.  Probability Theory and Related Fields, 115(2):153–195.
  * Durmus et al., (2018) Durmus, A., Moulines, E., and Pereyra, M. (2018).  Efficient Bayesian computation by proximal Markov chain Monte Carlo: when Langevin meets Moreau.  SIAM Journal on Imaging Sciences, 11(1):473–506.
  * Freidlin, (1985) Freidlin, M. I. (1985).  Functional Integration and Partial Differential Equations, volume 109.  Princeton University Press.
  * Garroni and Menaldi, (1992) Garroni, M. G. and Menaldi, J. L. (1992).  Green Functions for Second Order Parabolic Integro-Differential Problems, volume 275.  Chapman and Hall/CRC.
  * Gayrard et al., (2005) Gayrard, V., Bovier, A., and Klein, M. (2005).  Metastability in reversible diffusion processes II: Precise asymptotics for small eigenvalues.  Journal of the European Mathematical Society, 7(1):69–99.
  * Gelfand and Mitter, (1991) Gelfand, S. B. and Mitter, S. K. (1991).  Recursive stochastic algorithms for global optimization in $\mathbb{R}^{d}$.  SIAM Journal on Control and Optimization, 29(5):999–1018.
  * Geman and Hwang, (1986) Geman, S. and Hwang, C.-R. (1986).  Diffusions for global optimization.  SIAM Journal on Control and Optimization, 24(5):1031–1043.
  * Gilbarg and Trudinger, (1983) Gilbarg, D. and Trudinger, N. S. (1983).  Elliptic Partial Differential Equations of Second Order.  Springer.
  * Hinder and Ye, (2018) Hinder, O. and Ye, Y. (2018).  Worst-case iteration bounds for log barrier methods for problems with nonconvex constraints.  arXiv preprint arXiv:1807.00404.
  * Hsieh et al., (2018) Hsieh, Y.-P., Kavis, A., Rolland, P., and Cevher, V. (2018).  Mirrored Langevin dynamics.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 31, pages 2878–2887. Curran Associates, Inc.
  * Kingma and Ba, (2014) Kingma, D. P. and Ba, J. (2014).  Adam: a method for stochastic optimization.  arXiv preprint arXiv:1412.6980.
  * Lamperski, (2021) Lamperski, A. (2021).  Projected stochastic gradient Langevin algorithms for constrained sampling and non-convex learning.  In Belkin, M. and Kpotufe, S., editors, Proceedings of Thirty Fourth Conference on Learning Theory, volume 134 of Proceedings of Machine Learning Research, pages 2891–2937. PMLR.
  * Leimkuhler et al., (2020) Leimkuhler, B., Sharma, A., and Tretyakov, M. V. (2020).  Simplest random walk for approximating Robin boundary value problems and ergodic limits of reflected diffusions.  arXiv preprint arXiv:2006.15670.
  * Li and Lin, (2015) Li, H. and Lin, Z. (2015).  Accelerated proximal gradient methods for nonconvex programming.  In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 28, pages 379–387. Curran Associates, Inc.
  * Lieberman, (1986) Lieberman, G. M. (1986).  Intermediate Schauder estimates for oblique derivative problems.  Archive for Rational Mechanics and Analysis, 93(2):129–134.
  * Lin et al., (2019) Lin, Q., Ma, R., and Xu, Y. (2019).  Inexact proximal-point penalty methods for constrained non-convex optimization.  arXiv preprint arXiv:1908.11518.
  * Lions and Sznitman, (1984) Lions, P.-L. and Sznitman, A.-S. (1984).  Stochastic differential equations with reflecting boundary conditions.  Communications on Pure and Applied Mathematics, 37(4):511–537.
  * Masnadi-shirazi and Vasconcelos, (2009) Masnadi-shirazi, H. and Vasconcelos, N. (2009).  On the design of loss functions for classification: theory, robustness to outliers, and savageboost.  In Koller, D., Schuurmans, D., Bengio, Y., and Bottou, L., editors, Advances in Neural Information Processing Systems, volume 21, pages 1049–1056. Curran Associates, Inc.
  * Miranda, (1970) Miranda, C. (1970).  Partial Differential Equations of Elliptic Type, volume 2.  Springer.
  * Parpas et al., (2006) Parpas, P., Rustem, B., and Pistikopoulos, E. N. (2006).  Linearly constrained global optimization and stochastic differential equations.  Journal of Global Optimization, 36(2):191–217.
  * Patterson and Teh, (2013) Patterson, S. and Teh, Y. W. (2013).  Stochastic gradient Riemannian Langevin dynamics on the probability simplex.  In Burges, C. J. C., Bottou, L., Welling, M., Ghahramani, Z., and Weinberger, K. Q., editors, Advances in Neural Information Processing Systems, volume 26, pages 3102–3110. Curran Associates, Inc.
  * Peng et al., (2019) Peng, R., Zhang, G., and Zhou, M. (2019).  Asymptotic behavior of the principal eigenvalue of a linear second order elliptic operator with small/large diffusion coefficient.  SIAM Journal on Mathematical Analysis, 51(6):4724–4753.
  * Pereyra, (2016) Pereyra, M. (2016).  Proximal Markov chain Monte Carlo algorithms.  Statistics and Computing, 26(4):745–760.
  * Raginsky et al., (2017) Raginsky, M., Rakhlin, A., and Telgarsky, M. (2017).  Non-convex learning via stochastic gradient Langevin dynamics: a nonasymptotic analysis.  In Proceedings of the 2017 Conference on Learning Theory, volume 65 of Proceedings of Machine Learning Research, pages 1674–1703.
  * Salim et al., (2019) Salim, A., Kovalev, D., and Richtárik, P. (2019).  Stochastic proximal Langevin algorithm: Potential splitting and nonasymptotic rates.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32, pages 6653–6664. Curran Associates, Inc.
  * Salim and Richtarik, (2020) Salim, A. and Richtarik, P. (2020).  Primal dual interpretation of the proximal stochastic gradient Langevin algorithm.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 3786–3796. Curran Associates, Inc.
  * Signoretto et al., (2013) Signoretto, M., De Lathauwer, L., and Suykens, J. A. K. (2013).  Learning tensors in reproducing kernel Hilbert spaces with multilinear spectral penalties.  arXiv preprint arXiv:1310:4977.
  * Suzuki et al., (2016) Suzuki, T., Kanagawa, H., Kobayashi, H., Shimizu, N., and Tagami, Y. (2016).  Minimax optimal alternating minimization for kernel nonparametric tensor learning.  In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 29, pages 3783–3791. Curran Associates, Inc.
  * Tanaka, (1979) Tanaka, H. (1979).  Stochastic differential equations with reflecting boundary condition in convex regions.  Hiroshima Mathematical Journal, 9(1):163–177.
  * Vempala and Wibisono, (2019) Vempala, S. and Wibisono, A. (2019).  Rapid convergence of the unadjusted Langevin algorithm: Isoperimetry suffices.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32, pages 8094–8106. Curran Associates, Inc.
  * Wang et al., (2020) Wang, X., Lei, Q., and Panageas, I. (2020).  Fast convergence of Langevin dynamics on manifold: Geodesics meet log-Sobolev.  In Larochelle, H., Ranzato, M., Hadsell, R., Balcan, M. F., and Lin, H., editors, Advances in Neural Information Processing Systems, volume 33, pages 18894–18904. Curran Associates, Inc.
  * Welling and Teh, (2011) Welling, M. and Teh, Y. W. (2011).  Bayesian learning via stochastic gradient Langevin dynamics.  In Proceedings of the 28th International Conference on Machine Learning, pages 681–688.
  * Wen et al., (2018) Wen, F., Chu, L., Liu, P., and Qiu, R. C. (2018).  A survey on nonconvex regularization-based sparse and low-rank recovery in signal processing, statistics, and machine learning.  IEEE Access, 6:69883–69906.
  * Xie and Wright, (2021) Xie, Y. and Wright, S. J. (2021).  Complexity of proximal augmented Lagrangian for nonconvex optimization with nonlinear equality constraints.  Journal of Scientific Computing, 86(3):1–30.
  * Xu et al., (2018) Xu, P., Chen, J., Zou, D., and Gu, Q. (2018).  Global convergence of Langevin dynamics based algorithms for nonconvex optimization.  In Bengio, S., Wallach, H., Larochelle, H., Grauman, K., Cesa-Bianchi, N., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 31, pages 3122–3133. Curran Associates, Inc.
  * (55) Zhang, J., Oueslati, A., Shen, W., De Saxcé, G., Nguyen, A. D., Zhu, Q., and Shao, J. (2020a).  Shakedown analysis of a hollow sphere by interior-point method with non-linear optimization.  International Journal of Mechanical Sciences, 175:105515.
  * (56) Zhang, K. S., Peyré, G., Fadili, J., and Pereyra, M. (2020b).  Wasserstein control of mirror Langevin Monte Carlo.  In Abernethy, J. and Agarwal, S., editors, Proceedings of Thirty Third Conference on Learning Theory, volume 125 of Proceedings of Machine Learning Research, pages 3814–3841. PMLR.
  * Zhang et al., (2018) Zhang, X., Wang, L., and Gu, Q. (2018).  A unified framework for nonconvex low-rank plus sparse matrix recovery.  In Storkey, A. and Perez-Cruz, F., editors, Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pages 1097–1107. PMLR.
  * Zhang et al., (2017) Zhang, Y., Liang, P., and Charikar, M. (2017).  A hitting time analysis of stochastic gradient Langevin dynamics.  In Proceedings of the 2017 Conference on Learning Theory, volume 65 of Proceedings of Machine Learning Research, pages 1980–2022.
  * Zou et al., (2021) Zou, D., Xu, P., and Gu, Q. (2021).  Faster convergence of stochastic gradient Langevin dynamics for non-log-concave sampling.  In de Campos, C. and Maathuis, M. H., editors, Proceedings of the Thirty-Seventh Conference on Uncertainty in Artificial Intelligence, volume 161 of Proceedings of Machine Learning Research, pages 1152–1162. PMLR.

##  Appendix A Proof of Theorem 3.4

First, it is obvious that

$\displaystyle\min_{k\in[N]}f(X_{k})\leq\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k}),\quad\mathbb{E}\left[\min_{k\in[N]}f(X_{k})\right]-\min
f\leq\mathbb{E}\left[\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})\right]-\min f.$

Thus, it suffices to show

$\displaystyle\mathbb{E}\left[\frac{1}{N}\sum_{k=0}^{N-1}f(X_{k})\right]-\min
f\preceq\epsilon+\frac{d}{\beta}\log\beta$

to conclude that $\mathbb{E}\left[\min_{k\in[N]}f(X_{k})\right]-\min
f\leq\epsilon+d\beta^{-1}\log\beta$. From Lemma 3.6, it is enough to establish

$\displaystyle\frac{1}{\lambda_{*}N\eta}+\frac{\sqrt{\eta}\left(\beta\sqrt{\eta}+(\beta\eta+d)^{3/2}\right)}{\lambda_{*}}\preceq\epsilon.$

From the condition on $\eta$, we have
${\lambda_{*}}^{-1}{\beta\eta}\preceq\epsilon$,
${\lambda_{*}}^{-1}(d\eta)^{3/2}\preceq\epsilon$. Moreover, we have
$\lambda_{*}^{-1}\beta^{3/2}\eta^{2}\preceq\lambda_{*}\epsilon^{2}/\sqrt{\beta}$
and, noticing that $\lambda_{*},\epsilon,\beta^{-1/2}\preceq 1$, we have
$\lambda_{*}\epsilon^{2}/\sqrt{\beta}\preceq\epsilon$. Finally, from the
condition on $N$, we have $(\lambda_{*}N\eta)^{-1}\preceq\epsilon$.

The theorem follows by combining all the above.

##  Appendix B Proof of Lemma 3.9

From Assumptions 3.1.3 and 3.2.1, there exists a unique invariant measure of
$\left(X_{t}\right)$. Thus, it suffices to show that $\pi$ is stationary. It
holds if and only if

$\displaystyle\int_{K}\mathcal{A}u\differential\pi=0\quad\forall
u\in\mathrm{dom}\mathcal{A},$

where $\mathcal{A}$ is the infinitesimal generator of $(X_{t})$. The
infinitesimal generator $\mathcal{A}$ is represented by

$\displaystyle\mathcal{A}u(x)=\frac{1}{\beta}\Delta u(x)-\gradient
f(x)\cdot\gradient u(x)\quad x\in K,$

where $\Delta$ is the Laplace operator, and the domain of $\mathcal{A}$ is

$\displaystyle\mathrm{dom}\mathcal{A}=\left\\{u\mid\gradient u(x)\cdot
v(x)=0,\quad x\in\partial K\right\\}.$

Then, one can derive

$\displaystyle\int_{K}\mathcal{A}u\differential\pi$
$\displaystyle=\int_{K}\left(\frac{1}{\beta}\Delta u(x)-\gradient
f(x)\cdot\gradient u(x)\right)\frac{\exp(-\beta f(x))}{\int_{K}\exp(-\beta
f(y))\differential y}\differential x$
$\displaystyle=\int_{K}\divergence(\frac{\exp(-\beta f(x))\gradient
u(x)}{\beta\int_{K}\exp(-\beta f(y))\differential y})\differential x$
$\displaystyle=\frac{1}{\beta\int_{K}\exp(-\beta f(y))\differential
y}\int_{\partial K}\exp(-\beta f(x))\gradient u(x)\cdot v(x)\differential
x=0,$

from the definition of $\pi$ and $\mathcal{A}$, the divergence theorem, and
the condition of $u$.

##  Appendix C Proof of Lemma 3.7

Lemma 3.7Lemma C.1Lemma C.2Lemma C.3Lemma C.4 Figure 8: Relations between
Statements Used to Prove Lemma 3.7

###  C.1 Error Analysis via Taylor Expansion

In this part, we aim to establish a relationship between the Taylor expansion
of the solution $u$ of (7) and the error analysis of the invariant measure
approximation. First, from Garroni and Menaldi, (1992, Proposition 1.17), it
is known that $u\in C^{4}(K)$ can be extended to a function in $C^{4}(K\cup
K_{-r^{\prime}})$, where $r^{\prime}$ is the supremum of positive numbers such
that the reflection of any $x\in K_{-r^{\prime}}$ is well-defined. This
implies that $u$ and its derivatives up to fourth order are uniformly bounded
on $K\cup K_{-r^{\prime}}$. We often identify the $i$-th derivative of $u$
with a multilinear operator on $\mathbb{R}^{d}$, $D^{i}u$, by using the Riesz
representation theorem. For example, we identify $\Delta u(x)$ at $x\in K\cup
K_{-r^{\prime}}$ with the bilinear operator
$D^{2}u(x):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$. We write
applying the operator to $y,z\in\mathbb{R}^{d}$ as $D^{2}u(x)(y,z)$. We also
define the operator norm of $D^{i}u(x)$ by

$\displaystyle\norm{D^{i}u(x)}=\sup\frac{D^{i}u(x)(x_{1},\ldots,x_{i})}{\prod_{j=1}^{i}\norm{x_{j}}}.$

For simplicity, we let $\norm{D^{i}u}=\sup_{x\in K\cup
K_{-r^{\prime}}}\norm{D^{i}u(x)}$. Analogically, we let $D^{0}f$ to be $f$
itself and $\norm{D^{0}f}$ to be sup norm of f, $\norm{f}_{\infty}$.

Then, we discuss the Taylor expansion of $u(X_{k+1})$ of (7) around $X_{k}$.
Following Leimkuhler et al., (2020), we write
$u_{k}:=u(X_{k}),u^{\prime}_{k}:=u(X^{\prime}_{k}),v_{k}^{\pi}:=v(X_{k}^{\pi})$
and $f_{k}:=f(X_{k})$, where $X^{\prime}_{k}$ is introduced in 2.1 and
$X_{k}^{\pi}:=\mathcal{P}_{K}(X^{\prime}_{k})$. Following this notation, we
write

$\displaystyle
u_{k+1}-u_{k}=\left(u_{k+1}-u_{k+1}^{\prime}\right)+\left(u_{k+1}^{\prime}-u_{k}\right).$

For the first term, since $X_{k+1}^{\pi}$ is the midpoint between
$X^{\prime}_{k+1}$ and $X_{k+1}$ by (2), we can write

$\displaystyle X_{k+1}$
$\displaystyle=\mathcal{R}_{K}(X^{\prime}_{k+1})=X^{\prime}_{k+1}-2r_{k+1}v_{k+1}^{\pi}=X^{\pi}_{k+1}-r_{k+1}v_{k+1}^{\pi},$
$\displaystyle X^{\prime}_{k+1}$
$\displaystyle=X^{\prime}_{k+1}-r_{k+1}v_{k+1}^{\pi}+r_{k+1}v_{k+1}^{\pi}=X^{\pi}_{k+1}+r_{k+1}v_{k+1}^{\pi},$

where $r_{k}=\norm{X_{k}^{\pi}-X_{k}}$. Then, we obtain

$\displaystyle u_{k+1}-u_{k+1}^{\prime}=$ $\displaystyle
u(X_{k+1}^{\pi}-r_{k+1}v_{k+1}^{\pi})-u(X_{k+1}^{\pi}+r_{k+1}v_{k+1}^{\pi}).$

Using a Taylor expansion up to third order, we obtain

$\displaystyle\begin{split}&u(X_{k+1}^{\pi}-r_{k+1}v_{k+1}^{\pi})-u(X_{k+1}^{\pi}+r_{k+1}v_{k+1}^{\pi})\\\
&\qquad\qquad=\left(u_{k+1}^{\pi}-r_{k+1}Du_{k+1}^{\pi}\left(v_{k+1}^{\pi}\right)+\frac{r_{k+1}^{2}}{2}D^{2}u_{k+1}^{\pi}\left(v_{k+1}^{\pi},v_{k+1}^{\pi}\right)\right.\\\
&\qquad\qquad\qquad\left.-\frac{r_{k+1}^{3}}{6}D^{3}u(X_{k+1}^{\pi}-\alpha_{1}r_{k+1}v_{k+1}^{\pi})\left(v_{k+1}^{\pi},v_{k+1}^{\pi},v_{k+1}^{\pi}\right)\right)\\\
&\qquad\qquad\qquad-\left(u_{k+1}^{\pi}+r_{k+1}Du_{k+1}^{\pi}\left(v_{k+1}^{\pi}\right)+\frac{r_{k+1}^{2}}{2}D^{2}u_{k+1}^{\pi}\left(v_{k+1}^{\pi},v_{k+1}^{\pi}\right)\right.\\\
&\qquad\qquad\qquad\left.+\frac{r_{k+1}^{3}}{6}D^{3}u(X_{k+1}^{\pi}+\alpha_{2}r_{k+1}v_{k+1}^{\pi})\left(v_{k+1}^{\pi},v_{k+1}^{\pi},v_{k+1}^{\pi}\right)\right),\end{split}$
(9)

where $\alpha_{1},\alpha_{2}\in[0,1]$ are constants. Here, we have used the
mean value theorem in the third derivative term. The second derivative term is
canceled out and thus we have

$\displaystyle\begin{split}u_{k+1}-u_{k+1}^{\prime}=&-\frac{r_{k+1}^{3}}{6}D^{3}u(X_{k+1}^{\pi}-\alpha_{1}r_{k+1}v_{k+1}^{\pi})\left(v_{k+1}^{\pi},v_{k+1}^{\pi},v_{k+1}^{\pi}\right)\\\
&-\frac{r_{k+1}^{3}}{6}D^{3}u(X_{k+1}^{\pi}+\alpha_{2}r_{k+1}v_{k+1}^{\pi})\left(v_{k+1}^{\pi},v_{k+1}^{\pi},v_{k+1}^{\pi}\right)\\\
\geq&-\frac{r_{k+1}^{3}}{3}\norm{D^{3}u},\end{split}$ (10)

since $Du_{k+1}^{\pi}\left(v_{k+1}^{\pi}\right)=\gradient
u(X^{\pi}_{k+1})\cdot v_{k+1}^{\pi}=0$ from (7). The last inequality comes
from the definition of the operator norm and
$\norm{v_{k+1}^{\pi}}\in\left\\{0,1\right\\}$.

In the same way, we obtain

$\displaystyle\begin{split}u^{\prime}_{k+1}-u_{k}=&u\left(X_{k}-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)-u_{k}\\\
=&Du_{k}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\\\
&+\frac{1}{2}D^{2}u_{k}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\\\
&+\frac{1}{6}D^{3}u\left(X_{k}+\alpha_{3}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right)\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},\right.\\\ &\left.-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right),\end{split}$ (11)

where $\alpha_{3}\in[0,1]$ is some constant. Since $\mathbb{E}[\xi_{k+1}]=0$,
taking the expectation yields

$\displaystyle\begin{split}\mathbb{E}\left[u^{\prime}_{k+1}-u_{k}\right]=&\mathbb{E}\left[-\eta\gradient
f_{k}\dot{\gradient}u_{k}+\frac{\eta}{\beta}\Delta
u_{k}+\frac{\eta^{2}}{2}D^{2}u_{k}\left(\gradient f_{k},\gradient
f_{k}\right)\right.\\\
&\left.+\frac{1}{6}D^{3}u\left(X_{k}+\alpha_{3}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right)\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},\right.\right.\\\
&\left.\left.-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right].\end{split}$ (12)

Then, by the definition of $\mathcal{A}$, we obtain

$\displaystyle\begin{split}\mathbb{E}\left[u^{\prime}_{k+1}-u_{k}\right]=&\mathbb{E}\left[\eta\mathcal{A}u_{k}+\frac{\eta^{2}}{2}D^{2}u_{k}\left(\gradient
f_{k},\gradient f_{k}\right)\right.\\\
&\left.+\frac{1}{6}D^{3}u\left(X_{k}+\alpha_{3}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right)\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},\right.\right.\\\
&\left.\left.-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right]\\\
=&\mathbb{E}\left[\eta\left(f_{k}-\mathbb{E}_{\pi}f\right)+\frac{\eta^{2}}{2}D^{2}u(X_{k})\left(\gradient
f_{k},\gradient f_{k}\right)\right.\\\
&\left.+\frac{1}{6}D^{3}u\left(X_{k}+\alpha_{3}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right)\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},\right.\right.\\\
&\left.\left.-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right],\end{split}$ (13)

We apply (10), (13) to

$\displaystyle\mathbb{E}\left[u_{k+1}^{\prime}-u_{k}\right]=\mathbb{E}\left[u_{k+1}-u_{k}\right]-\mathbb{E}\left[u_{k+1}-u_{k+1}^{\prime}\right],$

to obtain

$\displaystyle\mathbb{E}\left[\eta\left(f_{k}-\mathbb{E}_{\pi}f\right)+\frac{\eta^{2}}{2}D^{2}u(X_{k})\left(\gradient
f_{k},\gradient f_{k}\right)\right.$
$\displaystyle\qquad\left.+\frac{1}{6}D^{3}u\left(X_{k}+\alpha_{3}\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right)\left(-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},\right.\right.$
$\displaystyle\qquad\left.\left.-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1},-\eta\gradient
f_{k}+\sqrt{\frac{2\eta}{\beta}}\xi_{k+1}\right)\right]$
$\displaystyle\qquad\qquad\leq\mathbb{E}\left[u_{k+1}-u_{k}\right]+\frac{r_{k+1}^{3}}{3}\norm{D^{3}u}.$

Since $\norm{\xi_{k+1}}=\sqrt{d}$, and $\norm{\gradient f(X_{k})}$ is
uniformly bounded, we have

$\displaystyle\mathbb{E}\left[\eta\left(f_{k}-\mathbb{E}_{\pi}f\right)\right]\leq\mathbb{E}\left[u_{k+1}-u_{k}\right]+\frac{r_{k+1}^{3}}{3}\norm{D^{3}u}+C\eta^{2}\norm{D^{2}u}+C\norm{D^{3}u}\left(\eta+\sqrt{\frac{d\eta}{\beta}}\right)^{3}.$
(14)

To obtain the upper bounds of the right-hand side, we will use the following
lemma, which is proved in Section C.2.

######  Lemma C.1.

There exists a constant $C_{f}>0$ depending only on $\norm{D^{i}f}\>(i=0,2,3)$
such that, for any $\beta>0$,

$\displaystyle\forall x\in K\cup
K_{-r^{\prime}},\quad\norm{D^{i}u(x)}<\frac{C_{f}\beta^{i/2}}{\lambda_{*}}\quad(i=0,1,2,3),$

where $r^{\prime}$ is the supremum of positive numbers such that the
reflection of any $x\in K_{-r^{\prime}}$ is well-defined.

##### Completing the Proof of Lemma 3.7

Therefore, combining (14) with Lemma C.1 yields the following inequality:

$\displaystyle\mathbb{E}\left[\eta\left(f(X_{k})-\mathbb{E}_{\pi}f\right)\right]<\mathbb{E}\left[u(X_{k+1})-u(X_{k})\right]+\frac{C}{\lambda_{*}}\left(\beta^{3/2}r_{k+1}^{3}+\beta\eta^{2}+\beta^{3/2}\eta^{3}+(d\eta)^{3/2}\right).$

Since $r_{k+1}=O(\eta+\sqrt{d\eta/\beta})$, we obtain

$\displaystyle\mathbb{E}\left[\eta\left(f(X_{k})-\mathbb{E}_{\pi}f\right)\right]<\mathbb{E}\left[u(X_{k+1})-u(X_{k})\right]+\frac{C\eta^{3/2}\left(\beta\sqrt{\eta}+(\beta\eta+d\right)^{3/2})}{\lambda_{*}}.$

Thus,

$\displaystyle\frac{1}{N}\sum_{k=1}^{N}\mathbb{E}\left[f(X_{k})-\mathbb{E}_{\pi}f\right]$
$\displaystyle<\mathbb{E}\left[\frac{u(X_{N})-u(X_{0})}{N\eta}\right]+\frac{C\sqrt{\eta}\left(\beta\sqrt{\eta}+(\beta\eta+d)^{3/2}\right)}{\lambda_{*}}$
$\displaystyle<\frac{C}{\lambda_{*}N\eta}+\frac{C\sqrt{\eta}\left(\beta\sqrt{\eta}+(\beta\eta+d)^{3/2}\right)}{\lambda_{*}},$

by Lemma C.1 again.

###  C.2 Proof of Lemma C.1

In this subsection, we prove Lemma C.1. From Benchérif-Madani and Pardoux,
(2009, Theorem 4), we have the following probabilistic representation of $u$:

$\displaystyle
u(x)=\int_{0}^{\infty}\mathbb{E}\left[f(X_{t}^{x})-\mathbb{E}_{\pi}f\right]\differential
t,$

where $X_{t}^{x}$ denotes the solution of (3) at time $t$ starting from $x$.
Let $\phi(t,x):=\mathbb{E}\left[f(X_{t}^{x})-\mathbb{E}_{\pi}f\right]$; we
will prove the following.

######  Lemma C.2.

There exists a constant $C>0$ with the same constant as in Lemma C.1 such
that, for any $\beta>0$, $t\geq 0$, and $x\in K$,

$\displaystyle\norm{D^{i}\phi(t,x)}<C\beta^{i/2}\mathrm{e}^{-\lambda_{*}t}\quad(i=1,2,3),$

where $\lambda_{*}$ is the spectral gap defined in (5).

Integrating over $t$ completes the proof. Thus, it suffices to show Lemma C.2.
First, we state the following lemma using the Bismut-Elworthy-Li formula.

######  Lemma C.3 (Cerrai, (1998, Propositions 6.1, 6.3, 6.4)).

For any $g(x)\in C^{4}(K\cup K_{-r^{\prime}})$, there exists a constant $C>0$
such that for any $\beta>0$, $0<t\leq 1$, and $x\in K\cup K_{-r^{\prime}}$,

$\displaystyle\norm{D^{i}\Phi(t,x)}<C\left(\frac{\beta}{t}\right)^{i/2}\norm{g}_{\infty},\quad(i=1,2,3),$

where $\Phi(t,x)=\mathbb{E}\left[g(X_{t}^{x})\right]$. $C$ depends
polynomially on the exponential of the supremum of
$\norm{D^{2}f},\norm{D^{3}f},\norm{D^{4}f}$.

###### Proof.

Most of the proof is given in Cerrai, (1998). The statement is rectified by
noticing that $N=\sqrt{\beta/2}$. ∎

Then, the convergence rate of $\phi$ is represented by the spectral gap
$\lambda_{*}$.

######  Lemma C.4.

(Freidlin,, 1985, Chapter II) There exists a constant $C$ such that for any
$t\geq 0$,

$\displaystyle\absolutevalue*{\phi(t,x)}<C\mathrm{e}^{-\lambda_{*}t},$

where $\lambda_{*}$ is the spectral gap defined in (5).

Let us see that we can prove Lemma C.2 from Lemmas C.3 and C.4. We define

$\displaystyle S:=\\{t\geq 0\>|\>X_{t}(x)\in\partial
K\\},\;r(t):=\sup(S\cap[0,t]),$

where $x\in K$. As preparation, we define three equations and their solutions.
For $h\in\mathbb{R}^{d}$ and $t\geq 0$, let $Y_{t}^{h}(x)$ be the solution of
the derivative equation,

$\displaystyle\begin{dcases}Y_{t}^{h}=h-\int_{0}^{t}\Delta
f(X_{s}^{x})Y_{s}^{h}\differential s,&\text{if}\>t<\inf S,\\\
Y_{t}^{h}=\mathcal{P}_{X_{r(t)}^{x}}(Y_{r(t)-}^{h})-\int_{r(t)}^{t}\Delta
f(X_{s}^{x})Y_{s}^{h}\differential s,&\text{if}\>t\geq\inf S,\end{dcases}$

where $\mathcal{P}_{x}$ is defined by

$\displaystyle\mathcal{P}_{x}(z)=z-\langle z,v(x)\rangle v(x),\quad\forall
x\in\partial K,\>z\in\mathbb{R}^{d},$

and $Y_{r(t)-}^{h}$ is the left-hand limit of $Y_{t}^{h}$ as $t\uparrow r(t)$.
It is possible to show (see Andres, (2011)) that $Y_{t}^{h}(x)$ is the first-
order derivative in probability of $X_{t}^{x}$ with respect to $x$ along the
direction $h\in\mathbb{R}^{d}$. We let $Z_{t}^{h,k}(x)$ be the solution of the
problem,

$\displaystyle\begin{dcases}Z_{t}^{h,k}=-\int_{0}^{t}\Delta
f(X_{s}^{x})Z_{s}^{h,k}+D^{3}f(X_{s}^{x})\left(Y_{s}^{h},Y_{s}^{k}\right)\differential
s,&\text{if}\>t<\inf S,\\\ Z_{t}^{h,k}=-\int_{r(t)}^{t}\Delta
f(X_{s}^{x})Z_{s}^{h,k}+D^{3}f(X_{s}^{x})\left(Y_{s}^{h},Y_{s}^{k}\right)\differential
s,&\text{if}\>t\geq\inf S.\end{dcases}$

As with the first derivative, it is possible to show (see Cerrai, (1998)) that
$Z_{t}^{h,k}(x)$ is the second-order derivative in probability of $X_{t}^{x}$
with respect to $x$ along the directions $h$ and $k$. Finally, let
$V_{t}^{l,m,n}(x)$ be the solution of the following equation:

$\displaystyle\begin{dcases}V_{t}^{l,m,n}=-\int_{0}^{t}\Delta
f(X_{s}^{x})V_{s}^{l,m,n}+D^{4}f(X_{s}^{x})\left(Y_{s}^{l},Y_{s}^{m},Y_{s}^{n}\right)\\\
\qquad\qquad\qquad\qquad+\frac{1}{2}\sum_{\sigma\in\mathcal{S}_{3}}D^{3}f(X_{s}^{x})\left(Y_{s}^{\sigma(l)},Z_{s}^{\sigma(m),\sigma(n)}\right)\differential
s,&\text{if}\>t<\inf S,\\\ V_{t}^{l,m,n}=-\int_{r(t)}^{t}\Delta
f(X_{s}^{x})V_{s}^{l,m,n}+D^{4}f(X_{s}^{x})\left(Y_{s}^{l},Y_{s}^{m},Y_{s}^{n}\right)\\\
\qquad\qquad\qquad\qquad+\frac{1}{2}\sum_{\sigma\in\mathcal{S}_{3}}D^{3}f(X_{s}^{x})\left(Y_{s}^{\sigma(l)},Z_{s}^{\sigma(m),\sigma(n)}\right)\differential
s,&\text{if}\>t\geq\inf S,\end{dcases}$

where $\mathcal{S}_{3}$ denotes the set of all permutations of three elements.
As with the first two derivatives, $V_{t}^{l,m,n}(x)$ is the third-order
derivative in probability of $X_{t}^{x}$ with respect to $x$ along the
directions $l$, $m$ and $n$.

###### Proof of Lemma C.2.

In the case of $t_{0}\geq 1$, applying Lemma C.3 to $t:=1$,
$g(x):=\phi(t_{0}-1,x)$ and evaluating $\norm{g}_{\infty}$ by Lemma C.4 yield

$\displaystyle\norm{D^{i}\mathbb{E}\left[\phi(t_{0}-1,X_{1}^{x})\right]}<C\beta^{i/2}\mathrm{e}^{-\lambda_{*}(t_{0}-1)}\Leftrightarrow\norm{D^{i}\phi(t_{0},x)}<C\beta^{i/2}\mathrm{e}^{-\lambda_{*}(t_{0}-1)}<C\beta^{i/2}\mathrm{e}^{-\lambda_{*}t_{0}}.$

Replacing $t_{0}$ with $t$ proves the statement. Otherwise, we start by
proving that $Y_{t}$, $Z_{t}$ and $V_{t}$ are bounded for any $0<t\leq 1$.

$\displaystyle\norm{Y_{t}^{h}}$
$\displaystyle\leq\norm{h}+\int_{0}^{t}\norm{\Delta
f(X_{s}^{x})Y_{s}^{h}}\differential s$
$\displaystyle\leq\norm{h}+C\int_{0}^{t}\norm{Y_{s}^{h}}\differential
s,\quad\text{if}\>t<\inf S.$

Then, from Gronwall’s inequality,

$\displaystyle\norm{Y_{t}^{h}}\leq\norm{h}\exp(\int_{0}^{t}C\differential
s)\leq\norm{h}\exp(\int_{0}^{1}C\differential
s)\leq\norm{h}\mathrm{e}^{C}<C\norm{h}.$

If $t\geq\inf S$, by making a recursion argument, we have

$\displaystyle\norm{Y_{t}^{h}}<C\norm{\mathcal{P}_{X_{r(t)}^{x}}(Y_{r(t)-}^{h})}<C\norm{Y_{r(t)-}^{h}}<C\norm{h}.$

As for $Z_{t}$, from the result on $Y_{t}^{h}$, we have

$\displaystyle\norm{Z_{t}^{h,k}}\leq\int_{0}^{t}C\norm{Z_{s}^{h,k}}+C\norm{h}\norm{k}\differential
s\leq C\norm{h}\norm{k}+\int_{0}^{t}C\norm{Z_{s}^{h,k}}\differential s.$

Moreover, from Gronwall’s inequality, we get $\norm{Z_{t}^{h,k}}\leq
C\norm{h}\norm{k}$. By the same argument, it follows that $V_{t}^{l,m,n}\leq
C\norm{l}\norm{m}\norm{n}$. By the way, $D^{i}\phi(t,x)(i=1,2,3)$ can be
calculated in the following way:

$\displaystyle D\phi(t,x)\cdot h=$
$\displaystyle\mathbb{E}\left[Df(X_{t}^{x})\cdot Y_{t}^{h}(x)\right],$
$\displaystyle D^{2}\phi(t,x)\left(h,k\right)=$
$\displaystyle\mathbb{E}\left[D^{2}f(X_{t}^{x})\left(Y_{t}^{h},Y_{t}^{k}\right)+Df(X_{t}^{x})\cdot
Z_{t}^{h,k}\right],$ $\displaystyle D^{3}\phi(t,x)\left(l,m,n\right)=$
$\displaystyle\mathbb{E}\left[D^{3}f(X_{t}^{x})\left(Y_{t}^{l},Y_{t}^{m},Y_{t}^{n}\right)+Df(X_{t}^{x})\cdot
V_{t}^{l,m,n}\right]$
$\displaystyle+\mathbb{E}\left[\frac{1}{2}\sum_{\sigma\in\mathcal{S}_{3}}D^{2}f(X_{t}^{x})\left(Y_{t}^{\sigma(l)},Z_{t}^{\sigma(m),\sigma(n)}\right)\right].$

Combining these formulas with the above statement gives

$\displaystyle\norm{D\phi(t,x)\cdot h}$ $\displaystyle\leq C\norm{h},$
$\displaystyle\norm{D^{2}\phi(t,x)\left(h,k\right)}$ $\displaystyle\leq
C\norm{h}\norm{k},$ $\displaystyle\norm{D^{3}\phi(t,x)\left(l,m,n\right)}$
$\displaystyle\leq C\norm{l}\norm{m}\norm{n}.$

In other words, $\norm{D^{i}\phi(t,x)}$ is bounded by some constant for any
$0<t\leq 1$. Thus, with a change of constant again, we get the desired result.
∎

##  Appendix D Multi-Dimensional Problems

In this section, we demonstrate the performance of RGLD on classical non-
convex test problems. In particular, we consider the Rosenbrock function (15)
and Rastrigin function (16),

$\displaystyle f(x)$
$\displaystyle=\sum_{i=1}^{d-1}\left\\{100(x_{i+1}-x_{i}^{2})^{2}+(1-x_{i})^{2}\right\\},$
(15) $\displaystyle f(x)$
$\displaystyle=10d+\sum_{i=1}^{d}\left\\{x_{i}^{2}-10\cos(2\pi
x_{i})\right\\}.$ (16)

The feasible region of the Rosenbrock function is set to the region between a
sphere of radius $0.5\sqrt{d}$ and one of radius $2\sqrt{d}$ centered at the
origin, and that of the Rastrigin function is set to the region between a
sphere of radius 0.9 and one of radius 5.12 centered at the origin. Figure 9
shows a 3D plot of each function. We can only show a 3D plot when $d=2$, but
the Rosenbrock function has two local minima when $d\geq 4$. It has a global
minimum at $(1,\ldots,1)$, with $f(1,\ldots,1)=0$, and a local minimum around
$(-1,1,\ldots,1)$, with $f(-1,1,\ldots,1)=4$. Meanwhile, the Rastrigin
function has many local optima, which makes global optimization extremely
challenging. It has a global minimum at $(0,\ldots,0)$, with
$f(0,\ldots,0)=0$.

Rosenbrock function

Rastrigin function

Figure 9: 3D plot of the objective function $f$ when $d=2$.

Let us examine how RGLD optimizes these two functions under different
dimensionalities. We scale $\beta$ linearly with $d$ according to our main
theorem and generate the initial values randomly. Figure 10 shows the results
for a constrained Rosenbrock function when $d=4,10,20$. RGLD outperforms PG
only when $d=4$. When $d\in\\{10,20\\}$, RGLD performs worse than PG, which
always stagnates around 4. This is partly because the feasible region gets
larger as $d$ grows and it also implies that the curse of dimensionality
affects RGLD, which may cause an explosion in iteration complexity. Figure 11
shows the results for constrained Rastrigin problems when $d=2,3,5,10,20,30$.
All of the experiments fail to find the global minimum but the figure
demonstrates that RGLD outperforms PG for small dimensionalities except
$d=10$. The failure at $d=10$ is partly because the direction of the
stochastic noise sometimes diverts away from deeper valleys or better local
solutions.

$d=4$

$d=10$

$d=20$

Figure 10: Convergence of the iterates by PG and RGLD for constrained
Rosenbrock problems. $\eta=5\times 10^{-4},\>\beta=0.25\times d$. The initial
values are randomly generated.

$d=2$

$d=3$

$d=5$

$d=10$

$d=20$

$d=30$

Figure 11: Convergence of the iterates by PG and RGLD for constrained
Rastrigin problems. $\eta=5\times 10^{-4},\>\beta=0.05\times d$. The initial
values are randomly generated.

###  D.1 Dependence on Dimensionality

In this subsection, we arrange the plots above to visualize how higher
dimensionality degrades performance. Figure 12 shows that, on the Rosenbrock
problems, RGLD performs worse as $d$ grows. Indeed, RGLD outperforms PG only
when $d=4$. We observe the same tendency on the Rastrigin function (Figure
12).

Rosenbrock

Rastrigin

Figure 12: Convergence of the iterates by RGLD for constrained problems under
different $d$.


# The Effect of Multiple Replies for Natural Language Generation Chatbots

Eason Chen  0000-0003-1486-8559 National Taiwan Normal UniversityTaipeiTaiwan
<EMAIL_ADDRESS>
(2022)

###### Abstract.

In this research, by responding to users’ utterances with multiple replies to
create a group chat atmosphere, we alleviate the problem that Natural Language
Generation chatbots might reply with inappropriate content, thus causing a bad
user experience. Because according to our findings, users tend to pay
attention to appropriate replies and ignore inappropriate replies. We
conducted a 2 (single reply vs. five replies) $\times$ 2 (anonymous avatar vs.
anime avatar) repeated measures experiment to compare the chatting experience
in different conditions. The result shows that users will have a better
chatting experience when receiving multiple replies at once from the NLG model
compared to the single reply. Furthermore, according to the effect size of our
result, to improve the chatting experience for NLG chatbots which is single
reply and anonymous avatar, providing five replies will have more benefits
than setting an anime avatar.

Natural Language Generation, Chatbot, Multiple Replies, Improve Chatting
Experience

††journalyear: 2022††copyright: rightsretained††conference: CHI Conference on
Human Factors in Computing Systems Extended Abstracts; April 29-May 5, 2022;
New Orleans, LA, USA††booktitle: CHI Conference on Human Factors in Computing
Systems Extended Abstracts (CHI ’22 Extended Abstracts), April 29-May 5, 2022,
New Orleans, LA, USA††doi: 10.1145/3491101.3516800††isbn:
978-1-4503-9156-6/22/04††ccs: Human-centered computing Laboratory experiments
Figure 1. The Screenshot of 2 (single reply vs. five replies) $\times$ 2
(anonymous avatar vs. anime avatar) experimental conditions. From left to
right is Five-Replies Anime, Five-Replies Anonymous, Single-Reply Anime,
Single-Reply Anonymous, respectively.

##  1\. Introduction

“Chatbot can be a software which interacts with people using natural language
by building a conversation” (Gunasekara and Vidanage, 2019). Chat-oriented
chatbots are becoming popular in providing chatting services (Gunasekara and
Vidanage, 2019) thanks to the recently developed neural language modeling
(Bengio et al., 2000), especially Natural Language Generation model, like
GPT-2 (Radford et al., 2019).

To generate text, Natural Language Generation (NLG) model adopts a decoding
strategy, such as the top-k or Nucleus method (Holtzman et al., 2020), to
generate one token at a time based on the input text and currently generated
tokens until a certain stop criterion is met. Due to the randomness nature of
the decoding strategy (Holtzman et al., 2020), the generated token sequence
will be different even when the same input sequence is given. Therefore, it is
hard to predict the quality of a single output based on the input.
Consequently, during the conversation between the user and the chatbot,
especially on open-domain topics, the current NLG chatbots inevitably make
errors (Higashinaka et al., 2015) and affect the chatting experience.

We notice that most NLG chatbot researches focus on providing a single reply
which is the best output by the best model (Gunasekara and Vidanage, 2019;
Yang and Tseng, 2020; Rahman et al., 2017). However, since the output from the
NLG model can’t be perfect, we consider another approach to let users maintain
a satisfactory chatting experience: providing more replies. In this research,
we want to know how multiple replies influence the chatting experience for NLG
Chatbot. Therefore, we recruited 102 participants to join a 2 (single reply
vs. five replies) $\times$ 2 (anonymous avatar vs. anime avatar) repeated-
measures experiment and analyzed with two-way ANOVA to compare the chatting
experience in different conditions (see Figure 1).

In sum, the main contributions for this paper are:

  1. (1)

According to our survey, we find users tend to pay attention to appropriate
replies and ignore inappropriate replies when chatting with the chatbot.

  2. (2)

Compared to a chatbot with single reply, the result from our experiment shows
that users have a better chatting experience when the chatbot responds to
users’ one utterance with multiple replies.

  3. (3)

According to the effect size of our result, to improve the chatting experience
for an NLG chatbot which is single reply and anonymous avatar, providing
multiple replies will have more benefits than setting an anime avatar.

##  2\. Theoretical Framework and Hypotheses

When using a chat-oriented chatbot, the key to having a good chatting
experience is to avoid dialog breakdown. Dialog breakdown means users are
feeling difficult to continue their conversation (Higashinaka et al., 2015).
The main reason causing dialog breakdown is inappropriate messages, which
includes ignoring the user’s input, providing a not understandable reply,
repetition, or contradictory content (Higashinaka et al., 2015). Therefore, we
suppose the frequency of suitable and unsuitable messages can predict the
chatting experience. H1: The frequency of appropriate replies can
significantly predict the chatting experience score.

Because the NLG model’s output is highly randomized (Holtzman et al., 2020),
it is hard to maintain the user’s chatting experience. Once users receive an
unacceptable answer, they might stop chatting. However, we can view it
differently: users will continue chatting if they receive an acceptable reply.
So, how to increase the likelihood that users receive a proper response? We
can increase the number of replies. In purely rational analysis, this is a
question of probability: assuming that given the user’s input x, the
probability that the model generates a suitable response will be p(x) (ex:
0.4). Then the probability that the model can’t generate a suitable response
will be 1 – p(x) (ex: 0.6). However, if the model generates five parallel
responses at once, the probability that the model can’t generate any fitting
response will be (1 – p(x))5 (ex: 0.65 = 0.0776). Compared with the single
reply, the chance that the user will not receive any proper reply in the multi
reply condition will be greatly reduced (ex: from 0.6 to 0.0776). In other
words, with more replies, the probability that the model will generate at
least one suitable answer will be greatly increased (ex: from 0.4 to 0.9224).
Thus, users might have a better chatting experience in the multi replies
condition than the single reply condition. This leads to the following
hypothesis: H2: Multiple replies condition has a better user experience than
the single reply condition.

The hypothesis mentioned above be valid or not highly dependent on users’
mental processes. According to the Sensemaking theory (Weick et al., 2005),
people will interpret the information they receive in a way that they can
understand. “Explicit efforts at sensemaking tend to occur when the current
state of the world is perceived to be different from the expected state of the
world” (Weick et al., 2005). Hence, users will be trying to make sense of what
the chatbot wants to say when interacting. Accordingly, suppose we present
more replies to create a group chat atmosphere. Then, it might be easier for
users to comprehend what the chatbot is saying and increase tolerance for
unfitting responses. For example, users might not take those unfitting
messages seriously since they assume the chatbot may not be talking to them
but to other chatbots in the same chatroom. Therefore, we propose the
following hypothesis: H3: Users tend to pay attention to appropriate replies
and ignore inappropriate replies when interacting with chatbots.

Figure 2. The Screenshot of the experimental platform at five-replies anime
condition. The left is the chatroom, and the right is the instruction of the
test.

Nonetheless, when providing multiple replies, we received feedback from pre-
study participants who said that it is hard to distinguish the replies as
‘many from one person’ or ‘many from many people’ without avatar. Therefore,
we provide avatars for each reply as a control condition to help participants
differentiate the multiple replies. We suppose that at the chatting experience
score, there is an interaction effect between the avatar type and the number
of replies. H4: There is an interaction effect between avatar type and the
number of replies at chatting experience score.

##  3\. Research Method

The experiment used a 2 (single reply vs. five replies) $\times$ 2 (anonymous
avatar vs. anime avatar) repeated measures design. Each participant was
required to undergo four randomly ordered experimental conditions. Apart from
that, we choose five replies here because it is best to fit the chatroom
window without overflow while maximizing the replies. And the anime avatars
for chatbots are selected from the popular anime “Attack on Titan” (see in
Figure 1 and 2).

###  3.1. Participants and Procedure

We recruited 102 participants in Taiwan from the internet through convenience
sampling with a mean age of 22.3 years (range 19 – 50, SD = 3.89). At the
experiment, participants will first login to the platform with their own
computer and read the instruction of the experiment. Following the
instruction, participants will have ten times of chat warmup to get familiar
with the platform. After that, participants will go through one of the four
conditions. Each condition participant will first chat ten rounds (Figure 2)
then fill the form for the chat experience of the round. After all four rounds
are finished, the participant will fill the final form for suitable response
frequency and their attention tendency.

###  3.2. Measures

All questionnaires are measured by Google Form.

####  3.2.1. Chatting Experience

We translated The Chatbot Usability Questionnaire (CUQ) (Holmes et al., 2019)
into Chinese to measure the chatting experience in each round. The scale is
designed to measure the user experience of using chatbots. The participant
must select 1 (strongly disagree) to 5 (strongly agree) for sixteen statements
(such as The chatbot understood me well, or Chatbot responses were
irrelevant). According to the author (Holmes et al., 2019), the score of CUQ
has a significant positive correlation with the score from System Usability
Scale (SUS) (Brooke, 1996) and the User Experience Questionnaire (UEQ)
(Laugwitz et al., 2008). Furthermore, based on the data we collected, the
value for Cronbach’s Alpha for the CUQ was $\alpha=.94$.

####  3.2.2. The frequency of appropriate replies

To measure the frequency of appropriate replies, we use a custom five-point
scale (Table 1). Participants must select 1 (strongly disagree) to 5 (strongly
agree) for the following two statements. The frequency of suitable replies
will be obtained by the score of the first question, plus six minus the score
of the second question. These statements were verified with three experts from
relative fields, and based on the data we collected, the value for Cronbach’s
Alpha for this survey was $\alpha=.83$. Participants will review the
definition of appropriate and inappropriate replies before measurement begins.

Table 1. Questions to measure the frequency of suitable replies. Nr. |  Statement
---|---
1 |  In the previous chat, there are many ‘appropriate’ replies
2 |  In the previous chat, there are many ‘inappropriate’ replies

####  3.2.3. Attention/Ignore tendency for appropriate replies and
inappropriate replies

To measure participants’ attention tendency, we use a custom five-point scale
(Table 2). The participants must select 1 (strongly disagree) to 5 (strongly
agree) for the following four statements. These statements were verified with
three experts from relative fields, and based on the data we collected, the
value for Cronbach’s Alpha for this survey was $\alpha=.81$. To compare users’
attention tendency, we will use the paired samples t-test for questions 1, 2,
and 3, 4, respectively.

Table 2. Questions to measure attention/ignore tendency for appropriate replies and inappropriate replies. Nr. |  Statement
---|---
1 |  In the previous chat, I often pay attention to the appropriate replies
2 |  In the previous chat, I often pay attention to inappropriate replies
3 |  In the previous chat, I often ignore appropriate replies
4 |  In the previous chat, I often ignore inappropriate replies

###  3.3. Experimental Platform

The experimental platform is a WEB service with front-end and back-end. The
front-end build using Vue.js and Bootstrap. The back-end uses Python FastAPI
with MongoDB to process the data. We connect Google Sheet API with the
platform so participants can fill the questionnaire at Google Form, and the
platform will check their progress via Sheet API to make sure they finish and
can go to the next round. The platform’s source code is available at
GitHub111EasonC13, AI_Chatbot_experiment_backend:
https://github.com/EasonC13/AI_Chatbot_experiment_backend.

Table 3. The regression result for the frequency of suitable replies predicts CUQ score of the four groups. (Note: *** p ¡ .001). | Five-Replies-Anime | Five-Replies Anonymous | Single-Reply Anime | Single-Reply Anonymous
---|---|---|---|---
The frequency of suitable replies | .576 | 7.046*** | .514 | 5.997*** | .339 | 3.607*** | .360 | 3.860***
F | 49.646*** | 35.969*** |  | 13.014*** |  | 14.897*** |  |
Adjusted R2 | .325 | .257 |  | .106 |  | .121 |  |

###  3.4. NLG Model

About the NLG model, we use the design from Yang and Tseng (Yang and Tseng,
2020), which is a GPT-2 and BERT dual-model reply generation system. To be
specific, after receiving the input message, GPT-2 will first generate many
replies and use BERT to predict replies’ coherence score. Then, the system
will sort the coherence score and use top n replies as the output. The system
can successfully input various dialogues and adjust emotions to output
appropriate responses after learning about the 1.7 million corpus provided by
the Japanese NTCIR Chinese Emotional Conversation Generation (CECG) evaluation
task in 2019 (Zhang and Huang, 2019). To control the variance in the
experiment, we only use ‘like’ as the emotion type and generate one response
per request. Moreover, even though the model is Chinese-based, we connect the
Google Translate API so that the chat system can reply according to the
detected language. That is why some texts in screenshots are English.

##  4\. Result

The data analysis was performed using SPSS 23.0.

###  4.1. Attention/Ignore tendency for appropriate replies and inappropriate
replies

The result from the attention/ignore tendency (Figure 3) showed that the
participants tend to pay attention to the appropriate replies ($M=3.78$,
$\textit{SD}=1.07$) more than inappropriate replies ($M=2.88$,
$\textit{SD}=1.14$) ($t(97)=4.521$, $p<.001$, $d=0.457$) and ignore
inappropriate replies ($M=2.33$, $\textit{SD}=0.95$) more than appropriate
replies ($M=3.40$, $\textit{SD}=1.182$) ($t(97)=-5.890$, $p<.001$, $d=0.595$).
H3 Accepted.

Figure 3. Attention/Ignore tendency survey result for appropriate replies and
inappropriate replies

###  4.2. The relation between the frequency of suitable replies and the CUQ
Score

Four simple linear regressions were carried out to test if the frequency of
suitable replies ($M=3.88$, $\textit{SD}=1.58$) predicted the CUQ score of the
four groups. All four results of the regressions are significant (Table 3).
Indicated that the frequency of suitable replies significantly predicted
chatting experience. Moreover, the result shows that the frequency of suitable
replies can predict the CUQ score more on five-replies conditions. (H1
Accepted)

###  4.3. Comparison of the CUQ scores in the four groups

CUQ scores are subjected to a two-way analysis of variance for the following
four conditions: 2 (single reply vs. five replies) $\times$ 2 (anonymous
avatar vs. anime avatar). The interaction effect for CUQ score between reply
count and avatar type is non-significant ($F(1,97)=2.13$, $p=.148$).
Therefore, H4 can be rejected. Regards to the CUQ scores for the four groups
(Figure 4), from high to low are Five-Replies-Anime ($M=37.5$,
$\textit{SD}=20$), Five-Replies Anonymous ($M=30.9$, $\textit{SD}=18.4$),
Single-Reply Anime ($M=30.5$, $\textit{SD}=18.8$), Single-Reply Anonymous
($M=29.0$, $\textit{SD}=20.1$). Furthermore, the main effect of the number of
replies indicates that for CUQ score, the five replies conditions are
significantly higher than the single reply conditions ($F(1,97)=10.09$,
$p=.002$, $\eta^{2}=.094$). (H2 Accepted) Apart from that, the main effect of
avatar type indicates that for CUQ score, the anime avatar conditions are
significantly higher than the anonymous avatar conditions ($F(1,97)=4.57$,
$p=.035$, $\eta^{2}=.045$).

Figure 4. The difference between the four groups.

##  5\. General Discussion

The result indicates the following: 1\. When users chat with an NLG Chatbot,
they tend to pay attention to appropriate replies and ignore inappropriate
replies. 2\. When designing an NLG chatbot, if responding to a user’s
utterance with multiple replies, simulating a multi-person group chat will
bring a better chatting experience. Therefore, it is easier to satisfy users
with multiple responses, which can reduce the influence of occasional poor
answers of the system. 3\. There is no interaction effect between the avatar
and the number of replies on the chatting experience. Thus, increasing the
number of answers while providing multiple virtual avatars can bring a higher
chatting experience score. Also, according to the effect size ($\eta^{2}$), to
improve the chatting experience for an NLG chatbot which is single reply and
anonymous avatar, providing multiple replies ($\eta^{2}=.094$) might benefit
more than setting an anime avatar ($\eta^{2}=.045$). Our finding demonstrates
the potential of letting chatbot to response to user’s utterance with multiple
replies. At present, most NLG chat-oriented chatbots only provide one reply.
Therefore, we hope that developers can refer to this research and put multiple
replies into design considerations when developing a chat-oriented chatbot in
the future.

###### Acknowledgements.

This work was supported by the Ministry of Science and Technology of Taiwan
(R.O.C.) under Grants 110-2813-C-003-033-E. We thank Yuen-Hsien Tseng, Tzung-
Jin Lin, and Ching-Lin Wu for commenting on our manuscript and suggestion from
Guo-Li Chiou, Chien Wen Tina Yuan, Tsung-Ren Huang and Liang-Yi Li.

## References

  * (1)
  * Bengio et al. (2000) Yoshua Bengio, Réjean Ducharme, and Pascal Vincent. 2000\.  A neural probabilistic language model.  _Advances in Neural Information Processing Systems_ 13 (2000).
  * Brooke (1996) John Brooke. 1996\.  Sus: a “quick and dirty’usability.  _Usability evaluation in industry_ 189, 3 (1996).
  * Gunasekara and Vidanage (2019) Lakindu Gunasekara and Kaneeka Vidanage. 2019.  Uniontbot: Semantic natural language generation based api approach for chatbot communication. In _2019 National Information Technology Conference (NITC)_. IEEE, 1–8.  https://doi.org/10.1109/NITC48475.2019.9114440
  * Higashinaka et al. (2015) Ryuichiro Higashinaka, Masahiro Mizukami, Kotaro Funakoshi, Masahiro Araki, Hiroshi Tsukahara, and Yuka Kobayashi. 2015.  Fatal or not? Finding errors that lead to dialogue breakdowns in chat-oriented dialogue systems. In _Proceedings of the 2015 Conference on Empirical Methods in Natural Language Processing_. 2243–2248.
  * Holmes et al. (2019) Samuel Holmes, Anne Moorhead, Raymond Bond, Huiru Zheng, Vivien Coates, and Michael Mctear. 2019\.  Usability Testing of a Healthcare Chatbot: Can We Use Conventional Methods to Assess Conversational User Interfaces?. In _Proceedings of the 31st European Conference on Cognitive Ergonomics_ (BELFAST, United Kingdom) _(ECCE 2019)_. Association for Computing Machinery, New York, NY, USA, 207–214.  https://doi.org/10.1145/3335082.3335094
  * Holtzman et al. (2020) Ari Holtzman, Jan Buys, Li Du, Maxwell Forbes, and Yejin Choi. 2020.  The curious case of neural text degeneration.  _arXiv preprint arXiv:1904.09751_ (2020).  arXiv:1904.09751 [cs.CL]
  * Laugwitz et al. (2008) Bettina Laugwitz, Theo Held, and Martin Schrepp. 2008\.  Construction and evaluation of a user experience questionnaire. In _Symposium of the Austrian HCI and usability engineering group_. Springer, 63–76.  https://doi.org/10.1007/978-3-540-89350-9_6
  * Radford et al. (2019) Alec Radford, Jeff Wu, Rewon Child, David Luan, Dario Amodei, and Ilya Sutskever. 2019\.  Language Models are Unsupervised Multitask Learners.
  * Rahman et al. (2017) AM Rahman, Abdullah Al Mamun, and Alma Islam. 2017\.  Programming challenges of chatbot: Current and future prospective. In _2017 IEEE Region 10 Humanitarian Technology Conference (R10-HTC)_. IEEE, 75–78.  https://doi.org/10.1109/R10-HTC.2017.8288910
  * Weick et al. (2005) KE Weick, KM Sutcliffe, and D Obstfeld. 2005.  Organizing and the Process of Sensemaking ‘, Organisation Science, 16 (4).
  * Yang and Tseng (2020) Te-Lun Yang and Yuen-Hsien Tseng. 2020.  Development and Evaluation of Emotional Conversation System Based on Automated Text Generation.  _Journal of Educational Media & Library Sciences_ 57, 3 (2020).  https://doi.org/10.6120/JoEMLS.202011_57(3).0048.RS.CM
  * Zhang and Huang (2019) Yaoqin Zhang and Minlie Huang. 2019.  Overview of the NTCIR-14 short text generation subtask: emotion generation challenge. In _Proceedings of the 14th NTCIR Conference_.