# A Proximal Algorithm for Sampling from Non-smooth Potentials

Jiaming Liang Yongxin Chen School of Industrial and Systems Engineering,
Georgia Institute of Technology, Atlanta, GA, 30332. (email:
jiaming.liang@gatech.edu).School of Aerospace Engineering, Georgia Institute
of Technology, Atlanta, GA, 30332. (email: yongchen@gatech.edu).

(Feburary 9, 2021)

###### Abstract

In this work, we examine sampling problems with non-smooth potentials. We
propose a novel Markov chain Monte Carlo algorithm for sampling from non-
smooth potentials. We provide a non-asymptotical analysis of our algorithm and
establish a polynomial-time complexity $\tilde{\cal O}(d\varepsilon^{-1})$ to
obtain $\varepsilon$ total variation distance to the target density, better
than most existing results under the same assumptions. Our method is based on
the proximal bundle method and an alternating sampling framework. This
framework requires the so-called restricted Gaussian oracle, which can be
viewed as a sampling counterpart of the proximal mapping in convex
optimization. One key contribution of this work is a fast algorithm that
realizes the restricted Gaussian oracle for any convex non-smooth potential
with bounded Lipschitz constant.

Key words. High-dimensional sampling, non-smooth potential, complexity
analysis, alternating sampling framework, rejection sampling, proximal bundle
method, restricted Gaussian oracle

##  1 Introduction

Core to many scientific and engineering problems that face uncertainty (either
physically or algorithmically) is the task of drawing samples from a given,
often unnormalized, probability density. Sampling plays a crucial role in many
applications such as statistical inference/estimation, operations research,
physics, biology, and machine learning, etc [3, 12, 14, 16, 18, 19, 20, 40].
For instance, in Bayesian statistics, we can sample from the posterior
distribution to infer its mean, covariance, or other important statistics.
Sampling is heavily used in molecular dynamics to discover new structures.
Sampling is also closely related to optimization. On the one hand,
optimization can be viewed as the limit of sampling when the temperature
parameter goes to 0. On the other hand, sampling can be viewed as an
optimization over the manifold of probability distributions [43, 45].

Over the years, many methods and algorithms have been developed for sampling
[1, 7, 14, 22, 23, 29, 38]. A very popular framework for sampling from high
dimensional complex distributions is the Markov chain Monte Carlo (MCMC)
algorithm [6, 8, 9, 11, 12]. In MCMC, a Markov chain is constructed so that
its invariant distribution is the given target distribution we want to sample
from. After running the Markov chain for sufficiently many iterations, the
state will follow this invariant distribution, generating samples from it.
Several widely used MCMC methods include Langevin Monte Carlo (LMC) [10, 17,
33, 36], Metropolis-adjusted Langevin algorithm (MALA) [4, 35, 36], and
Hamiltonian Monte Carlo (HMC) [32]. These three algorithms use gradient
information of the potential (log-density) to construct the Markov chain. They
resemble the gradient-based algorithms in optimization and can be viewed as
the sampling counterparts of them. Over the last few years, many theoretical
results (see [7, 10, 11, 13, 22, 23, 34, 36] and references therein) have been
established to understand the computational complexities of these MCMC
algorithms.

Most existing gradient-based MCMC methods are only applicable to settings with
smooth potentials [10, 22, 43] whose gradient is Lipschitz continuous.
However, non-smooth sampling is also an important problem as many applications
of sampling involve non-smooth potentials. For instance, in Bayesian
inference, the prior is naturally non-smooth when a compact support is
considered. Many problems in deep learning are also non-smooth, not only due
to non-smooth activation functions like ReLU used in the neural networks, but
also due to intrinsic scaling symmetries. Nevertheless, the study of sampling
without smoothness is nascent. This is in sharp contrast to optimization where
a plethora of algorithms, e.g., subgradient method, proximal algorithm, bundle
method have been developed for non-smooth optimization [25, 26, 27, 28, 30,
37, 44].

Our contributions. The goal of this work is to establish an efficient
algorithm to draw samples from a distribution with non-smooth potential. We
focus on the case where the potential is convex and is Lipschitz continuous.
Our algorithm is based on the recent alternative sampling framework (ASF)
[22], which can be viewed as a sampling counterpart of the proximal point
method in optimization [37]. The key of the alternative sampling framework is
a step known as the restricted Gaussian oracle (RGO) (see Definition 2.1) to
draw samples from a potential regularized by a large isotropic quadratic term.
To utilize this framework to sample from general non-smooth potentials with
bounded Lipschitz constants, we develop an efficient realization of the RGO
through rejection sampling with a properly designed proposal. A non-smooth
optimization technique known as the proximal bundle method [25, 27, 28] is
used to compute the proposal. We establish a polynomial-time complexity
$\tilde{\cal O}(d\varepsilon^{-1})$ to obtain $\varepsilon$ total variation
distance to the target density, better than most existing results under the
same assumptions.

A key contribution of this paper is a fast algorithm for implementing RGO for
any convex (either smooth or non-smooth) function. When the potential $g$ is
decomposable, e.g., $g$ is an $\ell_{1}$ norm or an indicator function of an
orthant, there exists simple sampling algorithms for RGO [31]. In general, the
implementation of RGO is a difficult algorithmic task, which makes ASF [22] a
conceptual method without implementable algorithms in some cases. Our
algorithm of implementing the RGO for any convex function broads the
applicability of the ASF significantly. In fact, our algorithm for RGO can be
used for any framework, not only ASF, that requires sampling from
$\exp(-g(x)-\frac{1}{2\eta}\|x-y\|^{2})$ for any $y$ and some proper $\eta>0$.
From an optimization point of view, our algorithm for RGO provides an
efficient realization of the proximal oracle for a wide range of
functions/potentials, solidifying the connections between the ongoing research
at the interface of optimization and sampling [43].

Related Work. Over the last few years, several new algorithms and theoretical
results in sampling with non-smooth potentials have been established. In [31],
sampling for non-smooth composite potentials is considered. The algorithm
needs the proximal sampling oracle that samples from the target potential
regularized by a large isotropic quadratic term as well as computes the
corresponding partition function, which is not realistic for general
potentials. In [39], algorithms to sample from non-smooth composite potentials
are developed; both are based on the RGO which is similar to the proximal
sampling oracle but do not need to compute the partition function. In [5], the
authors developed an algorithm to sample from non-smooth potentials by running
LMC on the Gaussian smoothing of the potentials. In [24], the author developed
the projected LMC algorithm and analyzed its complexity for non-smooth
potentials. In [15], the authors developed a new analysis that leads to
dimension-free complexity for sampling from a composite density which contains
a non-smooth component. In [11], the authors presented an optimization
approach to analyze the complexity of sampling and established a complexity
result for sampling with non-smooth composite potentials. In [23], the authors
studied the complexity of the ball walk to sample from an isotropic logconcave
density from a warm start. In [2, 42], a proximal algorithm was proposed. This
algorithm resembles the ASF for sampling with a major difference that the RGO
is replaced by proximal point optimization step, which introduces bias for
sampling.

To compare our results with [2] and [15], consider sampling from
$\exp(-f(x)-\mu\|x\|^{2}/2)$ where $f$ is convex and $M$-Lipschitz continuous.
Our complexity (see Theorem 4.1) is $\tilde{\cal O}(M^{2}d/\mu)$, better than
${\cal O}(M^{2}/(\mu\varepsilon^{2}))$ in [15] and $\tilde{\cal
O}(M^{2}d/(\mu\epsilon^{4}))$ (albeit in Wasserstein distance) [2] when
$\varepsilon<d^{-1/2}$. For sampling from non-smooth potentials, compared with
[31, 39], our algorithm does not require any sampling oracle. Compared with
[24], we consider sampling from a distribution supported on $\mathbb{R}^{d}$
instead of a convex compact set. Compared with [5, 11], our algorithm has
better complexity in terms of total variation when the target error
$\varepsilon$ is small. Compared with [23], we consider a generic setting
where the target distribution can be anisotropic, and our complexity is in
general better in the low resolution region. See Table 1 for the detailed
complexity bounds. Note that complexity results obtained in [5, 23] and this
paper are for the last iterate, while the bound established in [11] (also
[15]) is for the average of all iterates.

Paper |  [5] |  [11] |  [23] |  this paper
---|---|---|---|---
Complexity |  $\tilde{\cal O}(M^{6}d^{5}{\cal M}_{4}^{3/2}\varepsilon^{-10})$ |  ${\cal O}(M^{2}W_{2}^{2}\varepsilon^{-4})$ |  ${\cal O}(d^{5/2}\log(\beta/\varepsilon))$ |  $\tilde{\cal O}(M^{2}d{\cal M}_{4}^{1/2}\varepsilon^{-1})$
Table 1: Complexity bounds for sampling from non-smooth densities.

In Table 1, ${\cal M}_{4}$ denotes the finite fourth moment of the target
distribution, $W_{2}$ denotes the Wasserstein distance between the initial and
target distributions, and $\beta$ denotes the warmness of the initial
distribution. More specifically, ${\cal M}_{4}\approx d^{2}$ in the isotropic
case, $\log\beta\approx d$ if the initial distribution is not warm started,
and $W_{2}\approx\sqrt{d}$ in general. Under these simplifications, our bound
is $\tilde{\cal O}(M^{2}d^{2}\varepsilon^{-1})$, better than $\tilde{\cal
O}(M^{6}d^{8}\varepsilon^{-10})$ in [5], ${\cal O}(M^{2}d\varepsilon^{-4})$ in
[11], and $\tilde{\cal O}(d^{7/2})$ in [23] when $d$ is large and $\epsilon$
is relatively small. Note that, for typical problems, $M={\cal O}(d^{1/2})$,
but sparsity maybe exploited to improve this dependence.

Organization. The rest of this paper is structured as follows. In Section 2 we
provide the problem formulation we are interested in. We also briefly review
ASF on which our algorithm is based. In Section 3 we present our key
contribution, an efficient realization of the RGO for general convex
potentials with bounded Lipschitz constants. This is then combined with the
alternating sampling framework in Section 4 to establish our results for
sampling without smoothness. In addition, we further apply the proposed
algorithm to sample from smooth densities in Section 5, and establish
iteration-complexity bounds. Finally, we present some concluding remarks and
possible extensions in Section 6.

##  2 Problem formulation and alternating sampling framework

The problem of interest is to sample from a distribution on $\mathbb{R}^{d}$
proportional to $\exp(-f(x))$ where the potential $f$ is convex and
$M$-Lipschitz continuous. Note that the potential $f$ does not need to be
smooth. This violates the smoothness assumption for most existing gradient-
based MCMC sampling methods [22, 43]. Before describing our approach to design
an efficient algorithm for sampling with non-smooth potentials, we introduce
two important algorithmic notions used in this paper.

Our method is built on the alternating sampling framework (ASF) introduced in
[22] (a similar method was developed in [41]), which is a generic framework
for sampling from a distribution $\exp(-g(x))$; ASF is itself a special case
of Gibbs sampling. For a given point $x\in\mathbb{R}^{d}$ and stepsize
$\eta>0$, the alternating sampling framework repeats the two steps as in
Algorithm 1.

Algorithm 1 Alternating Sampling Framework [22]

0\. sample $y\sim\pi_{x}(y)\propto\exp(-\frac{1}{2\eta}\|x-y\|^{2})$

1\. sample $x\sim\pi_{y}(x)\propto\exp(-g(x)-\frac{1}{2\eta}\|x-y\|^{2})$

In Algorithm 1, sampling $y$ given $x$ in step 1 can be easily done since
$\pi_{x}(y)={\cal N}(x,\eta I)$. Sampling $x$ given $y$ in step 2 corresponds
to the so-called restricted Gaussian oracle for $g$ introduced in [22], which
is the second crucial algorithmic notion used in this paper.

######  Definition 2.1.

Given a point $y\in\mathbb{R}^{d}$ and stepsize $\eta>0$, a restricted
Gaussian oracle (RGO) for convex $g:\mathbb{R}^{d}\to\mathbb{R}$ is a sampling
oracle that returns a sample from a distribution proportional to
$\exp(-g(\cdot)-\|\cdot-y\|^{2}/(2\eta))$.

RGO is an analogy of the proximal mapping in convex optimization, which is
heavily used in proximal point methods. RGO is a key algorithmic ingredient
used in [22] together with the alternating sampling framework to improve the
iteration-complexity bounds for various sampling algorithms. Examples of a
convex function $g$ that admits an computationally efficient RGO have been
presented in [31, 39], including coordinate-separable regularizers,
$\ell_{1}$-norm, and group Lasso.

We recall the main result of [21, 22], which gives the complexity of Algorithm
1 in terms of number of calls to RGO and is useful in this paper.

######  Theorem 2.2.

(Theorem 1 of [21]) Let $\pi$ be a distribution on $\mathbb{R}^{d}$ with
$\pi(x)\propto\exp\left(-f_{\text{oracle }}(x)\right)$ such that
$f_{\text{oracle }}$ is $\mu$-strongly convex, and let $\varepsilon\in(0,1)$.
Let $\eta\leq
1/\mu,T=\Theta\left(\frac{1}{\eta\mu}\log\frac{d}{\eta\mu\varepsilon}\right)$.
Algorithm 1, initialized at the minimizer111Actually, this minimizer can be an
approximate solution, as long as $\|x-x_{\text{opt}}\|^{2}\leq d/\mu$ where
$x$ and $x_{\text{opt}}$ are the approximate and exact solutions,
respectively. of $f_{\text{oracle}}$, runs in $T$ iterations, each querying
RGO for $f_{\text{oracle}}$ with parameter $\eta$ a constant number of times,
and obtains $\varepsilon$ total variation distance to $\pi$.

We are now ready to describe our approach to sample from non-smooth
potentials. We first consider a regularized density $\exp(-g(x))$ where
$g(x)=f(x)+\mu\|x-x^{0}\|^{2}/2$ and $x^{0}\in\mathbb{R}^{d}$ is an arbitrary
point but preferred to be close to the minimum set of $g$. Since $g$ is
$\mu$-strongly convex, the alternating sampling framework is applicable to it.
We then develop an efficient implementation of the RGO based on the proximal
bundle method and rejection sampling for an arbitrary convex potential with
bounded Lipschitz constant. This is the main contribution of this work.
Finally, we justify the sample generated from $\exp(-g(x))$ by using Algorithm
1 and the implementable RGO with a proper choice of $\mu$ is a sample within
$\varepsilon$ total variation distance to the target density $\exp(-f(x))$.

##  3 Key result: an implementable restricted Gaussian oracle

The bottleneck of applying the alternating sampling framework (Algorithm 1) to
sample from general log-concave distributions is the availability of RGO. In
this section, we focus on designing computationally efficient and
implementable RGO for $\mu$-strongly convex $g$ of the form

$g=f+\mu\|\cdot-x^{0}\|^{2}/2.$ (1)

Subsection 3.1 presents an implementation of the RGO for $g$ by using the
proximal mapping of $f$ and rejection sampling. In order to develop an
implementable RGO for $g$ in the cases where there are no efficient
optimization oracles for $f$, we resort to the proximal bundle method, which
is a standard method in convex non-smooth optimization. Subsection 3.2 briefly
reviews the proximal bundle method and its iteration-complexity bound.
Finally, Subsection 3.3 describes the implementation of RGO for $g$ based on
rejection sampling and the proximal bundle method, instead of the proximal
mapping of $f$.

Our algorithm designed for RGO can in fact be used for any convex and
Lipschitz continuous function $f$. More specifically, in both settings where
the proximal mapping of $f$ exists or not, replacing $\mu$ by $0$, the
implementations and results for RGO in this section are also applicable for
$f$.

###  3.1 Sampling with an optimization oracle

Assume that $f$ has a proximal mapping and let

$x^{*}=\underset{x\in\mathbb{R}^{d}}{\mathrm{argmin}\,}\left\\{g^{\eta}(x):=g(x)+\frac{1}{2\eta}\|x-y\|^{2}\right\\}.$
(2)

Here $y$ is the output of ASF (Algorithm 1) in the previous iteration. Note
that solving (2) is equivalent to invoking one proximal mapping of $f$ since
$g(x)=f(x)+\mu\|x-x^{0}\|^{2}/2$; the quadratic term can be combined with that
in (2).

The RGO in each iteration requires sampling from $g^{\eta}$. Our strategy is
rejection sampling with a proper Gaussian proposal centered at $x^{*}$. This
is summarized in Algorithm 2. The following result is useful in the complexity
analysis of Algorithm 2. It is a special case of a more general result (Lemma
3.5) we present later and thus the proof is omitted. Throughout, denote

$\eta_{\mu}:=\eta/(1+\eta\mu).$ (3)

######  Lemma 3.1.

Let

$h_{1}:=\frac{1}{2\eta_{\mu}}\|\cdot-x^{*}\|^{2}+g^{\eta}(x^{*}),\quad
h_{2}:=\frac{1}{2\eta_{\mu}}\|\cdot-x^{*}\|^{2}+2M\|\cdot-x^{*}\|+g^{\eta}(x^{*}).$

Then, for every $x\in\mathbb{R}^{d}$, we have $h_{1}(x)\leq g^{\eta}(x)\leq
h_{2}(x)$.

Algorithm 2 Implementation of the RGO with an optimization oracle

1\. Compute $x^{*}$ as in (2);

2\. Generate $X\sim\exp(-h_{1}(x))$;

3\. Generate $U\sim{\cal U}[0,1]$;

4\. If  $U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))},$ then accept
$\tilde{X}=X$; otherwise, reject $X$ and go to step 2.

The next proposition justifies the correctness and gives the complexity of
Algorithm 2. Its proof is postponed to Appendix B.

######  Proposition 3.2.

Assume $f$ is convex and $M$-Lipschitz continuous. Let
$g=f+\mu\|\cdot-x^{0}\|^{2}/2$ and

$p(x|y)\propto\exp\left(-g(x)-\frac{1}{2\eta}\|x-y\|^{2}\right)$

for a fixed $y$, then $\tilde{X}$ generated by Algorithm 2 is such that
$\tilde{X}\sim p(x|y)$. If $\eta_{\mu}\leq 1/(16M^{2}d)$, then the expected
number of iteration in Algorithm 2 is at most 2.

###  3.2 Review of the proximal bundle method

The proximal bundle method [27, 28] is an efficient algorithm for solving
convex non-smooth optimization problems. In this subsection, we briefly review
an approach to solve the subproblem considered in the proximal bundle method,
the properties of the solution to the subproblem, and the iteration-complexity
for solving the subproblem.

Consider the optimization subproblem (recall (1))

$g^{\eta}_{*}:=g^{\eta}(x^{*})=\min\left\\{g^{\eta}(x):=g(x)+\frac{1}{2\eta}\|x-y\|^{2}:x\in\mathbb{R}^{d}\right\\},$
(4)

and we aim at obtaining a $\delta$-solution (i.e., a point $\bar{x}$ such that
$g^{\eta}(\bar{x})-g^{\eta}_{*}\leq\delta$) to (4). The algorithm is
summarized in Algorithm 3.

We make some remarks about Algorithm 3. First, (6) shows the flexibility in
the choice of $C_{j+1}$. More specifically, choosing
$C_{j+1}=C_{j}\cup\\{x_{j}\\}$ results in the standard cutting-plane model
$f_{j}$ which is underneath $f$, and choosing $C_{j+1}=A_{j}\cup\\{x_{j}\\}$
gives a cutting-plane model $f_{j}$ with less cuts. Second, (5) can be
reformulated into a convex quadratic programming with affine constraints and
the number of constraints is equal to the cardinality of $C_{j}$. Since the
subproblem (5) becomes harder to solve as the size of $C_{j}$ grows, we are in
favor of choosing $C_{j+1}$ in (6) as lean as possible.

The following lemma contains technical results about Algorithm 3 that are
useful in the complexity analysis in Subsection 3.3. Its proof is postponed to
Appendix B.

Algorithm 3 Solving the Proximal Bundle Subproblem (4)

0\. Let $y$, $\eta>0$ and $\delta>0$ be given, and set $\tilde{x}_{0}=y$,
$C_{1}=\\{y\\}$ and $j=1$;

1\. Update $f_{j}=\max\left\\{f(x)+\langle f^{\prime}(x),\cdot-x\rangle:\,x\in
C_{j}\right\\}$;

2\. Define $g_{j}:=f_{j}+\mu\|\cdot-x^{0}\|^{2}/2$ and compute
$x_{j}=\underset{u\in\mathbb{R}^{n}}{\mathrm{argmin}\,}\left\\{g_{j}^{\eta}(u):=g_{j}(u)+\frac{1}{2\eta}\|u-y\|^{2}\right\\},$
(5)
$\tilde{x}_{j}\in\mathrm{Argmin}\,\left\\{g^{\eta}(u):u\in\\{x_{j},\tilde{x}_{j-1}\\}\right\\};$

3\. If $g^{\eta}(\tilde{x}_{j})-g_{j}^{\eta}(x_{j})\leq\delta$, then stop;
else, go to step 4;

4\. Choose $C_{j+1}$ such that  $A_{j}\cup\\{x_{j}\\}\subset C_{j+1}\subset
C_{j}\cup\\{x_{j}\\}$ (6) where $A_{j}:=\left\\{x\in C_{j}:f(x)+\langle
f^{\prime}(x),x_{j}-x\rangle=f_{j}(x_{j})\right\\}$;

5\. Set $j$ $\leftarrow$ $j+1$ and go to step 1.

######  Lemma 3.3.

Assume $f$ is convex and $M$-Lipschitz continuous. Let $j$ denote the last
iteration index, then the following statements hold:

  * a)

$f_{j}\leq f$, $g_{j}\leq g$ and
$g_{j}^{\eta}(x_{j})+\|x-x_{j}\|^{2}/(2\eta_{\mu})\leq g_{j}^{\eta}(x)$ for
every $x\in\mathbb{R}^{d}$;

  * b)

$g^{\eta}(\tilde{x}_{j})-g_{j}^{\eta}(x_{j})\leq\delta$;

  * c)

$\left\|\mu(x_{j}-x^{0})+(x_{j}-y)/\eta\right\|\leq M$;

  * d)

$\|x_{j}-\tilde{x}_{j}\|^{2}\leq 2\eta_{\mu}\delta$.

The following result states the iteration-complexity bound for Algorithm 3 to
obtain a $\delta$-solution to the subproblem (4). We have omitted the proof
since it is relatively technical and beyond the scope of this paper, however,
a complete proof can be found in Section 4 of [28].

######  Proposition 3.4.

Algorithm 3 takes $\tilde{\cal O}(\eta_{\mu}M^{2}/\delta+1)$ iterations to
terminate, and each iteration solves an affinely constrained convex quadratic
programming problem.

###  3.3 Sampling without an optimization oracle

Let $j$ denote the last iteration index of Algorithm 3, i.e., Lemma 3.3 holds
with $j$. Define

$\displaystyle h_{1}$ $\displaystyle:=\frac{1}{2\eta_{\mu}}\|\cdot-
x_{j}\|^{2}+g^{\eta}(\tilde{x}_{j})-\delta,$ (7) $\displaystyle h_{2}$
$\displaystyle:=\frac{1}{2\eta_{\mu}}\|\cdot-\tilde{x}_{j}\|^{2}+\left(2M+\frac{\sqrt{2\delta}}{\sqrt{\eta_{\mu}}}\right)\|\cdot-\tilde{x}_{j}\|+g^{\eta}(\tilde{x}_{j}).$
(8)

Algorithms 4 describes the implementation of RGO for $g$ based on Algorithm 3
and rejection sampling. It differs from Algorithm 2 in that: 1) it uses
Algorithm 3 to compute an approximate solution to (4) instead of calling the
proximal mapping $f$ as in (2).

Algorithm 4 Implementation of the RGO without an optimization oracle

1\. Compute $x_{j}$ and $\tilde{x}_{j}$ as in Algorithm 3;

2\. Generate $X\sim\exp(-h_{1}(x))$;

3\. Generate $U\sim{\cal U}[0,1]$;

4\. If  $U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))},$ then accept
$\tilde{X}=X$; otherwise, reject $X$ and go to step 2.

The following lemma is a counterpart of Lemma 3.1 in the context of RGO
without an optimization oracle and plays an important role in Proposition 3.6.
It reduces to Lemma 3.1 when $\delta=0$.

######  Lemma 3.5.

Assume $f$ is convex and M-Lipschitz continuous. Let
$g=f+\mu\|\cdot-x^{0}\|^{2}/2$ and $g^{\eta}$ be as in (2). Then, for every
$x\in\mathbb{R}^{d}$, we have

$h_{1}(x)\leq g^{\eta}(x)\leq h_{2}(x)$ (9)

where $h_{1}$ and $h_{2}$ are as in (7) and (8), respectively.

Proof: Using Lemma 3.3(a)-(b) and the definition of $g_{j}^{\eta}$, we have

$\displaystyle g(\tilde{x}_{j})-g(x)+\frac{1}{2\eta_{\mu}}\|x-x_{j}\|^{2}$
$\displaystyle\leq
g(\tilde{x}_{j})-g_{j}(x)+\frac{1}{2\eta_{\mu}}\|x-x_{j}\|^{2}$
$\displaystyle\leq
g(\tilde{x}_{j})-g_{j}^{\eta}(x_{j})+\frac{1}{2\eta}\|x-y\|^{2}$
$\displaystyle\leq\delta-\frac{1}{2\eta}\|\tilde{x}_{j}-y\|^{2}+\frac{1}{2\eta}\|x-y\|^{2}.$

The first inequality in (9) holds in view of the definition of $h_{1}$ in (7).
Using the definition of $g^{\eta}$ in (2) and the fact that $f$ is
$M$-Lipschtz, we have

$\displaystyle g^{\eta}(x)-g^{\eta}(\tilde{x}_{j})$
$\displaystyle=f(x)-f(\tilde{x}_{j})+\frac{\mu}{2}\|x-x^{0}\|^{2}-\frac{\mu}{2}\|\tilde{x}_{j}-x^{0}\|^{2}+\frac{1}{2\eta}\|x-y\|^{2}-\frac{1}{2\eta}\|\tilde{x}_{j}-y\|^{2}$
$\displaystyle\leq
M\|x-\tilde{x}_{j}\|+\frac{\mu}{2}\|x-\tilde{x}_{j}\|^{2}+\mu\langle
x-\tilde{x}_{j},\tilde{x}_{j}-x^{0}\rangle+\frac{1}{2\eta}\|x-\tilde{x}_{j}\|^{2}+\frac{1}{\eta}\langle
x-\tilde{x}_{j},\tilde{x}_{j}-y\rangle$
$\displaystyle=M\|x-\tilde{x}_{j}\|+\frac{1}{2\eta_{\mu}}\|x-\tilde{x}_{j}\|^{2}+\mu\langle
x-\tilde{x}_{j},\tilde{x}_{j}-x_{j}+x_{j}-x^{0}\rangle+\frac{1}{\eta}\langle
x-\tilde{x}_{j},\tilde{x}_{j}-x_{j}+x_{j}-y\rangle.$

The above inequality, the Cauchy-Schwarz inequality and Lemma 3.3(c)-(d) imply
that

$\displaystyle g^{\eta}(x)-g^{\eta}(\tilde{x}_{j})$ $\displaystyle\leq
M\|x-\tilde{x}_{j}\|+\frac{1}{2\eta_{\mu}}\|x-\tilde{x}_{j}\|^{2}+\frac{1}{\eta_{\mu}}\|x-\tilde{x}_{j}\|\|\tilde{x}_{j}-x_{j}\|+\|x-\tilde{x}_{j}\|\left\|\mu(x_{j}-x^{0})+\frac{x_{j}-y}{\eta}\right\|$
$\displaystyle\leq
M\|x-\tilde{x}_{j}\|+\frac{1}{2\eta_{\mu}}\|x-\tilde{x}_{j}\|^{2}+\frac{\sqrt{2\delta}}{\sqrt{\eta_{\mu}}}\|x-\tilde{x}_{j}\|+M\|x-\tilde{x}_{j}\|$
$\displaystyle=\left(2M+\frac{\sqrt{2\delta}}{\sqrt{\eta_{\mu}}}\right)\|x-\tilde{x}_{j}\|+\frac{1}{2\eta_{\mu}}\|x-\tilde{x}_{j}\|^{2}.$

It follows from the above inequality and the definition of $h_{2}$ in (8) that
the second inequality in (9) holds.

The next proposition is the main result of this subsection and shows that the
number of rejections in Algorithm 4 is small in expectation. Hence, the
implementation of RGO for $g$ is computationally efficient. Its proof is
postponed to Appendix B.

######  Proposition 3.6.

If

$\eta_{\mu}\leq\frac{1}{64M^{2}d},\quad\delta\leq\frac{1}{32d},$ (10)

then the expected number of iterations in the rejection sampling is at most
$3$.

##  4 Sampling from non-smooth potentials

We now combine our implementation of RGO (Algorithm 4) and the ASF (Algorithm
1) to sample from log-concave probability densities with non-smooth
potentials. This section contains two subsections. Subsection 4.1 presents the
iteration-complexity bound for sampling from $\exp(-g(x))$ where
$g=f+\mu\|\cdot-x^{0}\|^{2}/2$. Based on this, Subsection 4.2 provides the
iteration-complexity for sampling from $\exp(-f(x))$ where $f$ is convex and
Lipschitz continuous. Apart from being a transition step to our final result
for sampling without smoothness, the results in Subsection 4.1 provide an
efficient method to sample from composite potentials of the form
$f+\mu\|\cdot-x^{0}\|^{2}/2$ and may be of independent interest.

###  4.1 Total complexity for strongly convex potential

Using the efficient implementation of RGO for $g$ developed in Section 3 and
the alternating sampling framework Algorithm 1, we are now able to sample from
$\exp(-g(x))$ and establish the complexity for this sampling task. The
following theorem states the iteration-complexity bound for Algorithm 1 using
Algorithm 4 as the RGO to sample from $\exp(-g(x))$. Note that this iteration-
complexity bound is poly-logarithmic in the precision $\varepsilon$ in terms
of total variation. Its proof is postponed to Appendix B.

######  Theorem 4.1.

Let $x^{0}\in\mathbb{R}^{d}$, $\varepsilon>0$, $\delta>0$, $M>0$, $\mu>0$ and
$\eta>0$ satisfying

$\frac{\delta}{M^{2}}\leq\eta\leq\min\left\\{\frac{1}{64M^{2}d},\frac{1}{\mu}\right\\}$
(11)

be given. Let $\pi$ be a distribution on $\mathbb{R}^{d}$ satisfying
$\pi(x)\propto\exp(-g(x))=\exp(-f(x)-\mu\|x-x^{0}\|^{2}/2)$ where $f$ is
convex and $M$-Lipschitz continuous on $\mathbb{R}^{d}$. Consider Algorithm 1
using Algorithm 4 as an RGO for step 1, initialized at the minimizer of $g$,
then the iteration-complexity bound for obtaining $\varepsilon$ total
tolerance to $\pi$ in terms of total variation is

$\tilde{\cal
O}\left(\frac{M^{2}}{\mu\delta}\log\left(\frac{d}{\eta\mu\varepsilon}\right)+1\right),$
(12)

and each iteration queries one subgradient oracle of $f$ and solves a
quadratic programming problem. Moreover, the number of Gaussian distribution
sampling queries in Algorithm 1 is

$\Theta\left(\frac{1}{\eta\mu}\log\left(\frac{d}{\eta\mu\varepsilon}\right)+1\right).$
(13)

It is worth noting that if the proximal mapping of $f$ exists, then the
implementation of RGO only requires one call to the proximal mapping and a
number of rejection sampling. As a result, the total complexity for Algorithm
1 is the same as the complexity in Theorem 2.2, i.e.,
$\Theta\left(\frac{1}{\eta\mu}\log\frac{d}{\eta\mu\varepsilon}\right)$.

###  4.2 Total complexity for convex potential

This subsection studies the main problem of this paper, i.e., sampling from
$\exp(-f(x))$. Building upon Theorem 4.1 for sampling from
$\exp(-g(x))=\exp(-f(x)-\mu\|x-x^{0}\|^{2}/2)$ and a proper choice of $\mu$,
the following theorem establishes the iteration-complexity bound for Algorithm
1 to sample from $\exp(-f(x))$. Its proof is postponed to Appendix B.

######  Theorem 4.2.

Let $\pi$ be a distribution on $\mathbb{R}^{d}$ satisfying
$\pi(x)\propto\exp(-f(x))$ where $f$ is convex and $M$-Lipschitz continuous on
$\mathbb{R}^{d}$. Let $x^{0}\in\mathbb{R}^{d}$ and $\varepsilon>0$ be given
and

$\mu=\frac{\varepsilon}{\sqrt{2}\left(\sqrt{{\cal
M}_{4}}+\|x^{0}-x_{\text{min}}\|^{2}\right)}$ (14)

where ${\cal M}_{4}=\int_{x\in\mathbb{R}^{d}}\|x-x_{\text{min}}\|^{4}d\pi(x)$
and $x_{\text{min}}\in\mathrm{Argmin}\,\\{f(x):x\in\mathbb{R}^{d}\\}$. Choose
$\delta>0$ and $\eta>0$ such that (11) holds and consider Algorithm 1 using
Algorithm 4 as an RGO for step 1, applied to $g=f+\mu\|\cdot-x^{0}\|^{2}/2$,
and initialized at the minimizer of $g$. Then, the iteration-complexity bound
for obtaining $\varepsilon$ total tolerance to $\pi$ is

$\tilde{\cal O}\left(\frac{M^{2}\left(\sqrt{{\cal
M}_{4}}+\|x^{0}-x_{\text{min}}\|^{2}\right)}{\varepsilon\delta}\log\left(\frac{d\left(\sqrt{{\cal
M}_{4}}+\|x^{0}-x_{\text{min}}\|^{2}\right)}{\eta\varepsilon^{2}}\right)+1\right).$
(15)

Finally, we remark that (11) implies that $\delta\leq 1/(64d)$, and if we
choose $\delta=Cd^{-1}$ for some universal constant $C>0$, then the total
complexity (15) for sampling from non-smooth potentials becomes $\tilde{\cal
O}(M^{2}d\sqrt{{\cal M}_{4}}\varepsilon^{-1})$.

######  Remark 4.3.

The strategy we use to sample from a non-smooth potential $f$ by considering a
regularized one $g=f+\mu\|\cdot-x^{0}\|^{2}/2$ first is not the only option.
The reason we do so is that the complexity bound for ASF in [21] requires the
potential to be strongly convex. This convergence result for ASF can, however,
be extended. Following a similar argument as in [21], in particular
Proposition 2 and Lemma 2, one can establish the complexity bound (with
respect to total variation) ${\cal
O}(\frac{1}{\eta\psi^{2}}\log(\beta/\varepsilon))$ for ASF with convex
potential where $\psi$ is the isoperimetry constant of the target distribution
and $\beta$ is a warm start constant. Combining this with our RGO
implementation (Algorithm 4) yields a method to sample from non-smooth
potential $f$ with complexity ${\cal
O}(\frac{M^{2}d}{\psi^{2}}\log(\beta/\varepsilon))$. This is better than
${\cal O}(d^{5/2}\log(\beta/\varepsilon))$ in [23], even in the high accuracy
region, if $M$ scales slower than $\sqrt{d}$, which is typical if sparsity
exists.

##  5 Sampling from smooth potentials

Surprisingly, the exactly same algorithm (Algorithm 4) we developed is
applicable to smooth potentials. In this section, we extend the proximal
sampling algorithm developed for non-smooth sampling to its smooth
counterpart, and establish iteration-complexity results. Subsection 5.1
analyzes the RGO (Algorithm 4) without an optimization oracle for smooth
potentials. Subsection 5.2 provides the total complexity results for sampling
from log-concave probability densities with non-smooth potentials; it is the
counterpart of Section 4 for the smooth case.

We assume in this section that $f$ is convex and $L$-smooth, i.e., for every
$u,v\in\mathbb{R}^{d}$,

$f(u)-f(v)-\langle\nabla f(v),u-v\rangle\leq\frac{L}{2}\|u-v\|^{2}.$ (16)

Considering Algorithm 1 using Algorithm 4 as an RGO for step 1, the goal of
this section is to show that this algorithm, which is originally designed for
non-smooth sampling, is also able to sample from smooth potentials.

###  5.1 Analysis of RGO without an optimization oracle

For simplicity, we focus on the case where we do not have an optimization
oracle for $f$. The analysis of RGO with an optimization oracle can be
obtained by taking $\delta=0$. In such a case, we have
$x_{j}=\tilde{x}_{j}=x^{*}$ (Recall $j$ denotes the last iteration index of
Algorithm 3).

The following lemma is a version of Lemma 3.5 in the smooth case. Its proof is
postponed to Appendix B.

######  Lemma 5.1.

Assume $f$ is convex and $L$-smooth. Let $g=f+\mu\|\cdot-x^{0}\|^{2}/2$ and
$g^{\eta}$ be as in (2). Let $h_{1}$ be as in (7), and define

$h_{2}:=\frac{1}{2\eta_{\mu,L}}\|\cdot-x^{*}\|^{2}+g^{\eta}(x^{*})$ (17)

where

$\eta_{\mu,L}:=\frac{\eta}{1+\eta\mu+\eta L}.$ (18)

Then, for every $x\in\mathbb{R}^{d}$, we have $h_{1}(x)\leq g^{\eta}(x)\leq
h_{2}(x)$.

The next proposition gives the number of rejections in Algorithm 4 and is the
counterpart of Proposition 3.6 in the context of $f$ being $L$-smooth. Its
proof is postponed to Appendix B.

######  Proposition 5.2.

If $\eta_{\mu}\leq 1/(Ld)$, then the expected number of iterations in the
rejection sampling is at most $\exp(1/2+\delta)$.

The following result provides the iteration-complexity for obtaining a
$\delta$-solution by Algorithm 3 in the smooth setting. We omit the proof
since it is relatively technical, however, a complete proof can be found in
Section 4 of [28]. Note that the $\tilde{\cal O}$ notation hides
$\log(\delta^{-1})$.

######  Proposition 5.3.

Algorithm 3 takes $\tilde{\cal O}(\eta_{\mu}L+1)$ iterations to terminate, and
each iteration solves an affinely constrained convex quadratic programming
problem.

###  5.2 Total complexity

We establish two total iteration-complexity bounds to sample from log-concave
probability densities with smooth potentials.

The first theorem states the iteration-complexity bound for sampling from
$\exp(-g(x))$ where $g=f+\mu\|\cdot-x^{0}\|^{2}/2$, and $f$ is convex and
$L$-smooth. Its proof is postponed to Appendix B.

######  Theorem 5.4.

Let $x^{0}\in\mathbb{R}^{d}$, $\varepsilon>0$, $L>0$, $\mu>0$ and $\eta>0$
satisfying $\eta\leq 1/(Ld)$ be given. Let $\pi$ be a distribution on
$\mathbb{R}^{d}$ satisfying
$\pi(x)\propto\exp(-g(x))=\exp(-f(x)-\mu\|x-x^{0}\|^{2}/2)$ where $f$ is
convex and $L$-smooth on $\mathbb{R}^{d}$. Consider Algorithm 1 using
Algorithm 4 as an RGO for step 1, initialized at the minimizer of $g$, then
the iteration-complexity bound for obtaining $\varepsilon$ total tolerance to
$\pi$ in terms of total variation is

$\tilde{\cal
O}\left(\frac{Ld}{\mu}\log\left(\frac{Ld^{2}}{\mu\varepsilon}\right)+1\right),$
(19)

and each iteration queries a gradient $\nabla f$ and solves a quadratic
programming problem.

The next theorem is the main result of this section, which studies the total
iteration-complexity for sampling from $\exp(-f(x))$ where $f$ is convex and
$L$-smooth. We omit the proof since the theorem can be similarly proved by
following the proof of Theorem 4.2.

######  Theorem 5.5.

Let $\pi$ be a distribution on $\mathbb{R}^{d}$ satisfying
$\pi(x)\propto\exp(-f(x))$ where $f$ is convex and $L$-smooth on
$\mathbb{R}^{d}$. Let $x^{0}\in\mathbb{R}^{d}$ and $\varepsilon>0$ be given
and $\mu$ be as in (14). Choose $\eta\leq 1/(Ld)$ and consider Algorithm 1
using Algorithm 4 as an RGO for step 1, applied to
$g=f+\mu\|\cdot-x^{0}\|^{2}/2$, and initialized at the minimizer of $g$. Then,
the iteration-complexity bound for obtaining $\varepsilon$ total tolerance to
$\pi$ is

$\tilde{\cal O}\left(\frac{Ld\left(\sqrt{{\cal
M}_{4}}+\|x^{0}-x_{\text{min}}\|^{2}\right)}{\varepsilon}\log\left(\frac{Ld^{2}\left(\sqrt{{\cal
M}_{4}}+\|x^{0}-x_{\text{min}}\|^{2}\right)}{\varepsilon^{2}}\right)+1\right).$

##  6 Conclusion

This paper presents an algorithm based on the alternating sampling framework
for sampling from non-smooth potentials and establishes a complexity bound
$\tilde{\cal O}(d\varepsilon^{-1})$ to obtain $\varepsilon$ total variation
distance to the target density. Moreover, the algorithm is also applicable to
sample from smooth potentials and has a complexity bound $\tilde{\cal
O}(d\varepsilon^{-1})$. The key contribution of this paper is a
computationally efficient implementation of RGO for any convex (either smooth
or non-smooth) function. One direct extension of the paper is to apply the
proposed algorithm to sample from semi-smooth densities, which include smooth
and non-smooth densities as two extreme cases. Another possible extension of
our analysis in this paper is to consider sampling from composite densities
proportional to $\exp(-f(x)-h(x))$ where $f$ is convex and smooth, and $h$ is
convex and semi-smooth.

## Acknowledgement

This work was supported by NSF under grant 1942523 and 2008513.

## References

  * [1] David Applegate and Ravi Kannan.  Sampling and integration of near log-concave functions.  In Proceedings of the twenty-third annual ACM symposium on Theory of computing, pages 156–163, 1991.
  * [2] E. Bernton.  Langevin monte carlo and jko splitting.  In Conference On Learning Theory, pages 1777–1798. PMLR, 2018.
  * [3] Dimitris Bertsimas and Santosh Vempala.  Solving convex programs by random walks.  Journal of the ACM (JACM), 51(4):540–556, 2004.
  * [4] Nawaf Bou-Rabee and Martin Hairer.  Nonasymptotic mixing of the mala algorithm.  IMA Journal of Numerical Analysis, 33(1):80–110, 2013.
  * [5] Niladri Chatterji, Jelena Diakonikolas, Michael I Jordan, and Peter Bartlett.  Langevin monte carlo without smoothness.  In International Conference on Artificial Intelligence and Statistics, pages 1716–1726. PMLR, 2020.
  * [6] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu.  Fast mcmc sampling algorithms on polytopes.  The Journal of Machine Learning Research, 19(1):2146–2231, 2018\.
  * [7] Yuansi Chen, Raaz Dwivedi, Martin J Wainwright, and Bin Yu.  Fast mixing of metropolized hamiltonian monte carlo: Benefits of multi-step gradients.  J. Mach. Learn. Res., 21:92–1, 2020.
  * [8] Xiang Cheng and Peter Bartlett.  Convergence of langevin mcmc in kl-divergence.  In Algorithmic Learning Theory, pages 186–211. PMLR, 2018.
  * [9] Xiang Cheng, Niladri S Chatterji, Peter L Bartlett, and Michael I Jordan.  Underdamped langevin mcmc: A non-asymptotic analysis.  In Conference on learning theory, pages 300–323. PMLR, 2018.
  * [10] Arnak S Dalalyan.  Theoretical guarantees for approximate sampling from smooth and log-concave densities.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 79(3):651–676, 2017.
  * [11] Alain Durmus, Szymon Majewski, and Błażej Miasojedow.  Analysis of langevin monte carlo via convex optimization.  The Journal of Machine Learning Research, 20(1):2666–2711, 2019\.
  * [12] Alain Durmus, Eric Moulines, and Marcelo Pereyra.  Efficient bayesian computation by proximal Markov Chain Monte Carlo: when Langevin meets Moreau.  SIAM Journal on Imaging Sciences, 11(1):473–506, 2018.
  * [13] Raaz Dwivedi, Yuansi Chen, Martin J Wainwright, and Bin Yu.  Log-concave sampling: Metropolis-hastings algorithms are fast!  In Conference on learning theory, pages 793–797. PMLR, 2018.
  * [14] Martin Dyer, Alan Frieze, and Ravi Kannan.  A random polynomial-time algorithm for approximating the volume of convex bodies.  Journal of the ACM (JACM), 38(1):1–17, 1991.
  * [15] Yoav Freund, Yi-An Ma, and Tong Zhang.  When is the convergence time of langevin algorithms dimension independent? a composite optimization viewpoint.  arXiv preprint arXiv:2110.01827, 2021.
  * [16] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin.  Bayesian Data Analysis.  CRC press, 2013.
  * [17] Ulf Grenander and Michael I Miller.  Representations of knowledge in complex systems.  Journal of the Royal Statistical Society: Series B (Methodological), 56(4):549–581, 1994.
  * [18] Adam Tauman Kalai and Santosh Vempala.  Simulated annealing for convex optimization.  Mathematics of Operations Research, 31(2):253–266, 2006.
  * [19] Ravi Kannan, László Lovász, and Miklós Simonovits.  Random walks and an ${O}^{*}(n^{5})$ volume algorithm for convex bodies.  Random Structures & Algorithms, 11(1):1–50, 1997.
  * [20] Werner Krauth.  Statistical mechanics: algorithms and computations, volume 13.  OUP Oxford, 2006.
  * [21] Yin Tat Lee, Ruoqi Shen, and Kevin Tian.  Structured logconcave sampling with a restricted gaussian oracle.  arXiv preprint arXiv:2010.03106, 2020.
  * [22] Yin Tat Lee, Ruoqi Shen, and Kevin Tian.  Structured logconcave sampling with a restricted gaussian oracle.  In Conference on Learning Theory, pages 2993–3050. PMLR, 2021.
  * [23] Yin Tat Lee and Santosh Srinivas Vempala.  Eldan’s stochastic localization and the kls hyperplane conjecture: an improved lower bound for expansion.  In 2017 IEEE 58th Annual Symposium on Foundations of Computer Science (FOCS), pages 998–1007. IEEE, 2017.
  * [24] Joseph Lehec.  The langevin monte carlo algorithm in the non-smooth log-concave case.  Available on arXiv:2101.10695, 2021.
  * [25] Claude Lemaréchal.  An extension of davidon methods to non differentiable problems.  In Nondifferentiable optimization, pages 95–109. Springer, 1975\.
  * [26] Claude Lemaréchal.  Nonsmooth optimization and descent methods.  1978\.
  * [27] Jiaming Liang and Renato D. C. Monteiro.  A proximal bundle variant with optimal iteration-complexity for a large range of prox stepsizes.  SIAM Journal on Optimization, 31(4):2955–2986, 2021.
  * [28] Jiaming Liang and Renato D. C. Monteiro.  A unified analysis of a class of proximal bundle methods for hybrid convex composite optimization problems.  Available on arXiv:2110.01084, 2021.
  * [29] László Lovász and Santosh Vempala.  Fast algorithms for logconcave functions: Sampling, rounding, integration and optimization.  In 2006 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS’06), pages 57–68. IEEE, 2006.
  * [30] Robert Mifflin.  A modification and an extension of Lemaréchal’s algorithm for nonsmooth minimization.  In Nondifferential and variational techniques in optimization, pages 77–90. Springer, 1982.
  * [31] Wenlong Mou, Nicolas Flammarion, Martin J Wainwright, and Peter L Bartlett.  An efficient sampling algorithm for non-smooth composite potentials.  Available on arXiv:1910.00551, 2019.
  * [32] Radford M Neal.  MCMC using hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 2(11):2, 2011.
  * [33] Giorgio Parisi.  Correlation functions and computer simulations.  Nuclear Physics B, 180(3):378–384, 1981.
  * [34] Gareth O Roberts and Jeffrey S Rosenthal.  Optimal scaling of discrete approximations to langevin diffusions.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 60(1):255–268, 1998.
  * [35] Gareth O Roberts and Osnat Stramer.  Langevin diffusions and metropolis-hastings algorithms.  Methodology and Computing in Applied Probability, 4(4):337–357, 2002.
  * [36] Gareth O Roberts and Richard L Tweedie.  Exponential convergence of langevin distributions and their discrete approximations.  Bernoulli, pages 341–363, 1996.
  * [37] Ralph Tyrrell Rockafellar.  Monotone operators and the proximal point algorithm.  SIAM Journal on Control and Optimization, 14(5):877–898, 1976.
  * [38] Ruoqi Shen and Yin Tat Lee.  The randomized midpoint method for log-concave sampling.  In Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [39] Ruoqi Shen, Kevin Tian, and Yin Tat Lee.  Composite logconcave sampling with a restricted gaussian oracle.  Available on arXiv:2006.05976, 2020.
  * [40] Jack W Sites Jr and Jonathon C Marshall.  Delimiting species: a renaissance issue in systematic biology.  Trends in Ecology & Evolution, 18(9):462–470, 2003.
  * [41] Maxime Vono, Daniel Paulin, and Arnaud Doucet.  Efficient MCMC sampling with dimension-free convergence rate using ADMM-type splitting.  Journal of Machine Learning Research, 23(25):1–69, 2022.
  * [42] A. Wibisono.  Proximal langevin algorithm: Rapid convergence under isoperimetry.  arXiv preprint arXiv:1911.01469, 2019.
  * [43] Andre Wibisono.  Sampling as optimization in the space of measures: The langevin dynamics as a composite optimization problem.  In Conference on Learning Theory, pages 2093–3027. PMLR, 2018.
  * [44] Philip Wolfe.  A method of conjugate subgradients for minimizing nondifferentiable functions.  In Nondifferentiable optimization, pages 145–173. Springer, 1975\.
  * [45] Zhuoran Yang, Yufeng Zhang, Yongxin Chen, and Zhaoran Wang.  Variational transport: A convergent particle-based algorithm for distributional optimization.  Available on arXiv:2012.11554, 2020.

##  Appendix A Technical results

######  Lemma A.1.

Useful Gaussian integrals:

  * a)

for any ${\lambda}>0$,

$\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2{\lambda}}\|x\|^{2}\right)dx=(2\pi{\lambda})^{d/2};$

  * b)

for any $c>0$ and $n\geq 1$,

$\int_{0}^{\infty}\exp(-cx^{2})x^{n}dx=\left\\{\begin{array}[]{ll}\frac{(n-1)!!}{2^{n/2+1}c^{n/2}}\sqrt{\frac{\pi}{c}},&\text{
for }n\text{ even},\\\\[4.26773pt]
\frac{\left(\frac{n-1}{2}\right)!}{2c^{(n+1)/2}},&\text{ for }n\text{
odd}.\end{array}\right.$

######  Lemma A.2.

Facts about the Gamma function $\Gamma$: for every $k\geq 1$,

$\Gamma(k)=(k-1)!,\quad\Gamma\left(k+\frac{1}{2}\right)=\left(k-\frac{1}{2}\right)\left(k-\frac{3}{2}\right)\cdots\frac{1}{2}\sqrt{\pi},$
(20)

and

$\sqrt{k}<\frac{\Gamma(k+1)}{\Gamma(k+\frac{1}{2})}<\sqrt{k+\frac{1}{2}}.$
(21)

Proof: Identities in (20) are well-known and hence their proofs are omitted.
We give an elementary proof of (21). Let

$I_{k}=\int_{0}^{\pi/2}\sin^{k}xdx,\quad\forall k\geq 0,$

then integration by parts gives the recursive formula

$I_{k}=\frac{k-1}{k}I_{k-2},\quad\forall k\geq 2.$ (22)

Applying the above identity recursively and using the facts that $I_{0}=\pi/2$
and $I_{1}=1$, we have for every $k\geq 1$,

$I_{2k}=\frac{2k-1}{2k}\frac{2k-3}{2k-2}\cdots\frac{1}{2}\frac{\pi}{2},\quad
I_{2k+1}=\frac{2k}{2k+1}\frac{2k-2}{2k-1}\cdots\frac{2}{3},$

and hence,

$\frac{I_{2k}}{I_{2k+1}}=\frac{2k+1}{2}\frac{((2k)!)^{2}}{4^{2k}(k!)^{4}}\pi.$
(23)

It follows from the fact that $\sin x<1$ for $x\in(0,\pi/2)$ that
$I_{k+1}<I_{k}$ for every $k\geq 0$. This observation together with (22)
implies that

$1<\frac{I_{2k}}{I_{2k+1}}<\frac{I_{2k-1}}{I_{2k+1}}=\frac{2k+1}{2k}.$

Using the above inequality and (23), we have

$\frac{1}{\sqrt{k+\frac{1}{2}}}<\frac{(2k)!}{4^{k}(k!)^{2}}\sqrt{\pi}=\frac{1}{4^{k}}{2k\choose
k}\sqrt{\pi}<\frac{1}{\sqrt{k}}.$ (24)

Finally, it follows from (20) that

$\frac{\Gamma\left(k+\frac{1}{2}\right)}{\Gamma(k+1)}=\frac{1}{4^{k}}{2k\choose
k}\sqrt{\pi},$

which together with (24) implies (21).

######  Proposition A.3.

For $a\geq 0$ and $d\geq 1$, if ${\lambda}\leq 1/(16a^{2}d)$, then

$\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2{\lambda}}\|x\|^{2}-2a\|x\|\right)dx\geq\frac{(2\pi{\lambda})^{d/2}}{2}.$
(25)

Proof: Let $r=\|x\|$ and note that

$dx=r^{d-1}dS^{d-1}dr$

where $S^{d-1}$ is the surface area of the $(d-1)$-dimensional unit sphere.
Define

$F_{d,{\lambda}}(a):=\int_{0}^{\infty}\exp\left(-\frac{1}{2{\lambda}}r^{2}-2ar\right)r^{d}dr,$
(26)

then we have

$\displaystyle\int\exp\left(-\frac{1}{2{\lambda}}\|x\|^{2}-2a\|x\|\right)dx$
$\displaystyle=\int\exp\left(-\frac{1}{2{\lambda}}r^{2}-2ar\right)r^{d-1}drdS^{d-1}$
$\displaystyle=\frac{2\pi^{d/2}}{\Gamma\left(\frac{d}{2}\right)}F_{d-1,{\lambda}}(a).$
(27)

In the above, we have used the fact that the total surface area of a
$(d-1)$-dimensional unit sphere is
$2\pi^{d/2}/\Gamma\left(\frac{d}{2}\right)$. It follows from the definition of
$F_{d,{\lambda}}$ in (26) that

$\frac{dF_{d-1,{\lambda}}(a)}{da}=\int_{0}^{\infty}\exp\left(-\frac{1}{2{\lambda}}r^{2}-2ar\right)(-2r)r^{d-1}dr=-2F_{d,{\lambda}}(a)\geq-2F_{d,{\lambda}}(0),$

and hence that

$F_{d-1,{\lambda}}(a)\geq F_{d-1,{\lambda}}(0)-2aF_{d,{\lambda}}(0).$ (28)

We now consider two cases: $d=2k+1$ for $k\geq 0$ and $d=2k$ for $k\geq 1$.

Case 1: $d=2k+1$. It follows from the definition of $F_{d,{\lambda}}$ in (26)
and Lemma A.1(b) that

$\displaystyle
F_{2k,{\lambda}}(0)=\int_{0}^{\infty}\exp\left(-\frac{1}{2{\lambda}}r^{2}\right)r^{2k}dr=\frac{(2k-1)!!{\lambda}^{k}}{2}\sqrt{2{\lambda}\pi},$
(29) $\displaystyle
F_{2k+1,{\lambda}}(0)=\int_{0}^{\infty}\exp\left(-\frac{1}{2{\lambda}}r^{2}\right)r^{2k+1}dr=\frac{k!(2{\lambda})^{k+1}}{2}.$

Using the above two identities and (28) with $d=2k+1$, we have

$\displaystyle F_{2k,{\lambda}}(a)$ $\displaystyle\geq
F_{2k,{\lambda}}(0)-2aF_{2k+1,{\lambda}}(0)$
$\displaystyle=\frac{(2k-1)!!{\lambda}^{k}}{2}\sqrt{2{\lambda}\pi}-ak!(2{\lambda})^{k+1}$
$\displaystyle={\lambda}^{k}(2k-1)!!\left(\sqrt{\frac{{\lambda}\pi}{2}}-2a{\lambda}\frac{(2k)!!}{(2k-1)!!}\right).$

The above inequality, the fact that

$\frac{(2k)!!}{(2k-1)!!}=\frac{\Pi_{i=1}^{k}2i}{\Pi_{i=1}^{k}(2i-1)}=\frac{\Pi_{i=1}^{k}i}{\Pi_{i=1}^{k}(i-\frac{1}{2})}=\frac{\Gamma(k+1)}{\Gamma(k+\frac{1}{2})}\sqrt{\pi}$

and (21) imply that

$\displaystyle F_{2k,{\lambda}}(a)$
$\displaystyle\geq{\lambda}^{k}(2k-1)!!\left(\sqrt{\frac{{\lambda}\pi}{2}}-2a{\lambda}\sqrt{\left(k+\frac{1}{2}\right)\pi}\right)$
$\displaystyle={\lambda}^{k}(2k-1)!!\left(\sqrt{\frac{{\lambda}\pi}{2}}-a{\lambda}\sqrt{2d\pi}\right)$
$\displaystyle\geq{\lambda}^{k}(2k-1)!!\frac{1}{2}\sqrt{\frac{{\lambda}\pi}{2}}$
$\displaystyle=(2{\lambda})^{k}\Gamma\left(k+\frac{1}{2}\right)\frac{\sqrt{{\lambda}}}{2\sqrt{2}}=\frac{(2{\lambda})^{d/2}\Gamma(\frac{d}{2})}{4}$

where the second inequality is due to the assumption that ${\lambda}\leq
1/(16a^{2}d)$, and the first identity is due to the fact that $d=2k+1$. Using
the above inequality, (27) and the fact that $d=2k+1$, we conclude that (25)
holds for case 1.

Case 2: $d=2k$. It follows from the definition of $F_{d,{\lambda}}$ in (26)
and Lemma A.1(b) that

$F_{2k-1,{\lambda}}(0)=\int_{0}^{\infty}\exp\left(-\frac{1}{2{\lambda}}r^{2}\right)r^{2k-1}dr=\frac{(k-1)!(2{\lambda})^{k}}{2}=(2k-2)!!{\lambda}^{k}.$

Using the above identity, (29), and (28) with $d=2k$, we have

$\displaystyle F_{2k-1,{\lambda}}(a)$ $\displaystyle\geq
F_{2k-1,{\lambda}}(0)-2aF_{2k,{\lambda}}(0)$
$\displaystyle=(2k-2)!!{\lambda}^{k}\left(1-a\sqrt{2{\lambda}\pi}\frac{(2k-1)!!}{(2k-2)!!}\right)$
$\displaystyle=(2k-2)!!{\lambda}^{k}\left(1-2ak\sqrt{2{\lambda}}\frac{\Gamma(k+\frac{1}{2})}{\Gamma(k+1)}\right)$
$\displaystyle\geq(2k-2)!!{\lambda}^{k}\left(1-2a\sqrt{2{\lambda}k}\right)$

where the second identity is due to (20), and the second inequality is due to
(21). It follows from the above inequality, the fact that $d=2k$ and the
assumption that ${\lambda}\leq 1/(16a^{2}d)$ that

$\displaystyle F_{d-1,{\lambda}}(a)$
$\displaystyle\geq(2k-2)!!{\lambda}^{k}\left(1-2a\sqrt{{\lambda}d}\right)\geq\frac{(2k-2)!!{\lambda}^{k}}{2}$
$\displaystyle=\frac{(k-1)!!(2{\lambda})^{k}}{4}=\frac{\Gamma(k)(2{\lambda})^{k}}{4}=\frac{\Gamma(\frac{d}{2})(2{\lambda})^{d/2}}{4}$

where the second identity is due to (20). Using the above inequality, (27) and
the fact that $d=2k$, we conclude that (25) holds for case 2.

##  Appendix B Missing proofs

Proof of Proposition 3.2: It is a well-known result for rejection sampling
that $\tilde{X}\sim p(x)$ and the probability that $\tilde{X}$ is accepted is

$\mathbb{P}\left(U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))}\right)=\frac{\int\exp(-g^{\eta}(x))dx}{\int\exp(-h_{1}(x))dx}.$

Using the assumption that $\eta_{\mu}\leq 1/(16M^{2}d)$, Lemma 3.1 and
Proposition A.3 that

$\int_{\mathbb{R}^{d}}\exp(-g^{\eta}(x))dx\geq\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))dx\geq\exp(-g^{\eta}(x^{*}))\frac{(2\pi\eta_{\mu})^{d/2}}{2}.$

Moreover, it follows from the definition of $h_{1}$ and Lemma A.1(a) that

$\int_{\mathbb{R}^{d}}\exp(-h_{1}(x))dx=\exp(-g^{\eta}(x^{*}))(2\pi\eta_{\mu})^{d/2}.$

The above three relations immediately imply that

$\mathbb{P}\left(U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))}\right)\geq\frac{1}{2},$

and hence that the expected number of iterations is bounded above by 2.

Proof of Lemma 3.3: a) The first two inequalities directly follow from the
convexity of $f$, and the definitions of $g$ and $g_{j}$. The third inequality
follows from (5).

b) This statement immediately follows from step 3 of Algorithm 3.

c) It follows from the optimality condition of (5) and the definition of
$g_{j}$ that

$-\mu(x_{j}-x^{0})-\frac{x_{j}-y}{\eta}\in\partial f_{j}(x_{j}).$

This inclusion and the fact that $\|f^{\prime}(x)\|\leq M$ imply that c)
holds.

d) The last inequality in (a) with $x=\tilde{x}_{j}$ and (b) imply this
statement.

Proof of Proposition 3.6: We first observe that the assumption (10) implies
that

$2\sqrt{\eta_{\mu}}M+\sqrt{2\delta}\leq\frac{1}{2\sqrt{d}},$

which satisfies the assumption in Proposition A.3 with

${\lambda}=\eta_{\mu},\quad a=M+\frac{\sqrt{\delta}}{\sqrt{2\eta_{\mu}}}.$

Using the definition of $h_{2}$ in (8) and Lemma A.3, we have

$\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))dx\geq\frac{1}{2}\exp(-g^{\eta}(\tilde{x}_{j}))(2\pi\eta_{\mu})^{d/2}.$

It follows from the definition of $h_{1}$ in (7) and Lemma A.1(a) with
${\lambda}=\eta_{\mu}$ that

$\int\exp(-h_{1}(x))dx=\exp\left(-g^{\eta}(\tilde{x}_{j})+\delta\right)(2\pi\eta_{\mu})^{d/2}.$

We conclude that

$\displaystyle\mathbb{P}\left(U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))}\right)$
$\displaystyle=\frac{\int\exp(-g^{\eta}(x))dx}{\int\exp(-h_{1}(x))dx}$
$\displaystyle\geq\frac{\int\exp(-h_{2}(x))dx}{\exp\left(-g^{\eta}(\tilde{x}_{j})+\delta\right)(2\pi\eta_{\mu})^{d/2}}\geq\frac{1}{2}\exp(-\delta),$

and the expected number of the iterations is

$\frac{1}{\mathbb{P}\left(U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-f(X))}\right)}\leq
2\exp(\delta)\leq 2(1+2\delta)\leq 2\left(1+\frac{1}{16d}\right)\leq 3$

where the last two inequalities are due to the second inequality in (10).

Proof of Theorem 4.1: It follows from Theorem 2.2 that the iteration-
complexity for Algorithm 1 to obtain $\varepsilon$ total tolerance to $\pi$ is
(13), which together with Proposition 3.4 implies that the total iteration-
complexity is

$\tilde{\cal
O}\left(\left[\frac{\eta_{\mu}M^{2}}{\delta}+1\right]\left[\frac{1}{\eta\mu}\log\left(\frac{d}{\eta\mu\varepsilon}\right)+1\right]\right).$
(30)

Let

$a=\frac{\eta_{\mu}M^{2}}{\delta},\quad b=\frac{1}{\eta\mu}.$

In view of the above definitions of $a$ and $b$, the total iteration-
complexity (30) becomes $\tilde{\cal O}((a+1)(b+1))$. It is easy to see from
(11) that $a\geq 1/2$ and $b\geq 1$, and hence $\tilde{\cal O}((a+1)(b+1))$ is
equal to $\tilde{\cal O}(ab+1)$. Since $\eta_{\mu}\leq\eta$, the total
iteration-complexity $\tilde{\cal O}(ab+1)$ is (12). Moreover, it follows from
(11) that (10) is satisfied, and hence that Proposition 3.6 holds. This
conclusion together with the iteration-complexity (13) for Algorithm 1 implies
the last conclusion of the theorem.

Proof of Theorem 4.2: Let $\rho$ denote the distribution of the points
generated by Algorithm 1 using Algorithm 4 as an RGO, and let $\hat{\pi}$
denote the distribution proportional to $\exp(-g(x))$ Following the proof of
Corollary 4.1 of [5], we similarly have

$\|\rho-\pi\|_{\text{TV}}\leq\|\rho-\hat{\pi}\|_{\text{TV}}+\|\hat{\pi}-\pi\|_{\text{TV}}$

and

$\displaystyle\|\hat{\pi}$
$\displaystyle-\pi\|_{\text{TV}}\leq\frac{1}{2}\left(\int_{\mathbb{R}^{d}}[f(x)-g(x)]^{2}d\pi(x)\right)^{1/2}=\frac{1}{2}\left(\int_{\mathbb{R}^{d}}\left(\frac{\mu}{2}\|x-x^{0}\|^{2}\right)^{2}d\pi(x)\right)^{1/2}$
$\displaystyle\leq\frac{\mu}{2}\left(\int_{\mathbb{R}^{d}}\left(\|x-x_{\text{min}}\|^{2}+\|x_{\text{min}}-x^{0}\|^{2}\right)^{2}d\pi(x)\right)^{1/2}$
$\displaystyle\leq\frac{\mu}{2}\left(\int_{\mathbb{R}^{d}}\left(2\|x-x_{\text{min}}\|^{4}+2\|x_{\text{min}}-x^{0}\|^{4}\right)d\pi(x)\right)^{1/2}$
$\displaystyle=\frac{\sqrt{2}\mu}{2}\left({\cal
M}_{4}+\|x_{\text{min}}-x^{0}\|^{4}\right)^{1/2}\leq\frac{\sqrt{2}\mu}{2}\left(\sqrt{{\cal
M}_{4}}+\|x_{0}-x_{\text{min}}\|^{2}\right)=\frac{\varepsilon}{2}$

where the last identity is due to the definition of $\mu$ in (14). Hence, it
suffices to derive the iteration-complexity bound for Algorithm 1 to obtain
$\|\rho-\hat{\pi}\|_{\text{TV}}\leq\varepsilon/2$, which is (15) in view of
Theorem 4.1 with $\mu$ as in (14).

Proof of Lemma 5.1: The first inequality $h_{1}(x)\leq g^{\eta}(x)$
immediately follows Lemma 3.5. We observe that the optimality condition of (4)
is

$\mu(x^{*}-x^{0})+\frac{1}{\eta}(x^{*}-y)=-\nabla f(x^{*}).$ (31)

It follows from the definition of $g^{\eta}$ in (2) and the observation (31)
that

$\displaystyle g^{\eta}(x)-g^{\eta}(x^{*})$ $\displaystyle=$ $\displaystyle
f(x)-f(x^{*})+\frac{\mu}{2}\|x-x^{0}\|^{2}-\frac{\mu}{2}\|x^{*}-x^{0}\|^{2}+\frac{1}{2\eta}\|x-y\|^{2}-\frac{1}{2\eta}\|x^{*}-y\|^{2}$
$\displaystyle=$ $\displaystyle
f(x)-f(x^{*})+\frac{\mu}{2}\|x-x^{*}\|^{2}+\mu\langle
x-x^{*},x^{*}-x^{0}\rangle+\frac{1}{2\eta}\|x-x^{*}\|^{2}+\frac{1}{\eta}\langle
x-x^{*},x^{*}-y\rangle$ $\displaystyle=$ $\displaystyle
f(x)-f(x^{*})+\langle\mu(x^{*}-x^{0})+\frac{1}{\eta}(x^{*}-y),x-x^{*}\rangle+\frac{1}{2\eta_{\mu}}\|x-x^{*}\|^{2}$
$\displaystyle=$ $\displaystyle f(x)-f(x^{*})-\langle\nabla
f(x^{*}),x-x^{*}\rangle+\frac{1}{2\eta_{\mu}}\|x-x^{*}\|^{2}.$ (32)

The above inequality and (16) with $(u,v)=(x,x^{*})$ imply that

$g^{\eta}(x)-g^{\eta}(x^{*})\leq\frac{L}{2}\|x-x^{*}\|^{2}+\frac{1}{2\eta_{\mu}}\|x-x^{*}\|^{2}.$

Using the above inequality and the definitions of $\eta_{\mu,L}$ and $h_{2}$
in (18) and (17), we conclude that the second inequality $g^{\eta}(x)\leq
h_{2}(x)$ holds.

Proof of Proposition 5.2: Using the definition of $h_{2}$ in (17) and Lemma
A.1 with ${\lambda}=\eta_{\mu,L}$, we have

$\displaystyle\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))dx$
$\displaystyle=\exp(-g^{\eta}(x^{*}))\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta_{\mu,L}}\|x-x^{*}\|^{2}\right)dx$
$\displaystyle=\exp(-g^{\eta}(x^{*}))(2\pi\eta_{\mu,L})^{d/2}.$

The above identity, the fact that $g^{\eta}(\tilde{x}_{j})\geq
g^{\eta}(x^{*})$, and the definitions of $\eta_{\mu}$ and $\eta_{\mu,L}$ imply
that

$\displaystyle\mathbb{P}\left(U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))}\right)$
$\displaystyle\geq\frac{\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))dx}{\exp\left(-g^{\eta}(\tilde{x}_{j})+\delta\right)(2\pi\eta_{\mu})^{d/2}}$
$\displaystyle\geq\exp(g^{\eta}(\tilde{x}_{j})-g^{\eta}(x^{*})-\delta)\left(\frac{\eta_{\mu,L}}{\eta_{\mu}}\right)^{d/2}\geq\exp(-\delta)\left(\frac{1}{1+\eta_{\mu}L}\right)^{d/2}.$

The above inequality and the assumption that $\eta_{\mu}\leq 1/(Ld)$ imply
that the expected number of the iterations is

$\frac{1}{\mathbb{P}\left(U\leq\frac{\exp(-g^{\eta}(X))}{\exp(-h_{1}(X))}\right)}\leq\exp(\delta)(1+\eta_{\mu}L)^{d/2}\leq\exp(\delta)\left(1+\frac{1}{d}\right)^{d/2}\leq\exp(1/2+\delta),$

where the last inequality is due to the fact that $(1+1/d)^{d}\leq e$.

Proof of Theorem 5.4: It follows from Theorem 2.2 that the iteration-
complexity (i.e., number of calls to RGO) for Algorithm 1 to obtain
$\varepsilon$ total tolerance to $\pi$ is (13). Using (13) with $\eta\leq
1/(Ld)$, we have the number of calls to RGO is (19). Moreover, it follows from
Proposition 5.3 with $\eta\leq 1/(Ld)$ that Algorithm 3 has iteration-
complexity $\tilde{\cal O}(1)$. As a consequence, the total iteration-
complexity bound for obtaining $\varepsilon$ total tolerance to $\pi$ in terms
of total variation is (19).


# Unsupervised Motor Imagery Saliency Detection Based on Self-Attention
Mechanism

Navid Ayoobi∗ and Elnaz Banan Sadeghian∗ ∗Navid Ayoobi and Elnaz Banan
Sadeghian are with the Department of Electrical and Computer Engineering,
Stevens Institute of Technology, Hoboken, USA. (E-mail:
nayoobi<EMAIL_ADDRESS>
###### Abstract

Detecting the salient parts of motor-imagery electroencephalogram (MI-EEG)
signals can enhance the performance of the brain-computer interface (BCI)
system and reduce the computational burden required for processing lengthy MI-
EEG signals. In this paper, we propose an unsupervised method based on the
self-attention mechanism to detect the salient intervals of MI-EEG signals
automatically. Our suggested method can be used as a preprocessing step within
any BCI algorithm to enhance its performance. The effectiveness of the
suggested method is evaluated on the most widely used BCI algorithm, the
common spatial pattern (CSP) algorithm, using dataset 2a from BCI competition
IV. The results indicate that the proposed method can effectively prune MI-EEG
signals and significantly enhance the performance of the CSP algorithm in
terms of classification accuracy.

##  I INTRODUCTION

Brain-computer interface (BCI) bypasses the need for peripheral nervous and
neuromuscular systems for communicating with surroundings and enables
individuals to control external devices by directly translating
neurophysiological brain signals into commands [1, 2]. Among all signal
acquisition methods in BCI, electroencephalography is employed in the vast
majority of non-invasive systems. Further, motor imagery electroencephalogram
(MI-EEG) patterns captured when an individual imagines performing a movement
are often used as the most discriminatory patterns for classification.
However, these MI-EEG signals have a poor signal-to-noise ratio (SNR).
Moreover, numerous forms of artifacts like electromyogram (EMG), power line
interference, electrooculogram (EOG), eye movement, and blinks decline the SNR
even further [3]. Due to the low SNR, unprocessed MI-EEG signals result in
poor detection and classification performance. Therefore, enhancing SNR is
necessary for developing a reliable BCI system.

Artifact removal is one of the most popular preprocessing methods to increase
the SNR of EEG signals. In [4], a method is proposed based on blind source
separation to remove EMG artifacts from EEG data with a limited number of
electrodes. Chavez et al. in [5] adopt the wavelet decomposition and a
resampling method to detect and remove ocular and muscular artifacts from a
single-channel EEG across different scales. However, these methods attempt to
remove the artifacts from the entire EEG signal where the user performance can
vary while performing an MI task varies during a trial. Therefore, using the
entire trial is not optimally suitable for classification and degrades the
classification accuracy. Moreover, processing a lengthy EEG is computationally
inefficient. Numerous efforts have been conducted to address this issue by
using a time window to detect the most salient interval in a trial for MI
classification [1, 6]: Gaur et al. in [1] use a $2$-second sliding window and
split the EEG signal into nine intervals where the CSP algorithm is used to
extract features and nine linear discriminant analysis (LDA) classifiers are
trained separately for each individual interval. The longest consecutive
repetition is then used to choose the best label among the nine labels
generated by different LDAs. Hsu et al. in [6] propose a method based on
continuous wavelet transform and Student’s two-sample t-statistics to detect
one salient interval within every $4$-second trial using a $1$-second
searching window. However, these approaches do not detect the entirety of
salient intervals which can be several and of varying lengths. Furthermore,
they require labeled data, and their effectiveness depends on the selected
classifier’s ability to separate them.

This paper proposes a method based on the self-attention mechanism [7] to
automatically extract several salient intervals in every MI-EEG trial.
Further, the method is unsupervised and therefore it does not require any
labeled data. The proposed method can be applied to MI-EEG signals as a
preprocessing step within any BCI algorithm to improve the SNR and enhance its
performance. The proposed architecture is comprised of an input embedding
layer, an encoder-decoder network, a self-attention layer, and a dense layer.
We utilize the encoder-decoder structure to reconstruct the input, yielding an
unsupervised architecture. The self-attention layer assigns a weight to each
input time sample and aids the decoder to precisely reconstruct the input
based on those time samples with the highest weights representing the salient
parts of each trial. Consequently, the proposed method can detect different
numbers of salient intervals with varying lengths. After the training phase,
the decoder part is removed from the network. The trained encoder and self-
attention layer is then used to extract salient intervals of unseen testing
trials.

To evaluate our approach, we use dataset 2a from BCI competition IV to extract
the salient intervals of the MI trials. We compare the performance of the most
widely used algorithm in BCI, the common spatial pattern (CSP) [8] algorithm,
with and without our proposed method. The results show that the proposed
method can effectively extract the salient intervals of MI-EEG signals and
significantly improve the CSP performance in terms of classification accuracy.

The remainder of this article is organized as follows. Section II introduces
the dataset. Section III elaborates the proposed method and architecture.
Section IV presents the simulation results and discussion. Section V concludes
this article.

##  II Dataset

To evaluate our proposed method, we use dataset 2a from BCI competition IV [9]
that contains the EEG data of nine healthy subjects. This dataset consists of
movement imaginations of left and right hands, feet, and tongue. The EEG
signals are recorded using twenty-two Ag/AgCl electrodes and sampled at $250$
Hz. The timing scheme of each trial is depicted in Fig. 1. In this paper, we
only utilize the left-hand and right-hand trials. We also extract all our
trials from the second $2$ after the cue to the second $6$.

Figure 1: Timing scheme of each trial in dataset 2a from BCI competition IV.

##  III Methods

In this section, we elaborate on our proposed unsupervised method for
detecting the salient parts of EEG trials to improve classification accuracy.
The suggested architecture consists of four main parts: input embedding, long
short-term memory (LSTM) encoder-decoder, self-attention, and dense layer. The
encoder-decoder structure helps the self-attention layer learn the salient
parts of the EEG trials in an unsupervised method by reconstructing the input.
The proposed architecture is shown in Fig. 2.

###  III-A The proposed architecture

Input embedding: Given a multivariate signal as input, the embedding layer
captures the dependencies between different channels [10] in EEG signals and
creates a representation vector for each time sample. We employ a 1D
convolutional layer with $m$ numbers of $d$-dimensional kernels. This layer
produces an $m$-dimensional vector for each time sample by looking at current
and $d{-}1$ next time samples. We utilize zero padding to ensure that the
input and output sequences have the same length.

LSTM encoder: LSTM as a variant of recurrent neural networks can capture long
temporal dependencies in sequential data [11]. We employ LSTM units to build
our encoder and decoder. The encoder processes the input sequence of length
$T$ and summarizes the temporal information of each time sample in two
$h$-dimensional vectors, the hidden state $h_{e}$ and the cell state $c_{e}$.
The produced hidden states are then fed to a self-attention layer. After $T$
recurrences, the whole sequence is summarized in final encoder hidden state
$h_{eT}$ and cell state $c_{eT}$. These two vectors are fed to the decoder’s
initial hidden state $h_{d0}$ and cell state $c_{d0}$.

Figure 2: The block diagram of the proposed method.

Self-attention: Attention is originally proposed in [7] to improve the
performance of an encoder-decoder in a machine translation. It alleviates the
bottleneck problem of using the fixed-length vector from which the decoder
produces the output sentence by allowing the model to automatically seek
source sentence parts that are significant for predicting a target word. Self-
attention, also known as intra-attention, is an attention mechanism that links
distinct points in a single sequence to calculate a representation of the
sequence [12]. In this work, we adopt the self-attention mechanism to detect
the salient parts of each MI-EEG trial. In our self-attention layer, the
$i$-th encoder hidden state $h_{ei}{=}\\{h_{ei,1},...\,,h_{ei,h}\\}$ is mapped
to an $n_{k}$-dimensional query vector $q_{i}$, an $n_{k}$-dimensional key
vector $k_{i}$, and an $n_{v}$-dimensional value vector $v_{i}$ by multiplying
$h_{ei}$ with three matrices $W_{q}$, $W_{k}$ and $W_{v}$, respectively:

$q_{i}=W_{q}h_{ei}$ (1) $k_{i}=W_{k}h_{ei}$ (2) $v_{i}=W_{v}h_{ei}$ (3)

These vectors are used to calculate a context vector in each position in the
sequential data. Here, dot-product attention [12] is used to compute the
context vectors. The context vector at $i$-th position is a weighted sum of
value vectors in all positions

$c_{i}=\sum_{t=1}^{T}\alpha_{i,t}v_{t},$ (4)

where the attention weight $\alpha_{i,t}$ measures how well the $i$-th encoder
hidden state aligns with the encoder hidden state at position $t$. The
attention weights are computed as the inner product of the query with the
corresponding key. After computing the weights associated with all keys, a
softmax function is used to make the weights sum up to one. If all query
vectors, key vectors, and value vectors are packed together into $\mathbf{Q}$,
$\mathbf{K}$, and $\mathbf{V}$ matrices, the context matrix $\mathbf{C}$ is
calculated as follows

$\mathbf{C}=softmax(\mathbf{Q}\mathbf{K}^{T})\mathbf{V}.$ (5)

LSTM decoder and dense layer: The LSTM decoder accepts the final hidden state
and cell state of the encoder as its initial hidden state and cell state. The
context vectors generated by the self-attention layer are fed to the decoder
inputs. The decoder outputs are then passed through a dense layer to increase
the dimensionality of decoder outputs to the number of EEG channels. The
decoder and dense layer learn to reconstruct the original MI trial by
compelling the self-attention layer to pay more attention to certain parts of
the MI trial that are more significant for reconstructing the input.

###  III-B Training and Testing stages

Training stage: To improve the performance of the proposed network in
reconstructing the input EEG, $x$, we randomly select $p_{1}$ time samples of
$x$ and set $p_{2}$ feature values to zero to obtain $\tilde{x}$. In the
training stage, $\tilde{x}$ is used as the training trial, and its
corresponding original EEG signal $x$ is used as the ground truth. The mean
squared error is used as the loss function

$Q=\frac{1}{Tn_{c}}\sum_{t=1}^{T}\sum_{i=1}^{n_{c}}(\hat{x}_{t}^{(i)}-x_{t}^{(i)})^{2},$
(6)

where $n_{c}$, $\hat{x}$ and $x$ are the number of EEG channels, the
reconstructed EEG, and the original EEG, respectively.

Testing stage: The decoder and dense layer are removed from the network for
the testing stage. A test trial is passed through the embedding layer, the
encoder, and the self-attention layer. The computed attention weights
$\\{\alpha\\}$ are used to detect the salient parts of the test trial. Assume
that the attention weights of a test trial are represented in a matrix
$\Lambda_{T\times T}$. The attention at position $t$, $a_{t}$, is computed as
follows

$a_{t}=\frac{1}{T}\sum_{j=1}^{T}\Lambda_{j,t}.$ (7)

The computed attentions form an attention vector $A=[a_{1},\,...\,,a_{T}]$. We
segment the attention vector into $n$ intervals for separating salient
intervals. $r$ out of $n$ intervals with the highest average attention are
selected. These two hyperparameters are tuned in the cross validation
procedure. The time samples corresponding to the selected attention intervals
are then concatenated to form the pruned trial.

##  IV Results and discussion

Figure 3: Comparison of the performance of the CSP algorithm for different
subjects in the dataset before and after applying our proposed method in terms
of classification accuracy.

In all our experiments, the MI-EEG signals are bandpass filtered between $4$
Hz and $40$ Hz with a fifth-order Butterworth filter. For all subjects, the
number of kernels $m$, kernel dimension $d$, LSTM hidden state dimension $h$,
key-query dimension $n_{k}$, the dimension of the value vector $n_{v}$,
parameter $p_{1}$, and parameter $p_{2}$ are set to $5$, $4$, $4$, $4$, $4$,
$0.6$, and $0.4$, respectively. A deep neural network with three hidden layers
(with 5,10, and 15 nodes) is used for the dense layer. The parameters $r$ and
$n$ are tuned by 5-fold cross validation for each subject. In addition, we use
three pairs of spatial filters in the CSP algorithm [8].

###  IV-A The effect of the proposed method on the classification performance

To evaluate the proposed method, we compare the performance of the CSP
algorithm under two scenarios: 1) with and 2) without applying our method to
the MI-EEG trials.

In the first scenario, we train our proposed architecture with the $4$-second
training trials. After the training stage, we use the encoder and self-
attention layer to extract the salient intervals from the training and testing
trials. The spatial filters of the CSP algorithm are then designed using the
pruned training trials. The accuracy of the CSP algorithm is computed using
the pruned testing trials. In the second scenario, the $4$-second training
trials and testing trials are directly used to design the spatial filters and
obtain the accuracy of the CSP algorithm, respectively. In both scenarios,
linear discriminant analysis is used as the classifier.

The results are shown in Fig. 3 for all the subjects. We observe that the
classification accuracy is improved approximately by $14.2\%$, $19.4\%$,
$6.6\%$, $31.9\%$, $24.7\%$, $9.2\%$, $10.4\%$, $1.5\%$, and $7.4\%$ for
subject one through nine, respectively. On average, the proposed method
enhances the classification accuracy of the CSP by about $13.9\%$ across all
subjects.

In addition, the ratio of the length of the pruned trials to the length of the
unpruned trials is $0.7$, $0.65$, $0.45$, $0.7$, $0.55$, $0.7$, $0.75$,
$0.75$, and $0.75$ for subject 1 through 9, respectively. Hence, the proposed
method reduces the computational burden required for processing EEG signals
and increases the classification accuracy simultaneously.

###  IV-B The effect of segment length on attention vector splitting

This section evaluates the effect of different segment lengths $T/n$ used for
splitting the attention vector to extract the salient intervals. Fig. 4
compares the results in terms of classification accuracy for subject one in
the dataset versus the number of time samples in each interval, $T/n$. $r$
intervals with the highest average attention are concatenated to form a trial
containing $\ell{=}rT/n{=}750$ (red dashed line), $\ell{=}550$ (green dashed
line), $\ell{=}350$ (orange dashed line), and $\ell{=}250$ (blue dashed line)
time samples. The black horizontal line shows the classification accuracy of
the CSP algorithm without our method. We observe that restricting the
concatenated signal to have a very short length ($\ell{=}350$ and
$\ell{=}250$) results in lower accuracy compared to the CSP method alone. As
shown the best curves correspond to $\ell{=}750$ and $\ell{=}550$. These
curves prove that using the entire lengths of the trials does not always
achieve the best obtainable accuracy.

In addition, the results suggest that shorter segment lengths $T/n$ can
generally provide higher accuracy than longer segment lengths. In other words,
using multiple time windows instead of a single time window [1, 6] results in
detecting more salient intervals in a trial. The proposed method may be
employed as a prepossessing step in MI-BCI studies for enhancing motor imagery
classification.

Figure 4: The effect of different segment lengths on attention vector
splitting for subject one.

##  V Conclusion

In this paper, we present an unsupervised method based on the self-attention
mechanism to automatically detect the salient intervals of MI-EEG trials. We
test our proposed method using the CSP algorithm on the dataset 2a from BCI
competition IV. The results show that pruning the MI-EEG trials by the
proposed method can effectively improve the SNR and reduce the computational
burden in BCI systems. The results also show that the average classification
accuracy of the CSP algorithm is improved by approximately $13.9\%$ across all
subjects. The suggested approach may be used as a preprocessing step before
any BCI algorithm to improve its performance.

## References

  * [1] P. Gaur, H. Gupta, A. Chowdhury, K. McCreadie, R. B. Pachori, and H. Wang, “A sliding window common spatial pattern for enhancing motor imagery classification in EEG-BCI,” _IEEE Transactions on Instrumentation and Measurement_ , vol. 70, pp. 1–9, 2021.
  * [2] E. Banan Sadeghian and M. H. Moradi, “Continuous detection of motor imagery in a four-class asynchronous BCI,” in _2007 29th Annual International Conference of the IEEE Engineering in Medicine and Biology Society_ , 2007, pp. 3241–3244.
  * [3] Q. Wei, S. Zhu, Y. Wang, X. Gao, H. Guo, and X. Wu, “Maximum signal fraction analysis for enhancing signal-to-noise ratio of EEG signals in SSVEP-based BCIs,” _IEEE Access_ , vol. 7, pp. 85 452–85 461, 2019\.
  * [4] L. Zou, X. Chen, G. Dang, Y. Guo, and Z. J. Wang, “Removing muscle artifacts from EEG data via underdetermined joint blind source separation: A simulation study,” _IEEE Transactions on Circuits and Systems II: Express Briefs_ , vol. 67, no. 1, pp. 187–191, 2019.
  * [5] M. Chavez, F. Grosselin, A. Bussalb, F. D. V. Fallani, and X. Navarro-Sune, “Surrogate-based artifact removal from single-channel EEG,” _IEEE Transactions on Neural Systems and Rrehabilitation Engineering_ , vol. 26, no. 3, pp. 540–550, 2018.
  * [6] W.-Y. Hsu, C.-C. Lin, M.-S. Ju, and Y.-N. Sun, “Wavelet-based fractal features with active segment selection: Application to single-trial EEG data,” _Journal of Neuroscience Methods_ , vol. 163, no. 1, pp. 145–160, 2007.
  * [7] D. Bahdanau, K. Cho, and Y. Bengio, “Neural machine translation by jointly learning to align and translate,” _arXiv preprint arXiv:1409.0473_ , 2014\.
  * [8] B. Blankertz, R. Tomioka, S. Lemm, M. Kawanabe, and K.-R. Muller, “Optimizing spatial filters for robust EEG single-trial analysis,” _IEEE Signal Processing Magazine_ , vol. 25, no. 1, pp. 41–56, 2007.
  * [9] C. Brunner, R. Leeb, G. Müller-Putz, A. Schlögl, and G. Pfurtscheller, “BCI competition 2008–graz data set a,” _Institute for Knowledge Discovery (Laboratory of Brain-Computer Interfaces), Graz University of Technology_ , vol. 16, pp. 1–6, 2008.
  * [10] H. Song, D. Rajan, J. J. Thiagarajan, and A. Spanias, “Attend and diagnose: Clinical time series analysis using attention models,” in _Thirty-Second AAAI Conference on Artificial Intelligence_ , 2018.
  * [11] K. Greff, R. K. Srivastava, J. Koutník, B. R. Steunebrink, and J. Schmidhuber, “LSTM: A search space odyssey,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 28, no. 10, pp. 2222–2232, 2016\.
  * [12] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser, and I. Polosukhin, “Attention is all you need,” in _Advances in Neural Information Processing Systems_ , 2017, pp. 5998–6008.