# Proximal Oracles for Optimization and Sampling

Jiaming Liang Yongxin Chen Goergen Institute for Data Science & Department of
Computer Science, University of Rochester, Rochester, NY 14620. (email:
jiaming.liang@rochester.edu).School of Aerospace Engineering, Georgia
Institute of Technology, Atlanta, GA, 30332. (email: yongchen@gatech.edu).

(April 2, 2024)

###### Abstract

We consider convex optimization with non-smooth objective function and log-
concave sampling with non-smooth potential (negative log density). In
particular, we study two specific settings where the convex
objective/potential function is either semi-smooth or in composite form as the
finite sum of semi-smooth components. To overcome the challenges caused by
non-smoothness, our algorithms employ two powerful proximal frameworks in
optimization and sampling: the proximal point framework for optimization and
the alternating sampling framework (ASF) that uses Gibbs sampling on an
augmented distribution. A key component of both optimization and sampling
algorithms is the efficient implementation of the proximal map by the
regularized cutting-plane method. We establish the iteration-complexity of the
proximal map in both semi-smooth and composite settings. We further propose an
adaptive proximal bundle method for non-smooth optimization. The proposed
method is universal since it does not need any problem parameters as input.
Additionally, we develop a proximal sampling oracle that resembles the
proximal map in optimization and establish its complexity using a novel
technique (a modified Gaussian integral). Finally, we combine this proximal
sampling oracle and ASF to obtain a Markov chain Monte Carlo method with non-
asymptotic complexity bounds for sampling in semi-smooth and composite
settings.

Key words. Non-smooth optimization, proximal point method, universal method,
high-dimensional sampling, Markov chain Monte Carlo, complexity analysis

##  1 Introduction

We are interested in convex optimization problems

$\min_{x\in\mathbb{R}^{d}}f(x)$ (1)

as well as log-concave sampling problems

${\rm sample}\quad\nu(x)\propto\exp(-f(x)),$ (2)

where $f:\mathbb{R}^{d}\to\mathbb{R}$ is convex but not necessarily smooth.

Optimization and sampling are two of the most important algorithmic tools at
the interface of data science and computation. Optimization has been
extensively studied in machine learning, communication, and supply chain.
Drawing samples from a given (often unnormalized) probability density plays a
crucial role in many scientific and engineering problems that face uncertainty
(either physically or algorithmically). Sampling algorithms are widely used in
many areas such as statistical inference/estimation, operations research,
physics, biology, and machine learning, etc [1, 8, 9, 12, 17, 18, 22, 46]. For
instance, in Bayesian inference, one draws samples from the posterior
distribution to infer its mean, covariance, or other important statistics.
Sampling is also heavily used in molecular dynamics to discover new molecular
structures.

This work is along the recent line of researches that lie in the interface of
sampling and optimization [7, 44]. Indeed, sampling is closely related to
optimization. On the one hand, optimization can be viewed as the limit of
sampling when the temperature parameter, or equivalently the randomness in the
problem, goes to $0$. On the other hand, sampling can be viewed as an
optimization over the manifold of probability distributions [48, 49]. The
popular gradient-based Markov chain Monte Carlo (MCMC) methods such as
Langevin Monte Carlo (LMC) [6, 15, 39, 41], Metropolis-adjusted Langevin
algorithm (MALA) [2, 40, 41], and Hamiltonian Monte Carlo (HMC) [34] resemble
the gradient-based algorithms in optimization and can be viewed as the
sampling counterparts of them.

The goal of this paper is to develop efficient proximal algorithms to solve
optimization problems (1) as well as to draw samples from potentials (2),
where both $f$ in (1) and (2) lack smoothness. In particular, we consider two
settings where the convex objective/potential function $f$ is either semi-
smooth (i.e., the (sub)gradient $f^{\prime}$ is Hölder-continuous with
exponent $\alpha\in[0,1]$) or a composite function with multiple semi-smooth
components. The core of both proximal optimization and sampling algorithms is
a proximal map of $f$. Based on an efficient implementation of the proximal
map, we develop an adaptive proximal bundle method to solve (1). Building upon
the proximal map of $f$ and rejection sampling, we propose a highly efficient
method to realize a proximal sampling oracle, which is used in a proximal
sampling framework [23, 5] in the same spirit as the proximal point method for
optimization. With those proximal oracles for optimization and sampling in
hand, we are finally able to establish the complexity to sample from densities
with non-smooth potentials.

We summarize our contributions as follows.

  * i)

We analyze the complexity bounds for implementing the proximal map of $f$
using the regularized cutting-plane method in both semi-smooth and composite
settings (Section 3).

  * ii)

We develop an adaptive proximal bundle method (APBM) using the regularized
cutting-plane method and a novel adaptive stepsize strategy in the proximal
point method, and establish the complexity bound for semi-smooth optimization
(Section 4). APBM is a universal method as it does not need any problem-
dependent parameters as input.

  * iii)

We propose an efficient scheme to realize the proximal sampling oracle that
lacks smoothness and establish novel techniques to bound its complexity.
Combining the proximal sampling oracle and the proximal sampling framework, we
obtain a general proximal sampling algorithm for convex semi-smooth and
composite potentials. Finally, we establish complexity bounds for the proximal
sampling algorithm in both cases (Section 5).

##  2 Proximal optimization and sampling

The proximal point framework (PPF), proposed in [30] and further developed in
[42, 43] (see [38] for a modern and comprehensive monograph), is a general
class of optimization algorithms that involve solving a sequence of
subproblems of the form

$x_{k+1}\leftarrow\mathrm{argmin}\,\left\\{f(x)+\frac{1}{2\eta}\|x-x_{k}\|^{2}:x\in\mathbb{R}^{d}\right\\},$
(3)

where $\eta>0$ is a prox stepsize and $\leftarrow$ means the subproblem can be
solved either exactly or approximately. When the exact solution is available,
we denote

$x_{k+1}={\rm prox}_{\eta f}(x_{k}),$

where ${\rm prox}_{f}(\cdot)$ is called a proximal map of $f$ and defined as

${\rm
prox}_{f}(y):=\mathrm{argmin}\,\left\\{f(x)+\frac{1}{2}\|x-y\|^{2}:x\in\mathbb{R}^{d}\right\\}.$
(4)

If the subproblem (3) does not admit a closed-form solution, it can usually be
solved with standard or specialized iterative methods.

Many classical first-order methods in optimization, such as the proximal
gradient method, the proximal subgradient method, the primal-dual hybrid
gradient method of [4] (also known as the Chambolle-Pock method), the
extragradient method of [21] are instances of PPF. It is worth noting that, by
showing that the alternating direction method of multipliers (ADMM) as an
instance of PPF, [32] gives the first iteration-complexity result of ADMM for
solving a class of linearly constrained convex programming problems. Inspired
by the PPF viewpoint, [27] develops a variant of the proximal bundle method
and establishes the optimal iteration-complexity, which is the first optimal
complexity result for proximal bundle methods. Recent works [19, 26, 20, 29]
have also applied PPF to solve weakly convex optimization and weakly convex-
concave min-max problems.

Proximal map in sampling. Sampling shares many similarities with optimization.
An interesting connection between the two problems is through the algorithm
design and analysis from the perspective of PPF. The alternating sampling
framework (ASF) introduced in [23] is a generic framework for sampling from a
distribution $\pi^{X}(x)\propto\exp(-f(x))$. Analogous to PPF in optimization,
ASF with stepsize $\eta>0$ repeats the two steps as in Algorithm 1.

Algorithm 1 Alternating Sampling Framework [23]

1\. Sample $y_{k}\sim\pi^{Y|X}(y\mid
x_{k})\propto\exp\left(-\frac{1}{2\eta}\|x_{k}-y\|^{2}\right)$

2\. Sample $x_{k+1}\sim\pi^{X|Y}(x\mid
y_{k})\propto\exp\left(-f(x)-\frac{1}{2\eta}\|x-y_{k}\|^{2}\right)$

ASF is a special case of Gibbs sampling [13] of the joint distribution

$\pi(x,y)\propto\exp\left(-f(x)-\frac{1}{2\eta}\|x-y\|^{2}\right).$

Starting from the original paper [23] that proposes ASF, a subsequent paper
[5] provides an improved analysis of ASF, and [50] also studies Gibbs sampling
based on ASF over bipartite graphs. In Algorithm 1, sampling $y_{k}$ given
$x_{k}$ in step 1 can be easily done since $\pi^{Y|X}(y\mid x_{k})={\cal
N}(x_{k},\eta I)$ is a simple Gaussian distribution. Sampling $x_{k+1}$ given
$y_{k}$ in step 2 is however a nontrivial task; it corresponds to the so-
called restricted Gaussian oracle (RGO) for $f$ introduced in [23], which is
defined as follows.

######  Definition 2.1.

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

RGO is an analog of the proximal map (4) in optimization. To use ASF in
practice, one needs to efficiently implement RGO. Some examples of $f$ that
admit a computationally efficient RGO have been presented in [33, 45]. These
instances of $f$ have simple structures such as coordinate-separable
regularizers, $\ell_{1}$-norm, and group Lasso. To apply ASF on a general
potential function $f$, developing an efficient implementation of the RGO is
essential.

A rejection sampling-based implementation of RGO for general convex nonsmooth
potential function $f$ with bounded Lipschitz constant is given in [24]. If
the stepsize $\eta$ is small enough, then it only takes a constant number of
rejection steps to generate a sample according to RGO in expectation. Another
exact realization of RGO is provided in [25] for nonconvex composite potential
$f$ satisfying Hölder continuous conditions. It is also shown that the
expected number of rejections to implement RGO is a small constant if $\eta$
is small enough. Other inexact realizations of RGO based on approximate
rejection sampling are studied in [14, 10]. In all these implementations, a
key step is realizing the proximal map (4). It is worth noting that [25] also
connects ASF with other well-known Langevin-type sampling algorithms such as
Langevin Monte Carlo (LMC) and Proximal Langevin Monte Carlo (PLMC) via RGO.
In a nutshell, [25] shows that both LMC and PLC are instances of ASF but with
approximate implementations of RGO, which always accept the sample from the
proposal distribution without rejection. Hence, this explains why the samples
generated by LMC is biased while ASF is unbiased.

Based on the cutting-plane method, this paper develops a generic and efficient
implementation of the proximal map (4) and applies the proximal map in both
optimization and sampling. For optimization, we use this proximal map and an
adaptive stepsize rule to design a universal bundle method. For sampling, we
combine this proximal map and rejection sampling to realize the RGO, and then
propose a practical and efficient proximal sampling algorithm based on it.

For both optimization and sampling, we consider two specific scenarios: 1) $f$
is semi-smooth, i.e., $f$ satisfies

$\|f^{\prime}(u)-f^{\prime}(v)\|\leq L_{\alpha}\|u-v\|^{\alpha},\quad\forall
u,v\in\mathbb{R}^{d},$ (5)

where $f^{\prime}$ denotes a subgradient of $f$, $\alpha\in[0,1]$, and
$L_{\alpha}>0$; and 2) $f$ is a composite function of semi-smooth components,
i.e., $f$ satisfies

$\|f^{\prime}(u)-f^{\prime}(v)\|\leq\sum_{i=1}^{n}L_{\alpha_{i}}\|u-v\|^{\alpha_{i}},\quad\forall
u,v\in\mathbb{R}^{d},$ (6)

where $\alpha_{i}\in[0,1]$ and $L_{\alpha_{i}}>0$ for every $1\leq i\leq n$.
When $\alpha=0$, (5) reduces to a Lipschitz continuous condition, and when
$\alpha=1$, it reduces to a smoothness condition. It follows from (5) and (6)
that for every $u,v\in\mathbb{R}^{d}$,

$f(u)-f(v)-\langle
f^{\prime}(v),u-v\rangle\leq\frac{L_{\alpha}}{\alpha+1}\|u-v\|^{\alpha+1},$
(7)

and

$f(u)-f(v)-\langle
f^{\prime}(v),u-v\rangle\leq\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|u-v\|^{\alpha_{i}+1}.$
(8)

The proof is given in Appendix A.

##  3 Algorithm and complexities for the proximal subproblem

The proximal subproblem (3) generally does not admit a closed-form solution.
We design an iterative method that approximately solves (3) and derive the
corresponding iteration-complexities for semi-smooth and composite $f$ in
Subsections 3.1 and 3.2, respectively.

Given a point $y\in\mathbb{R}^{d}$, we consider the optimization problem

$f_{y}^{\eta}(x^{*})=\min\left\\{f_{y}^{\eta}(x)=f(x)+\frac{1}{2\eta}\|x-y\|^{2}:x\in\mathbb{R}^{d}\right\\}$
(9)

and aim at obtaining a $\delta$-solution, i.e., a point $\bar{x}$ such that
$f_{y}^{\eta}(\bar{x})-f_{y}^{\eta}(x^{*})\leq\delta$. In both semi-smooth and
composite settings, we use a regularized cutting-plane method (Algorithm 2),
which is usually used in the proximal bundle method [27, 28] for solving
convex non-smooth optimization problems. We remark that though Algorithm 2 is
widely used in the proximal bundle method and is not new, the complexity
analyses (i.e., Theorems 3.5 and 3.9) for semi-smooth and composite functions
$f$ are lacking.

Algorithm 2 Regularized Cutting-plane Method

1\. Let $y\in\mathbb{R}^{d}$, $\eta>0$, and $\delta>0$ be given, and set
$x_{0}=\tilde{x}_{0}=y$, and $j=1$.

2\. Update $f_{j}(x)=\max\left\\{f(x_{i})+\langle
f^{\prime}(x_{i}),x-x_{i}\rangle:\,0\leq i\leq j-1\right\\}$.

3\. Compute
$x_{j}=\underset{x\in\mathbb{R}^{d}}{\mathrm{argmin}\,}\left\\{f_{j}^{\eta}(x):=f_{j}(x)+\frac{1}{2\eta}\|x-y\|^{2}\right\\},$
(10)
$\tilde{x}_{j}=\mathrm{argmin}\,\left\\{f_{y}^{\eta}(x):x\in\\{x_{j},\tilde{x}_{j-1}\\}\right\\}.$
(11)

4\. If $f_{y}^{\eta}(\tilde{x}_{j})-f_{j}^{\eta}(x_{j})\leq\delta$, then stop
and return $J=j,x_{J},\tilde{x}_{J}$; else, go to step 5.

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

The basic idea of Algorithm 2 is to approximate $f$ with piece-wise affine
functions constructed by a collection of cutting-planes and solve the
resulting simplified problem (10). As the approximation becomes more and more
accurate, the best approximate solution $\tilde{x}_{j}$ converges to the
solution $x^{*}$ to (9). Subproblem (10) can be reformulated into convex
quadratic programming with $j$ affine constraints and hence is solvable.

The following technical lemma summarizes basic properties of Algorithm 2. It
is useful in the complexity analysis for both optimization and sampling.

######  Lemma 3.1.

Assume $f$ is convex. For every $j\geq 1$, define

$\delta_{j}:=f_{y}^{\eta}(\tilde{x}_{j})-f_{j}^{\eta}(x_{j}).$ (12)

Let $J,x_{J},\tilde{x}_{J}$ be the outputs of Algorithm 2, then the following
statements hold:

  * a)

$f_{j}(x)\leq f_{j+1}(x)$ and $f_{j}(x)\leq f(x)$ for $\forall
x\in\mathbb{R}^{d}$ and $\forall j\geq 1$;

  * b)

$f_{j}^{\eta}(x_{j})+\|x-x_{j}\|^{2}/(2\eta)\leq f_{j}^{\eta}(x)$ for $\forall
x\in\mathbb{R}^{d}$ and $\forall j\geq 1$;

  * c)

$\delta_{J}\leq\delta$ and
$\delta_{j+1}+\frac{1}{2\eta}\|x_{j+1}-x_{j}\|^{2}\leq\delta_{j}$ for $\forall
j\geq 1$;

  * d)

$f_{y}^{\eta}(\tilde{x}_{J})-f_{y}^{\eta}(x)\leq\delta-\frac{1}{2\eta}\|x_{J}-x\|^{2}$
for $\forall x\in\mathbb{R}^{d}$;

  * e)

$-\frac{1}{\eta}(x^{*}-y)\in\partial f(x^{*})$ where $\partial f$ denotes the
subdifferential of $f$.

Proof: a) The first inequality follows from the definition of $f_{j}$ in step
2 of Algorithm 2. The second inequality directly follows from the definition
of $f_{j}$ and the convexity of $f$.

b) This statement follows from (10) and the fact that $f_{j}^{\eta}$ is
$(1/\eta)$-strongly convex.

c) This first inequality immediately follows from (12) and step 4 of Algorithm
2.

Using the first inequality in 3.1(a) and 3.1(b) with $x=x_{j+1}$, we obtain

$f_{j+1}^{\eta}(x_{j+1})\geq f_{j}^{\eta}(x_{j+1})\geq
f_{j}^{\eta}(x_{j})+\frac{1}{2\eta}\|x_{j+1}-x_{j}\|^{2}.$

This inequality, the definition of $\tilde{x}_{j}$ in (11), and the definition
of $\delta_{j}$ in (12) imply that

$\displaystyle\delta_{j+1}$
$\displaystyle=f_{y}^{\eta}(\tilde{x}_{j+1})-f_{j+1}^{\eta}(x_{j+1})\leq
f_{y}^{\eta}(\tilde{x}_{j})-f_{j}^{\eta}(x_{j})-\frac{1}{2\eta}\|x_{j+1}-x_{j}\|^{2}=\delta_{j}-\frac{1}{2\eta}\|x_{j+1}-x_{j}\|^{2}.$

d) Using the second inequality in (a), (b) with $j=J$, and the first
inequality in (c), we have

$\displaystyle f(\tilde{x}_{J})$
$\displaystyle-f(x)+\frac{1}{2\eta}\|x-x_{J}\|^{2}\stackrel{{\scriptstyle\text{(a)}}}{{\leq}}f(\tilde{x}_{J})-f_{J}(x)+\frac{1}{2\eta}\|x-x_{J}\|^{2}$
$\displaystyle\stackrel{{\scriptstyle\text{(b)}}}{{\leq}}f(\tilde{x}_{J})-f_{J}^{\eta}(x_{J})+\frac{1}{2\eta}\|x-y\|^{2}\stackrel{{\scriptstyle\text{(c)}}}{{\leq}}\delta-\frac{1}{2\eta}\|\tilde{x}_{J}-y\|^{2}+\frac{1}{2\eta}\|x-y\|^{2}.$

This statement then follows from rearranging the terms and the definition of
$f_{y}^{\eta}$ in (9)

e) This statement directly follows from the optimality condition of (9).

Clearly, when Algorithm 2 terminates, the output $\tilde{x}_{J}$ is a
$\delta$-solution to (9). To see this, note that, using the first inequality
in Lemma 3.1(c), (10), and the fact that $f_{J}^{\eta}(\cdot)\leq
f_{y}^{\eta}(\cdot)$, we have

$f_{y}^{\eta}(\tilde{x}_{J})\leq\delta+f_{J}^{\eta}(x_{J})\overset{\eqref{def:xj}}{\leq}\delta+f_{J}^{\eta}(x^{*})\leq\delta+f_{y}^{\eta}(x^{*}).$

It is also easy to see that $\delta_{j}$ is computable upper bound on the gap
$f_{y}^{\eta}(\tilde{x}_{j})-f_{y}^{\eta}(x^{*})$. Hence, Algorithm 2
terminates when $\delta_{j}\leq\delta$.

###  3.1 Complexity for semi-smooth optimization

This subsection is devoted to the complexity analysis of Algorithm 2 for
solving (9) where $f$ is semi-smooth, i.e., satisfying (5). The following
lemma provides basic recursive formulas and is the starting point of the
analysis of Algorithm 2.

######  Lemma 3.2.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth. Then, for every $j\geq 1$,
the following statements hold:

  * a)

$\delta_{j}\leq\frac{L_{\alpha}}{\alpha+1}\|x_{j}-x_{j-1}\|^{\alpha+1}$;

  * b)

$\delta_{j+1}+\frac{1}{2\eta}\left(\frac{\alpha+1}{L_{\alpha}}\delta_{j+1}\right)^{\frac{2}{\alpha+1}}\leq\delta_{j}$.

Proof: a) It follows from the definition of $\delta_{j}$ in (12) and the
definition of $\tilde{x}_{j}$ in (11) that

$\displaystyle\delta_{j}$ $\displaystyle\leq
f_{y}^{\eta}(x_{j})-f_{j}^{\eta}(x_{j})=f(x_{j})-f_{j}(x_{j})$
$\displaystyle\leq f(x_{j})-f(x_{j-1})-\langle
f^{\prime}(x_{j-1}),x_{j}-x_{j-1}\rangle$
$\displaystyle\leq\frac{L_{\alpha}}{\alpha+1}\|x_{j}-x_{j-1}\|^{\alpha+1},$

where the second inequality is due to the definition of $f_{j}$ in the step 2
of Algorithm 2, and the third inequality is due to (7) with
$(u,v)=(x_{j},x_{j-1})$.

b) This statement directly follows from a) and the second inequality in Lemma
3.1(c).

We know from Lemma 3.1(c) that $\\{\delta_{j}\\}_{j\geq 1}$ is non-increasing.
The next proposition gives a bound on $j$ so that $\delta_{j}\leq\delta$,
i.e., the termination criterion in step 4 of Algorithm 2 is satisfied.

######  Proposition 3.3.

Define

$\beta:=\frac{1}{2\eta}\left(\frac{\alpha+1}{L_{\alpha}}\right)^{\frac{2}{\alpha+1}}\delta^{\frac{1-\alpha}{\alpha+1}},\quad
j_{0}=1+\left\lceil\frac{1+\beta}{\beta}\log\left(\frac{\delta_{1}}{\delta}\right)\right\rceil.$
(13)

Then, the following statements hold:

  * a)

if $\delta_{j}>\delta$, then $(1+\beta)\delta_{j}\leq\delta_{j-1}$;

  * b)

$\delta_{j}\leq\delta$ for every $j\geq j_{0}$.

As a consequence, the iteration count $J$ in Algorithm 2 satisfies $J\leq
j_{0}$.

Proof: a) Using the definition of $\beta$ in (13), the assumption that
$\delta_{j}>\delta$, and Lemma 3.2(b), we obtain

$(1+\beta)\delta_{j}=\delta_{j}+\frac{1}{2\eta}\left(\frac{\alpha+1}{L_{\alpha}}\right)^{\frac{2}{\alpha+1}}\delta^{\frac{1-\alpha}{\alpha+1}}\delta_{j}\leq\delta_{j}+\frac{1}{2\eta}\left(\frac{\alpha+1}{L_{\alpha}}\delta_{j}\right)^{\frac{2}{\alpha+1}}\leq\delta_{j-1}.$

b) Since $\\{\delta_{j}\\}_{j\geq 1}$ is non-increasing, it suffices to prove
that $\delta_{j_{0}}\leq\delta$. We prove this statement by contradiction.
Suppose that $\delta_{j_{0}}>\delta$, then we have $\delta_{j}>\delta$ for
$j\leq j_{0}$. Hence, statement (a) holds for $j\leq j_{0}$. Using this
conclusion repeatedly and the fact that $\tau\leq\exp(\tau-1)$ with
$\tau=1/(1+\beta)$, we have

$\delta_{j_{0}}\leq\frac{1}{(1+\beta)^{j_{0}-1}}\delta_{1}\leq\exp\left(-\frac{\beta}{1+\beta}(j_{0}-1)\right)\delta_{1}\leq\delta,$

where the last inequality is due to the definition of $j_{0}$ in (13). This
contradicts with the assumption that $\delta_{j_{0}}>\delta$, and hence we
prove this statement.

The following result shows that $\delta_{1}$ is bounded from above, and hence
the bound in Proposition 3.3 is meaningful.

######  Lemma 3.4.

For a given $y\in\mathbb{R}^{d}$, we have

$\delta_{1}\leq\frac{L_{\alpha}\eta^{\alpha+1}}{\alpha+1}\|f^{\prime}(y)\|^{\alpha+1}.$

Proof: Following the optimality condition of (10) with $j=1$, we have
$x_{0}-x_{1}=\eta f^{\prime}(x_{0})=\eta f^{\prime}(y)$. This identity and
Lemma 3.2(a) with $j=1$ then imply that the lemma holds.

We now conclude the iteration-complexity bound for Algorithm 2.

######  Theorem 3.5.

Algorithm 2 takes $\tilde{\cal O}\left(\eta
L_{\alpha}^{\frac{2}{\alpha+1}}\left(\frac{1}{\delta}\right)^{\frac{1-\alpha}{\alpha+1}}+1\right)$
iterations to terminate.

Proof: This theorem follows directly from Proposition 3.3 and Lemma 3.4.

###  3.2 Complexity for composite optimization

This subsection is devoted to the complexity analysis of Algorithm 2 for
solving (9) where $f$ is a composite function satisfying (6). The following
lemma is an analogue of Lemma 3.2 and provides key recursive formulas for
$\delta_{j}$, which is defined in (12).

######  Lemma 3.6.

Assume $f$ is convex and satisfies (6). For $\delta>0$, define

$M=\sum_{i=1}^{n}\frac{L_{\alpha_{i}}^{\frac{2}{\alpha_{i}+1}}}{[(\alpha_{i}+1)\delta]^{\frac{1-\alpha_{i}}{\alpha_{i}+1}}}.$
(14)

Then, for every $j\geq 1$, the following statements hold:

  * a)

$\delta_{j}\leq\frac{M}{2}\|x_{j}-x_{j-1}\|^{2}+\sum_{i=1}^{n}(1-\alpha_{i})\frac{\delta}{2}$;

  * b)

$\left(1+\frac{1}{\eta
M}\right)\left(\delta_{j+1}-\sum_{i=1}^{n}(1-\alpha_{i})\frac{\delta}{2}\right)\leq\delta_{j}-\sum_{i=1}^{n}(1-\alpha_{i})\frac{\delta}{2}$.

Proof: a) Following a similar argument as in the proof of Lemma 3.2(a) with
(7) replaced by (8), we have

$\delta_{j}\leq\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|u-v\|^{\alpha_{i}+1}.$
(15)

Using the Young’s inequality $ab\leq a^{p}/p+b^{q}/q$ with

$a=\frac{L_{\alpha}}{(\alpha+1)\delta^{\frac{1-\alpha}{2}}}\|x_{j}-x_{j-1}\|^{\alpha+1},\quad
b=\delta^{\frac{1-\alpha}{2}},\quad p=\frac{2}{\alpha+1},\quad
q=\frac{2}{1-\alpha},$

we obtain

$\frac{L_{\alpha}}{\alpha+1}\|x_{j}-x_{j-1}\|^{\alpha+1}\leq\frac{L_{\alpha}^{\frac{2}{\alpha+1}}}{2[(\alpha+1)\delta]^{\frac{1-\alpha}{\alpha+1}}}\|x_{j}-x_{j-1}\|^{2}+\frac{(1-\alpha)\delta}{2}.$

Combining the above inequality and (15), and using the definition of $M$ in
(14), we prove the statement.

b) It immediately follows from (a) and the second inequality in Lemma 3.1(c)
that

$\delta_{j+1}+\frac{1}{\eta
M}\left(\delta_{j+1}-\sum_{i=1}^{n}(1-\alpha_{i})\frac{\delta}{2}\right)\leq\delta_{j+1}+\frac{1}{2\eta}\|x_{j+1}-x_{j}\|^{2}\leq\delta_{j},$

and hence the statement follows.

The following lemma gives an upper bound on $\delta_{1}$ similar to Lemma 3.4.

######  Lemma 3.7.

For a given $y\in\mathbb{R}^{d}$, we have

$\delta_{1}\leq\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|f^{\prime}(y)\|^{\alpha_{i}+1}.$

Proof: This lemma follows from a similar argument as in the proof of Lemma
3.4.

The following proposition is the key result in establishing the iteration-
complexity of Algorithm 2.

######  Proposition 3.8.

We have $\delta_{j}\leq\delta$, for every $j$ such that

$j\geq(1+\eta M)\log\left(\frac{2\delta_{1}}{\delta}\right).$ (16)

Proof: Let

$\tau=\frac{\eta M}{1+\eta M},$ (17)

then Lemma 3.6(b) becomes

$\delta_{j+1}-\sum_{i=1}^{n}(1-\alpha_{i})\frac{\delta}{2}\leq\tau\left(t_{j}-\sum_{i=1}^{n}(1-\alpha_{i})\frac{\delta}{2}\right).$

Using the above inequality repeatedly and the fact that
$\tau\leq\exp(\tau-1)$, we have for every $j\geq 1$,

$\delta_{j}-\frac{(1-\alpha)\delta}{2}\leq\tau^{j-1}\left(\delta_{1}-\frac{(1-\alpha)\delta}{2}\right)\leq\tau^{j-1}\delta_{1}\leq\exp\\{(\tau-1)(j-1)\\}\delta_{1}.$

Hence, it is easy to see that $\delta_{j}\leq\delta$ if
$j\geq\frac{1}{1-\tau}\log\left(\frac{2\delta_{1}}{\delta}\right)$. Using the
definition of $\tau$ in (17), we have if $j$ is as in (16), then
$\delta_{j}\leq\delta$.

We are ready to present the complexity bound for Algorithm 2.

######  Theorem 3.9.

Algorithm 2 takes $\tilde{\cal O}\left(\eta M+1\right)$ iterations to
terminate, where $M$ is as in (14).

Proof: This theorem follows directly from Proposition 3.8 and Lemma 3.7.

##  4 Adaptive proximal bundle method

As discussed in Section 3, the cutting-plane method (i.e., Algorithm 2) is
widely used in the proximal bundle method as a subroutine to repeatedly solve
the proximal subproblem (3). Since the proximal bundle method uses a more
accurate cutting-plane model $f_{j}$ rather than a linearization as an
approximation of the objective function $f$, it generalizes the subgradient
method and is able to work with weaker regularization, namely larger stepsize
$\eta$. This explains why both methods have optimal complexity bounds [27,
28], but the proximal bundle method is always more efficient in practice.

For both the subgradient method and the proximal bundle method to have the
optimal performance, one needs to carefully select the stepsize $\eta$,
namely, being small enough for the subgradient method and within a certain
(but relatively large) range for the proximal bundle method. In both cases, we
need to know the problem-dependent parameters such as $\alpha$ and
$L_{\alpha}$, which are unknown or hard to estimate in practice. In this
section, we develop the APBM based on an adaptive stepsize strategy for (3)
and using Algorithm 2 to solve each subproblem (3). We also discuss variants
of adaptive subgradient methods and compare them with APBM. For simplicity, we
only present the analysis of semi-smooth functions satisfying (5), while the
composite functions satisfying (6) can be similarly analyzed using results
from Subsection 3.2.

From practical observations, the proximal bundle method works well when the
number of inner iterations (i.e., those of Algorithm 2 to solve (3)) stays as
a constant much larger than 1 (i.e., that of the subgradient method), say
$10$. Recall from Theorem 3.5 that inner complexity is $\tilde{\cal
O}\left(\eta
L_{\alpha}^{\frac{2}{\alpha+1}}\left(\frac{1}{\delta}\right)^{\frac{1-\alpha}{\alpha+1}}+1\right)$.
Since we do not know $\alpha$ and $L_{\alpha}$, we cannot choose a constant
stepsize $\eta$ so that the number of inner iterations is close to a desired
number such as $10$. Hence, an adaptive stepsize rule is indeed needed.

By carefully examining Proposition 3.3 and Theorem 3.5, we find that the inner
complexity is $\tilde{\cal O}(\beta^{-1}+1)$ where $\beta$ is as in (13).
Suppose we want to prescribe the number of inner iterations to be close to
$\beta_{0}^{-1}$ for some $\beta_{0}\in(0,1]$, if $\beta_{0}\leq\beta$, then
by Proposition 3.3(a), we have

$(1+\beta_{0})\delta_{j}\leq\delta_{j-1}.$ (18)

Hence, it suffices to begin with a relatively large $\eta$, check (18) to
determine whether the $\eta$ is small enough (i.e., $\beta$ is large enough),
and adjust $\eta$ (if necessary) by progressively halving it.

Algorithm 3 below is a formal statement of APBM based on the above intuition.

Algorithm 3 Adaptive Proximal Bundle Method

1\. Let $y_{0}\in\mathbb{R}^{d}$, $\eta_{0}>0$, $\beta_{0}\in(0,1]$, and
$\varepsilon>0$ be given, and set $k=1$.

2\. Call Algorithm 2 with $(y,\eta,\delta)=(y_{k-1},\eta_{k-1},\varepsilon/2)$
and output $(y_{k},\tilde{y}_{k})=(x_{J},\tilde{x}_{J})$.

3\. In the execution of Algorithm 2, if (18) is always true, then set
$\eta_{k}=\eta_{k-1}$; otherwise, set $\eta_{k}=\eta_{k-1}/2$.

4\. Set $k$ $\leftarrow$ $k+1$ and go to step 2.

The following lemma provides basic results of Algorithm 2 and is the starting
point of the analysis of APBM.

######  Lemma 4.1.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth. The following statements
hold for APBM:

  * a)

for every $k\geq 1$ and $u\in\mathbb{R}^{d}$, we have

$2\eta_{k-1}[f(\tilde{y}_{k})-f(u)]\leq\|y_{k-1}-u\|^{2}-\|y_{k}-u\|^{2}+\eta_{k-1}\varepsilon;$
(19)

  * b)

for any $k\geq 1$, if

$\eta_{k-1}\leq\frac{1}{2\beta_{0}}\left(\frac{\alpha+1}{L_{\alpha}}\right)^{\frac{2}{\alpha+1}}\left(\frac{\varepsilon}{2}\right)^{\frac{1-\alpha}{\alpha+1}},$
(20)

then $\eta_{k}=\eta_{k-1}$;

  * c)

$\\{\eta_{k}\\}$ is a non-increasing sequence;

  * d)

for every $k\geq 0$,

$\eta_{k}\geq\underline{\eta}:=\min\left\\{\frac{1}{4\beta_{0}}\left(\frac{\alpha+1}{L_{\alpha}}\right)^{\frac{2}{\alpha+1}}\left(\frac{\varepsilon}{2}\right)^{\frac{1-\alpha}{\alpha+1}},\eta_{0}\right\\}.$
(21)

Proof: a) It follows from Lemma 3.1(d) that for every $u\in\mathbb{R}^{d}$

$2\eta[f(\tilde{x}_{J})-f(x)]\leq 2\eta\delta+\|u-y\|^{2}-\|x_{J}-u\|^{2}.$

Noting from step 2 of Algorithm 3 that
$(\delta,\eta,y,x_{J},\tilde{x}_{J})=(\varepsilon/2,\eta_{k-1},y_{k-1},y_{k},\tilde{y}_{k})$,
which together with the above inequality, implies that (19) holds.

b) It follows from Proposition 3.3(a) and (20) that (18) always holds in the
execution of Algorithm 2. In view of step 3 of Algorithm 3, there holds
$\eta_{k}=\eta_{k-1}$.

c) This statement clearly follows from step 3 of Algorithm 3.

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

The following theorem gives the total iteration-complexity of APBM.

######  Theorem 4.2.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth. If
$\eta_{0}\leq\|y_{0}-x_{*}\|^{2}/\varepsilon$, then the iteration-complexity
to obtain an $\varepsilon$-solution to (1) is given by

$\tilde{\cal
O}\left(\frac{L_{\alpha}^{\frac{2}{\alpha+1}}\|y_{0}-x_{*}\|^{2}}{\varepsilon^{\frac{2}{\alpha+1}}}+\eta_{0}L_{\alpha}^{\frac{2}{\alpha+1}}\left(\frac{1}{\varepsilon}\right)^{\frac{1-\alpha}{\alpha+1}}\log\left(\frac{\eta_{0}}{\underline{\eta}}\right)+1\right)$
(22)

where $\underline{\eta}$ is as in (21).

Proof: Noting from Lemma 4.1 (d) that $\eta_{k}$ is bounded from below, we
know there exists some $\tilde{\eta}\in[\underline{\eta},\eta_{0}]$ such that
for some $k_{0}\geq 1$, $\eta_{k}\equiv\tilde{\eta}$ for $k\geq k_{0}$. Thus,
it follows from the assumption that
$\eta_{0}\leq\|y_{0}-x_{*}\|^{2}/\varepsilon$ that

$\underline{\eta}\leq\tilde{\eta}\leq\frac{\|y_{0}-x_{*}\|^{2}}{\varepsilon}.$
(23)

We consider the worst-case scenario where APBM keeps halving the stepsize
until it is stable at $\tilde{\eta}$ and the convergence relies on the
conservative stepsize $\tilde{\eta}$. Summing (19) from $k=1$ to $n$, we have

$2\sum_{k=1}^{n}\eta_{k-1}\left(\min_{1\leq k\leq
n}f(\tilde{y}_{k})-f(u)\right)\leq
2\sum_{k=1}^{n}\eta_{k-1}[f(\tilde{y}_{k})-f(u)]\leq\|y_{0}-u\|^{2}-\|y_{n}-u\|^{2}+\varepsilon\sum_{k=1}^{n}\eta_{k-1}.$

The above inequality with $u=x_{*}$, the fact that $\eta_{k}\leq\eta_{0}$, and
the assumption that $\eta_{k}\equiv\tilde{\eta}$ for $k\geq k_{0}$ imply that

$\|y_{n}-x_{*}\|^{2}\leq\|y_{0}-x_{*}\|^{2}+n\eta_{0}\varepsilon$ (24)

and

$\min_{1\leq k\leq
n}f(\tilde{y}_{k})-f_{*}\leq\frac{\|y_{0}-x_{*}\|^{2}}{2\sum_{k=1}^{n}\eta_{k-1}}+\frac{\varepsilon}{2}\leq\frac{\|y_{0}-x_{*}\|^{2}}{2(n-k_{0})\tilde{\eta}}+\frac{\varepsilon}{2}.$

In order to have $\min_{1\leq k\leq n}f(\tilde{y}_{k})-f_{*}\leq\varepsilon$,
we need

$n-k_{0}={\cal
O}\left(\frac{\|y_{0}-x_{*}\|^{2}}{\tilde{\eta}\varepsilon}+1\right).$ (25)

Moreover, it follows from the way $\eta_{k}$ is updated in step 3 and Lemma
4.1(d) that

$k_{0}={\cal
O}\left(\log\left(\frac{\eta_{0}}{\tilde{\eta}}\right)+1\right)={\cal
O}\left(\log\left(\frac{\eta_{0}}{\underline{\eta}}\right)+1\right).$ (26)

Indeed, (24) holds with $n$ replaced by any $k\leq n$ and

$\|y_{k}-x_{*}\|^{2}\leq\|y_{0}-x_{*}\|^{2}+n\eta_{0}\varepsilon.$

It thus follows from (25) and (26) that $\\{y_{k}\\}$ is bounded. As a result,
using Lemma 3.4, we can derive a uniform bound on $\delta_{1}$ for every call
to Algorithm 2. Now, using Theorem 3.5, we have the iteration-complexity of
every call to Algorithm 2 is uniformly bounded by

$\tilde{\cal
O}\left(\tilde{\eta}L_{\alpha}^{\frac{2}{\alpha+1}}\left(\frac{1}{\varepsilon}\right)^{\frac{1-\alpha}{\alpha+1}}+1\right)$
(27)

for every cycle $k\geq k_{0}$ and by

$\tilde{\cal
O}\left(\eta_{0}L_{\alpha}^{\frac{2}{\alpha+1}}\left(\frac{1}{\varepsilon}\right)^{\frac{1-\alpha}{\alpha+1}}+1\right)$
(28)

for every cycle $k\leq k_{0}-1$. Hence, multiplying (25) and (27) and using
(23) and the definition of $\underline{\eta}$ in (21), we obtain the
iteration-complexity

$\tilde{\cal
O}\left(\frac{L_{\alpha}^{\frac{2}{\alpha+1}}\|y_{0}-x_{*}\|^{2}}{\varepsilon^{\frac{2}{\alpha+1}}}+1\right)$

for cycles $k\geq k_{0}$, and multiplying (26) and (28), we obtain the
iteration-complexity

$\tilde{\cal
O}\left(\eta_{0}L_{\alpha}^{\frac{2}{\alpha+1}}\left(\frac{1}{\varepsilon}\right)^{\frac{1-\alpha}{\alpha+1}}\log\left(\frac{\eta_{0}}{\underline{\eta}}\right)+1\right)$

for cycles $k\leq k_{0}-1$. Finally, the total iteration-complexity (22)
clearly follows from the above two bounds.

Discussion on other universal methods Several universal methods based on the
backtracking line-search procedure have been studied in the literature. Paper
[37] considers the same semi-smooth problem (with an additional composite
function $h$) as in this paper. To finds an $\varepsilon$-solution of (1), the
universal primal gradient method proposed in [37] starts from an initial pair
$(\hat{x}_{0},\eta_{0})$ and in the $(j+1)$-th iteration searches for a pair
$(x_{\eta},\eta)$ satisfying a condition

$f(x_{\eta})-\ell_{f}(x_{\eta};\hat{x}_{j})-\frac{1}{2\eta}\|x_{\eta}-\hat{x}_{j}\|^{2}\leq\frac{\varepsilon}{2},$
(29)

where $\ell_{f}(u;v)=f(v)+\langle f^{\prime}(v,u-v\rangle$ and

$x_{\eta}=\underset{u\in\mathbb{R}^{n}}{\mathrm{argmin}\,}\left\\{\ell_{f}(u;\hat{x}_{j})+h(u)+\frac{1}{2\eta}\|u-\hat{x}_{j}\|^{2}\right\\}.$
(30)

If the condition (29) is not satisfied, then the method rejects the pair, sets
$\eta\leftarrow\eta/2$, and updates $x_{\eta}$ as in (30) with the new $\eta$,
otherwise, it accepts the pair and sets
$(\hat{x}_{j+1},\eta_{j+1})=(x_{\eta},\eta)$. Paper [28] proposes an adaptive
composite subgradient (A-CS) method for solving (1) where $f$ satisfies

$\|f^{\prime}(x)-f^{\prime}(y)\|\leq 2M_{f}+L_{f}\|x-y\|,\quad\forall
x,y\in\mathbb{R}^{d}.$ (31)

It is shown in Proposition 2.1 of [28] that any function $f$ that satisfies

$\|f^{\prime}(x)-f^{\prime}(y)\|\leq
2M_{\alpha}+L_{\alpha}\|x-y\|^{\alpha},\quad\forall x,y\in\mathbb{R}^{d},$

for some $\alpha\in[0,1]$ also satisfies (31) with

$M_{f}(\theta):=M_{\alpha}+\frac{L_{\alpha}\theta}{2},\quad
L_{f}(\theta):=L_{\alpha}\alpha\left(\frac{1-\alpha}{\theta}\right)^{\frac{1-\alpha}{\alpha}}$

for any $\theta>0$. Hence, the semi-smooth functions (5) considered in this
paper are included in the class of functions satisfying (31). More
interestingly, a careful look at A-CS of [28] and the universal primal
gradient method of [37] reveals that the two methods are identical.

The universal primal gradient method is essentially an adaptive subgradient
method and the convergence of subgradient methods relies on small enough
stepsizes, so it is natural to enforce (29) to make the method adaptive.
However, the bundle method converges with any constant stepsize $\eta$ since
it guarantees the condition $\delta_{j}\leq\varepsilon/2$, which is in the
same spirit of (29), by the cutting-plane approach (i.e., Algorithm 2) but not
by small $\eta$. Therefore, it is not necessary to use a small $\eta$ in every
iteration of each call to Algorithm 2. Instead of frequently reducing $\eta$,
by the introduction of $\beta_{0}$, APBM develops a way to regulate the
complexity of Algorithm 2 and adjust $\eta$ only when (18) is not always true
in the previous call to Algorithm 2. Another difference between APBM and the
universal primal gradient method is that the latter rejects all the pairs
$(x_{\eta},\eta)$ until (29) is satisfied, but APBM always accepts the output
of Algorithm 2 even if (18) is not true for every iteration in Algorithm 2.
Therefore, APBM potentially employs a larger stepsize $\eta$ than the
universal primal gradient method and is thus a more relaxed adaptive method.

Discussion on optimal universal methods The lower complexity bound for solving
(1) is shown in [35] to be

${\cal
O}\left(\left(\frac{L_{\alpha}\|y_{0}-x_{*}\|^{1+\alpha}}{\varepsilon}\right)^{\frac{2}{1+3\alpha}}\right).$

The well-known Nesterov’s accelerated gradient method has been shown in [36]
to match the above complexity bound and hence is an optimal method. On the
other hand, the dominant term of bound (22) is its first term and it is only
optimal when $\alpha=0$, i.e., $f$ is $L_{0}$-Lipschitz continuous. Motivated
by [36], it is possible to develop optimal universal methods based on the
accelerated gradient method. This requires accelerated schemes in both PPF and
Algorithm 2. Paper [31] proposes an accelerated variant of PPF, which is
extended by [3, 16, 11] to obtain optimal $p$-th order methods with
convergence rate ${\cal O}(k^{-(3p+1)/2})$ for $p\geq 2$.

##  5 Proximal sampling algorithm

Assuming the RGO in the ASF can be realized, the ASF exhibits remarkable
convergence properties. It was shown in [23] that Algorithm 1 converges
linearly when $f$ is strongly convex. This convergence result is recently
improved in [5] under various weaker assumptions on the target distribution
$\pi^{X}\propto\exp(-f)$. Below we present several convergence results
established in [5] that will be used in this paper, under the assumptions that
$\pi^{X}$ is log-concave, or satisfies the log-Sobolev inequality or Poincaré
inequality (PI). Recall that a probability distribution $\nu$ satisfies PI
with constant $C_{\rm PI}>0$ ($1/C_{\rm PI}$-PI) if for any smooth bounded
function $\psi:\mathbb{R}^{d}\to\mathbb{R}$,

$\mathbb{E}_{\nu}[(\psi-\mathbb{E}_{\nu}(\psi))^{2}]\leq C_{\rm
PI}\mathbb{E}_{\nu}[\left\|\nabla\psi\right\|^{2}].$

To this end, for two probability distributions $\rho\ll\nu$, we denote by

$H_{\nu}(\rho):=\int\rho\log\frac{\rho}{\nu},\quad\chi_{\nu}^{2}(\rho):=\int\frac{\rho^{2}}{\nu}-1$

the KL divergence and the Chi-squared divergence, respectively. We denote by
$W_{2}$ the Wasserstein-2 distance

$W_{2}^{2}(\nu,\rho):=\min_{\gamma\in\Pi(\nu,\rho)}\int\|x-y\|^{2}{\mathrm{d}}\gamma(x,y),$

where $\Pi(\nu,\rho)$ represents the set of all couplings between $\nu$ and
$\rho$.

######  Theorem 5.1 ([5, Theorems 2 & 4]).

We denote by $\rho_{k}^{X}$ the law of $x_{k}$ of Algorithm 1 starting from
any initial distribution $\rho_{0}^{X}$. Then, the following statements hold:

  * a)

if $\pi^{X}\propto\exp(-f)$ is log-concave (i.e., $f$ is convex), then
$H_{\pi^{X}}(\rho_{k}^{X})\leq W_{2}^{2}(\rho_{0}^{X},\pi^{X})/(k\eta)$;

  * b)

if $\pi^{X}\propto\exp(-f)$ satisfies ${\lambda}$-PI, then
$\chi_{\pi^{X}}^{2}(\rho^{X}_{k})\leq\chi_{\pi^{X}}^{2}(\rho^{X}_{0})/{(1+{\lambda}\eta)}^{2k}$.

As discussed earlier, to use ASF in sampling problems, we need to realize the
RGO with efficient implementations. In the rest of this section, we develop
efficient algorithms for RGO associated with the two scenarios of sampling we
are interested in, and then combine them with the ASF to establish a proximal
algorithm for sampling. The complexity of the proximal algorithm can be
obtained by combining the above convergence results for ASF and the complexity
results we develop for RGO. The rest of the section is organized as follows.
In Subsection 5.1 we develop an efficient algorithm for RGO associated with
semi-smooth potentials via rejection sampling. This is combined with ASF to
obtain an efficient sampling algorithm from semi-smooth potentials. In
Subsection 5.2, we further extend results to the second setting, i.e.,
composite potentials.

###  5.1 Sampling from semi-smooth potentials

The bottleneck of using the ASF (Algorithm 1) in sampling tasks with general
distributions is the availability of RGO implementations. In this subsection,
we address this issue for convex semi-smooth potentials by developing an
efficient algorithm for the corresponding RGO.

Our algorithm of RGO for $f$ is based on rejection sampling. We use a special
proposal, namely a Gaussian distribution centered at the $\delta$-solution of
(9), which is obtained by invoking Algorithm 2. With this proposal and a
sufficiently small $\eta>0$, the expected number of rejection sampling steps
to obtain one effective sample turns out to be bounded from above by a
dimension-free constant. To bound the complexity of the rejection sampling, we
develop a novel technique to estimate a modified Gaussian integral (see
Proposition 5.3).

To this end, let $J,\tilde{x}_{J},x_{J}$ be the outputs of Algorithm 2 and
define

$\displaystyle h_{1}$ $\displaystyle:=\frac{1}{2\eta}\|\cdot-
x_{J}\|^{2}+f_{y}^{\eta}(\tilde{x}_{J})-\delta,$ (32a) $\displaystyle h_{2}$
$\displaystyle:=\frac{1}{2\eta}\|\cdot-x^{*}\|^{2}+\frac{L_{\alpha}}{\alpha+1}\|\cdot-x^{*}\|^{\alpha+1}+f_{y}^{\eta}(x^{*}).$
(32b)

Note that $h_{2}$ is only used for analysis and thus the fact it depends on
$x^{*}$ is not an issue. Algorithm 4 describes the implementation of RGO for
$f$ based on Algorithm 2 and rejection sampling.

Algorithm 4 RGO Implementation based on Rejection Sampling

1\. Run Algorithm 2 to compute $x_{J}$ and $\tilde{x}_{J}$.

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

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

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

######  Lemma 5.2.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth. Let $f_{y}^{\eta}$ be as in
(9) and $h_{1}$ and $h_{2}$ be as in (32). Then, for every
$x\in\mathbb{R}^{d}$, we have

$h_{1}(x)\leq f_{y}^{\eta}(x)\leq h_{2}(x).$ (33)

Proof: The first inequality in (33) immediately follows from Lemma 3.1(d) and
the definition of $h_{1}$ in (32a). By the definition of $f_{y}^{\eta}$ in (9)
we get

$\displaystyle f_{y}^{\eta}(x)-f_{y}^{\eta}(x^{*})=$ $\displaystyle
f(x)-f(x^{*})+\frac{1}{2\eta}\|x-y\|^{2}-\frac{1}{2\eta}\|x^{*}-y\|^{2}$
$\displaystyle=$ $\displaystyle
f(x)-f(x^{*})+\frac{1}{2\eta}\|x-x^{*}\|^{2}+\frac{1}{\eta}\langle
x-x^{*},x^{*}-y\rangle.$ (34)

It follows from Lemma 3.1(e) and (7) with $(u,v)=(x,x^{*})$ that

$f(x)-f(x^{*})+\frac{1}{\eta}\langle
x^{*}-y,x-x^{*}\rangle\leq\frac{L_{\alpha}}{\alpha+1}\|x-x^{*}\|^{\alpha+1},$

which together with (34) implies that

$f_{y}^{\eta}(x)-f_{y}^{\eta}(x^{*})\leq\frac{L_{\alpha}}{\alpha+1}\|x-x^{*}\|^{\alpha+1}+\frac{1}{2\eta}\|x-x^{*}\|^{2}.$

Using the above inequality and the definition of $h_{2}$ in (32b), we conclude
that the second inequality in (33) holds.

From the expression of $h_{1}$ in (32a), it is clear that the proposal
distribution $\exp(-h_{1}(x))$ is a Gaussian centered at $x_{J}$. To achieve a
tight bound on the expected runs of the rejection sampling, we use a function
$h_{2}$ which is not quadratic; the standard choice of quadratic function does
not give as tight results due to the lack of smoothness. To use this $h_{2}$
in the complexity analysis, we need to estimate the integral
$\int\exp(-h_{2})$, which turns out to be a highly nontrivial task. Below we
establish a technical result on a modified Gaussian integral, which will be
used later to bound the integral $\int\exp(-h_{2})$ and hence the complexity
of the RGO rejection sampling in Algorithm 4.

######  Proposition 5.3.

Let $\alpha\in[0,1]$, $\eta>0$, $a\geq 0$ and $d\geq 1$. If

$2a(\eta d)^{(\alpha+1)/2}\leq 1,$ (35)

then

$\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}-a\|x\|^{\alpha+1}\right)\mathrm{d}x\geq\frac{(2\pi\eta)^{d/2}}{2}.$
(36)

Proof: Denote $r=\|x\|$, then

$\mathrm{d}x=r^{d-1}\mathrm{d}r\mathrm{d}S^{d-1},$

where $\mathrm{d}S^{d-1}$ is the surface area of the $(d-1)$-dimensional unit
sphere. It follows that

$\displaystyle\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}-a\|x\|^{\alpha+1}\right)\mathrm{d}x$
$\displaystyle=\int_{0}^{\infty}\int\exp\left(-\frac{1}{2\eta}r^{2}-ar^{\alpha+1}\right)r^{d-1}\mathrm{d}r\mathrm{d}S^{d-1}$
$\displaystyle=\frac{2\pi^{d/2}}{\Gamma\left(\frac{d}{2}\right)}\int_{0}^{\infty}\exp\left(-\frac{1}{2\eta}r^{2}-ar^{\alpha+1}\right)r^{d-1}\mathrm{d}r.$
(37)

In the above equation, 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)$
where $\Gamma(\cdot)$ is the gamma function, i.e.,

$\Gamma(z)=\int_{0}^{\infty}t^{z-1}e^{-t}\mathrm{d}t.$ (38)

Defining

$F_{d,\eta}(a):=\int_{0}^{\infty}\exp\left(-\frac{1}{2\eta}r^{2}-ar^{\alpha+1}\right)r^{d}\mathrm{d}r,$
(39)

to establish (36), it suffices to bound $F_{d-1,\eta}(a)$ from below.

It follows directly from the definition of $F_{d,\eta}$ in (39) that

$\frac{\mathrm{d}F_{d-1,\eta}(a)}{\mathrm{d}a}=\int_{0}^{\infty}\exp\left(-\frac{1}{2\eta}r^{2}-ar^{\alpha+1}\right)(-r^{\alpha+1})r^{d-1}\mathrm{d}r=-F_{d+\alpha,\eta}(a).$

This implies $F_{d,\eta}$ is monotonically decreasing and thus
$F_{d+\alpha,\eta}(a)\leq F_{d+\alpha,\eta}(0)$. As a result,

$\frac{\mathrm{d}F_{d-1,\eta}(a)}{\mathrm{d}a}\geq-F_{d+\alpha,\eta}(0)$

and therefore,

$F_{d-1,\eta}(a)\geq F_{d-1,\eta}(0)-aF_{d+\alpha,\eta}(0).$ (40)

Setting $t=r^{2}/(2\eta)$, we can write

$\displaystyle F_{d,\eta}(0)$
$\displaystyle=\int_{0}^{\infty}\exp\left(-\frac{1}{2\eta}r^{2}\right)r^{d}\mathrm{d}r=\int_{0}^{\infty}e^{-t}(2\eta
t)^{\frac{d-1}{2}}\eta\mathrm{d}t$
$\displaystyle=2^{\frac{d-1}{2}}\eta^{\frac{d+1}{2}}\int_{0}^{\infty}e^{-t}t^{\frac{d-1}{2}}\mathrm{d}t.$
(41)

In view of the definition of the gamma function (38), we obtain

$F_{d,\eta}(0)=2^{\frac{d-1}{2}}\eta^{\frac{d+1}{2}}\Gamma\left(\frac{d+1}{2}\right).$
(42)

Applying the Wendel’s double inequality (53) yields

$\frac{\Gamma\left(\frac{d+\alpha+1}{2}\right)}{\Gamma\left(\frac{d}{2}\right)}\leq\left(\frac{d}{2}\right)^{\frac{\alpha+1}{2}}.$

Using (40), (42), the above inequality and the assumption (35), we have

$\displaystyle F_{d-1,\eta}(a)$ $\displaystyle\geq
F_{d-1,\eta}(0)-aF_{d+\alpha,\eta}(0)$
$\displaystyle=2^{\frac{d}{2}-1}\eta^{\frac{d}{2}}\Gamma\left(\frac{d}{2}\right)-a2^{\frac{d+\alpha-1}{2}}\eta^{\frac{d+\alpha+1}{2}}\Gamma\left(\frac{d+\alpha+1}{2}\right)$
$\displaystyle=2^{\frac{d}{2}-1}\eta^{\frac{d}{2}}\Gamma\left(\frac{d}{2}\right)\left(1-a2^{\frac{\alpha+1}{2}}\eta^{\frac{\alpha+1}{2}}\frac{\Gamma\left(\frac{d+\alpha+1}{2}\right)}{\Gamma\left(\frac{d}{2}\right)}\right)$
$\displaystyle\geq
2^{\frac{d}{2}-1}\eta^{\frac{d}{2}}\Gamma\left(\frac{d}{2}\right)\left(1-a(\eta
d)^{\frac{\alpha+1}{2}}\right)\geq\frac{1}{4}(2\eta)^{\frac{d}{2}}\Gamma\left(\frac{d}{2}\right).$

The result (36) then follows from the above inequality and (37).

We now proceed to show that the number of rejections in Algorithm 4 is bounded
from above by a small constant when $\delta$ is properly chosen. In
particular, as shown in Proposition 5.4, it only gets worse by a factor of
$\exp(\delta)$ and the factor does not depend on the dimension $d$. Hence, the
implementation of RGO for $f$ is computationally efficient in practice.

######  Proposition 5.4.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth. If

$\eta\leq\frac{(\alpha+1)^{\frac{2}{\alpha+1}}}{(2L_{\alpha})^{\frac{2}{\alpha+1}}d}\,,$
(43)

then the expected number of iterations in the rejection sampling of Algorithm
4 is at most $2\exp(\delta)$.

Proof: It is a well-known result for rejection sampling that
$X\sim\pi^{X|Y}(x\mid y)$ and the probability that $X$ is accepted is

$\mathbb{P}\left(U\leq\frac{\exp(-f_{y}^{\eta}(X))}{\exp(-h_{1}(X))}\right)=\frac{\int_{\mathbb{R}^{d}}\exp(-f_{y}^{\eta}(x))\mathrm{d}x}{\int_{\mathbb{R}^{d}}\exp(-h_{1}(x))\mathrm{d}x}.$
(44)

If follows directly from the definition of $h_{2}$ in (32b) that

$\displaystyle\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))\mathrm{d}x$
$\displaystyle=\exp(-f_{y}^{\eta}(x^{*}))\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x-x^{*}\|^{2}-\frac{L_{\alpha}}{\alpha+1}\|x-x^{*}\|^{\alpha+1}\right)\mathrm{d}x$

Applying Proposition 5.3 to the above yields

$\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))\mathrm{d}x\geq\exp(-f_{y}^{\eta}(x^{*}))\frac{(2\pi\eta)^{d/2}}{2}.$

Note that the condition (35) in Proposition 5.3 holds thanks to (43). By Lemma
5.2, the above inequality leads to

$\int_{\mathbb{R}^{d}}\exp(-f_{y}^{\eta}(x))\mathrm{d}x\geq\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))\mathrm{d}x\geq\exp(-f_{y}^{\eta}(x^{*}))\frac{(2\pi\eta)^{d/2}}{2}.$
(45)

Using the definition of $h_{1}$ in (32a) and Lemma A.1, we have

$\int_{\mathbb{R}^{d}}\exp(-h_{1}(x))\mathrm{d}x=\exp\left(-f_{y}^{\eta}(\tilde{x}_{J})+\delta\right)(2\pi\eta)^{d/2}.$
(46)

Using (44), (45) and the above identity, we conclude that

$\mathbb{P}\left(U\leq\frac{\exp(-f_{y}^{\eta}(X))}{\exp(-h_{1}(X))}\right)\geq\frac{1}{2}\exp(-f_{y}^{\eta}(x^{*})+f_{y}^{\eta}(\tilde{x}_{J})-\delta)\geq\frac{1}{2}\exp(-\delta),$

and the expected number of the iterations is

$\frac{1}{\mathbb{P}\left(U\leq\frac{\exp(-f_{y}^{\eta}(X))}{\exp(-h_{1}(X))}\right)}\leq
2\exp(\delta).$

We finally bound the total complexity to sample from a log-concave
distribution $\nu$ in (2) with a semi-smooth potential $f$. We combine our
efficient algorithm (Algorithm 4) of RGO for semi-smooth potentials and the
convergent results for ASF, namely Theorem 5.1, to achieve this goal.

######  Theorem 5.5.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth, then Algorithm 1,
initialized with $\rho_{0}^{X}$ and stepsize $\eta\asymp
1/(L_{\alpha}^{\frac{2}{\alpha+1}}d)$, using Algorithm 4 as an RGO has the
iteration-complexity bound

${\cal
O}\left(\frac{L_{\alpha}^{\frac{2}{\alpha+1}}dW_{2}^{2}(\rho_{0}^{X},\nu)}{\varepsilon}\right)$

to achieve $\varepsilon$ error to the target $\nu\propto\exp(-f)$ in terms of
KL divergence. Each RGO requires $\tilde{\cal
O}\left(\frac{1}{d}\left(\frac{1}{\delta}\right)^{\frac{1-\alpha}{\alpha+1}}+1\right)$
subgradient evaluations of $f$ and $2\exp(\delta)$ rejection steps in
expectation. Moreover, if $\nu$ satisfies PI with constant $C_{\rm PI}>0$,
then the iteration-complexity bound to achieve $\varepsilon$ error in terms of
Chi-squared divergence is

$\tilde{\cal O}\left(C_{\rm PI}L_{\alpha}^{\frac{2}{\alpha+1}}d\right).$

Proof: The results follow directly from Theorem 5.1, Theorem 3.5 and
Proposition 5.4 with the choice of stepsize $\eta\asymp
1/(L_{\alpha}^{\frac{2}{\alpha+1}}d)$.

###  5.2 Sampling from composite potentials

In this subsection, we consider sampling from a log-concave distribution
$\nu\propto\exp(-f(x))$ associated with a composite potential $f$ satisfying
(6). This setting is a generalization of the semi-smooth setting studied in
the previous sections. It turns out that both Algorithm 1 and the
implementation for RGO, Algorithm 4, developed for semi-smooth sampling can be
applied directly to this general setting with properly chosen stepsizes.
Below, we extend the analysis in Subsection 5.1 to the composite setting and
establish corresponding complexity results.

The following lemma is a counterpart of Lemma 5.2 in the composite setting.
Its proof is given in Appendix B.

######  Lemma 5.6.

Assume $f$ is convex and satisfies (6). Define

$h_{2}(x):=\frac{1}{2\eta}\|x-x^{*}\|^{2}+\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|x-x^{*}\|^{\alpha_{i}+1}+f_{y}^{\eta}(x^{*}).$
(47)

Then, $h_{2}(x)\geq f_{y}^{\eta}(x)$ for every $x\in\mathbb{R}^{d}$.

The next result is an analogue of the modified Gaussian integral in
Proposition 5.3.

######  Proposition 5.7.

Let $\alpha_{i}\in[0,1]$, $a_{i}\geq 0$, $\eta>0$, and $d\geq 1$. If

$\eta d\sum_{i=1}^{n}a_{i}^{\frac{2}{\alpha_{i}+1}}\leq 1,$ (48)

then

$\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}-\sum_{i=1}^{n}a_{i}\|x\|^{\alpha_{i}+1}\right)\mathrm{d}x\geq(2\pi\eta)^{\frac{d}{2}}\exp\left(-\frac{1}{2}+\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right).$
(49)

Proof: Using the Young’s inequality $st\leq s^{p}/p+t^{q}/q$ with

$s=a2^{\frac{1-\alpha}{2}}\|x\|^{\alpha+1},\quad
t=\frac{1}{2^{\frac{1-\alpha}{2}}},\quad p=\frac{2}{\alpha+1},\quad
q=\frac{2}{1-\alpha},$

we obtain

$a\|x\|^{\alpha+1}\leq(\alpha+1)a^{\frac{2}{\alpha+1}}2^{\frac{-2\alpha}{\alpha+1}}\|x\|^{2}+\frac{1-\alpha}{4}\leq
a^{\frac{2}{\alpha+1}}\|x\|^{2}+\frac{1-\alpha}{4},$

where the second inequality is due to the fact that for $\alpha\in[0,1]$,
$(\alpha+1)2^{\frac{-2\alpha}{\alpha+1}}\leq 1$. Hence, the above inequality
generalizes to

$\sum_{i=1}^{n}a_{i}\|x\|^{\alpha_{i}+1}\leq\sum_{i=1}^{n}a_{i}^{\frac{2}{\alpha_{i}+1}}\|x\|^{2}+\sum_{i=1}^{n}\frac{1-\alpha_{i}}{4}.$

This inequality and Lemma A.1 imply that

$\displaystyle\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}-\sum_{i=1}^{n}a_{i}\|x\|^{\alpha_{i}+1}\right)\mathrm{d}x$
$\displaystyle\geq$
$\displaystyle\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}-\sum_{i=1}^{n}a_{i}^{\frac{2}{\alpha_{i}+1}}\|x\|^{2}-\sum_{i=1}^{n}\frac{1-\alpha_{i}}{4}\right)\mathrm{d}x$
$\displaystyle=$
$\displaystyle\exp\left(\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right)\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\tilde{\eta}}\|x\|^{2}\right)\mathrm{d}x$
$\displaystyle=$
$\displaystyle\exp\left(\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right)(2\pi\tilde{\eta})^{\frac{d}{2}}$
(50)

where

$\frac{1}{\tilde{\eta}}=\frac{1}{\eta}+\sum_{i=1}^{n}a_{i}^{\frac{2}{\alpha_{i}+1}}.$
(51)

It follows from (48) that
$\tilde{\eta}\geq\left(1+\frac{1}{d}\right)^{-1}\eta$. Plugging this
inequality into (50), we have

$\displaystyle\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}-a\|x\|^{\alpha+1}\right)\mathrm{d}x\geq$
$\displaystyle(2\pi\eta)^{\frac{d}{2}}\left(1+\frac{1}{d}\right)^{-\frac{d}{2}}\exp\left(\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right)$
$\displaystyle\geq$
$\displaystyle(2\pi\eta)^{\frac{d}{2}}\exp\left(-\frac{1}{2}+\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right),$

where in the second inequality, we use the fact that

$\left(1+\frac{1}{d}\right)^{\frac{d}{2}}\leq\exp\left(\frac{1}{2}\right).$

With Lemma 5.6 and Proposition 5.7 in hand, we can bound the complexity of
Algorithm 4 as follows. The proof is postponed to Appendix B.

######  Proposition 5.8.

If stepsize $\eta$ satisfies

$\eta
d\sum_{i=1}^{n}\left(\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\right)^{\frac{2}{\alpha_{i}+1}}\leq
1,$ (52)

then rejection steps in Algorithm 4 take at most
$\exp\left(\delta+\frac{1}{2}+\frac{\sum_{i=1}^{n}(1-\alpha_{i})}{4}\right)$
iterations in expectation.

Through the above arguments, we show that Algorithm 4 designed for semi-smooth
potentials is equally effective for composite potentials satisfying (6).
Combining Proposition 5.8 and Theorem 3.9 with the convergence results for
ASF, we obtain the following iteration-complexity bounds for sampling from
composite potentials. The proof is similar to that of Theorem 5.5 and is thus
omitted.

######  Theorem 5.9.

Assume $f$ is a convex and satisfies (6). Consider Algorithm 1, initialized
with $\rho_{0}^{X}$ and stepsize $\eta$ satisfies (52), using Algorithm 4 as a
RGO. Each RGO requires $\tilde{\cal O}\left(\frac{1}{d\delta}+1\right)$
subgradient evaluations of $f$ and
$\exp\left(\delta+\frac{1}{2}+\frac{\sum_{i=1}^{n}(1-\alpha_{i})}{4}\right)$
rejection steps in expectation. The total complexity of Algorithm 1 to achieve
$\varepsilon$ error in terms of KL divergence is

${\cal
O}\left(\frac{\sum_{i=1}^{n}\left(\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\right)^{\frac{2}{\alpha_{i}+1}}dW_{2}^{2}(\rho_{0}^{X},\nu)}{\varepsilon}\right).$

Moreover, if $\nu$ satisfies PI with constant $C_{\rm PI}>0$, then total
complexity to achieve $\varepsilon$ error in terms of Chi-squared divergence
is

$\tilde{\cal
O}\left(\sum_{i=1}^{n}\left(\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\right)^{\frac{2}{\alpha_{i}+1}}dC_{\rm
PI}\right).$

##  6 Conclusion

In this paper, we study proximal algorithms for both optimization and sampling
lacking smoothness. We first establish the complexity bounds of the
regularized cutting-plane method for solving proximal subproblem (3), where
$f$ is convex and satisfies either (5) (semi-smooth) or (6) (composite). This
efficient implementation gives an approximate solution to the proximal map in
optimization, which is the core of both proximal optimization and sampling
algorithms.

For optimization, we develop APBM using a novel adaptive stepsize strategy in
the proximal point method and the approximate proximal map to solve each
proximal subproblem. The proposed APBM is a universal method as it does not
requuire any problem-dependent parameters as input.

For sampling, we propose an efficient method based on rejection sampling and
the approximate proximal map to realize the RGO, which is a proximal sampling
oracle. Finally, combining the sampling complexity of RGO and the complexity
bounds of ASF, which is a counterpart of the proximal point method in
sampling, we establish the complexity bounds of the proximal sampling
algorithm in both semi-smooth and composite settings.

This paper provides a unified perspective to study proximal optimization and
sampling algorithms, while many other interesting questions remain open.
First, APBM is only optimal when $\alpha=0$, i.e., $f$ is Lipschitz
continuous. We are interested in developing a universal method that is optimal
for any $\alpha\in[0,1]$. One possible direction is to incorporate the
acceleration technique into both the regularized cutting-plane method and the
PPF. Second, as acceleration methods are widely used in optimization to obtain
optimal performance, accelerated proximal sampling algorithms are less
explored. It is worth investigating a counterpart of the accelerated proximal
point method [31] in sampling. Finally, we develop APBM as a universal method
for non-smooth optimization, it would be equally important to design a
universal method for sampling.

## Acknowledgement

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

## References

  * [1] Dimitris Bertsimas and Santosh Vempala.  Solving convex programs by random walks.  Journal of the ACM (JACM), 51(4):540–556, 2004.
  * [2] Nawaf Bou-Rabee and Martin Hairer.  Nonasymptotic mixing of the MALA algorithm.  IMA Journal of Numerical Analysis, 33(1):80–110, 2013.
  * [3] Sébastien Bubeck, Qijia Jiang, Yin Tat Lee, Yuanzhi Li, and Aaron Sidford.  Near-optimal method for highly smooth convex optimization.  In Conference on Learning Theory, pages 492–507. PMLR, 2019.
  * [4] Antonin Chambolle and Thomas Pock.  A first-order primal-dual algorithm for convex problems with applications to imaging.  Journal of Mathematical Imaging and Vision, 40(1):120–145, 2011\.
  * [5] Yongxin Chen, Sinho Chewi, Adil Salim, and Andre Wibisono.  Improved analysis for a proximal algorithm for sampling.  In Conference on Learning Theory, pages 2984–3014. PMLR, 2022.
  * [6] 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.
  * [7] 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\.
  * [8] 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.
  * [9] 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.
  * [10] Jiaojiao Fan, Bo Yuan, and Yongxin Chen.  Improved dimension dependence of a proximal algorithm for sampling.  In The Thirty Sixth Annual Conference on Learning Theory, pages 1473–1521. PMLR, 2023.
  * [11] Alexander Gasnikov, Pavel Dvurechensky, Eduard Gorbunov, Evgeniya Vorontsova, Daniil Selikhanovych, and César A Uribe.  Optimal tensor methods in smooth convex and uniformly convexoptimization.  In Conference on Learning Theory, pages 1374–1391. PMLR, 2019.
  * [12] Andrew Gelman, John B Carlin, Hal S Stern, David B Dunson, Aki Vehtari, and Donald B Rubin.  Bayesian Data Analysis.  CRC press, 2013.
  * [13] Stuart Geman and Donald Geman.  Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images.  IEEE Transactions on pattern analysis and machine intelligence, (6):721–741, 1984.
  * [14] Sivakanth Gopi, Yin Tat Lee, and Daogao Liu.  Private convex optimization via exponential mechanism.  In Conference on Learning Theory, pages 1948–1989. PMLR, 2022.
  * [15] 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.
  * [16] Bo Jiang, Haoyue Wang, and Shuzhong Zhang.  An optimal high-order tensor method for convex optimization.  Mathematics of Operations Research, 46(4):1390–1412, 2021.
  * [17] Adam Tauman Kalai and Santosh Vempala.  Simulated annealing for convex optimization.  Mathematics of Operations Research, 31(2):253–266, 2006.
  * [18] 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.
  * [19] Weiwei Kong, Jefferson G Melo, and Renato DC Monteiro.  Complexity of a quadratic penalty accelerated inexact proximal point method for solving linearly constrained nonconvex composite programs.  SIAM Journal on Optimization, 29(4):2566–2593, 2019.
  * [20] Weiwei Kong and Renato DC Monteiro.  An accelerated inexact proximal point method for solving nonconvex-concave min-max problems.  SIAM Journal on Optimization, 31(4):2558–2585, 2021.
  * [21] Galina M. Korpelevič.  An extragradient method for finding saddle points and for other problems.  Èkonom. i Mat. Metody, 12(4):747–756, 1976.
  * [22] Werner Krauth.  Statistical mechanics: algorithms and computations, volume 13.  OUP Oxford, 2006.
  * [23] 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.
  * [24] Jiaming Liang and Yongxin Chen.  A proximal algorithm for sampling from non-smooth potentials.  In 2022 Winter Simulation Conference (WSC), pages 3229–3240. IEEE, 2022.
  * [25] Jiaming Liang and Yongxin Chen.  A proximal algorithm for sampling.  Transactions on Machine Learning Research, 2023.
  * [26] Jiaming Liang and Renato D. C. Monteiro.  A doubly accelerated inexact proximal point method for nonconvex composite optimization problems.  Available on arXiv:1811.11378, 2018.
  * [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 solving hybrid convex composite optimization problems.  Mathematics of Operations Research, 2023.
  * [29] Jiaming Liang, Renato D. C. Monteiro, and Honghao Zhang.  Proximal bundle methods for hybrid weakly convex composite optimization problems.  arXiv preprint arXiv:2303.14896, 2023.
  * [30] Bernard Martinet.  Regularisation d’inequations variationelles par approximations successives.  Revue Francaise d’informatique et de Recherche operationelle, 4:154–159, 1970.
  * [31] Renato D. C. Monteiro and B. F. Svaiter.  An accelerated hybrid proximal extragradient method for convex optimization and its implications to second-order methods.  SIAM Journal on Optimization, 23(2):1092–1125, 2013.
  * [32] Renato D. C. Monteiro and B. F. Svaiter.  Iteration-complexity of block-decomposition algorithms and the alternating direction method of multipliers.  SIAM J. Optim., 23(1):475–507, 2013.
  * [33] 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.
  * [34] Radford M Neal.  MCMC using Hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 2(11):2, 2011.
  * [35] Arkadi Nemirovski and David B. Yudin.  Problem complexity and method efficiency in optimization. Wiley, 1983.
  * [36] Arkaddii S Nemirovskii and Yu E Nesterov.  Optimal methods of smooth convex minimization.  USSR Computational Mathematics and Mathematical Physics, 25(2):21–30, 1985.
  * [37] Yu Nesterov.  Universal gradient methods for convex optimization problems.  Mathematical Programming, 152(1-2):381–404, 2015.
  * [38] Neal Parikh and Stephen Boyd.  Proximal algorithms.  Foundations and Trends in optimization, 1(3):127–239, 2014.
  * [39] Giorgio Parisi.  Correlation functions and computer simulations.  Nuclear Physics B, 180(3):378–384, 1981.
  * [40] Gareth O Roberts and Osnat Stramer.  Langevin diffusions and Metropolis-Hastings algorithms.  Methodology and Computing in Applied Probability, 4(4):337–357, 2002.
  * [41] Gareth O Roberts and Richard L Tweedie.  Exponential convergence of Langevin distributions and their discrete approximations.  Bernoulli, pages 341–363, 1996.
  * [42] R. T. Rockafellar.  Augmented Lagrangians and applications of the proximal point algorithm in convex programming.  Math. Oper. Res., 1(2):97–116, 1976.
  * [43] R Tyrrell Rockafellar.  Monotone operators and the proximal point algorithm.  SIAM journal on control and optimization, 14(5):877–898, 1976.
  * [44] Adil Salim and Peter Richtárik.  Primal dual interpretation of the proximal stochastic gradient Langevin algorithm.  Advances in Neural Information Processing Systems, 33:3786–3796, 2020.
  * [45] Ruoqi Shen, Kevin Tian, and Yin Tat Lee.  Composite logconcave sampling with a restricted Gaussian oracle.  Available on arXiv:2006.05976, 2020.
  * [46] 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.
  * [47] JG Wendel.  Note on the gamma function.  The American Mathematical Monthly, 55(9):563–564, 1948.
  * [48] 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.
  * [49] 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.
  * [50] Bo Yuan, Jiaojiao Fan, Jiaming Liang, Andre Wibisono, and Yongxin Chen.  On a class of gibbs sampling over networks.  In The Thirty Sixth Annual Conference on Learning Theory, pages 5754–5780. PMLR, 2023.

##  Appendix A Technical results

This section collects technical results that are useful in the paper.

######  Lemma A.1 (Gaussian integral).

For any $\eta>0$,

$\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x\|^{2}\right)\mathrm{d}x=(2\pi\eta)^{d/2}.$

The following lemma provides both lower and upper bounds on the ratio of gamma
functions. Its proof can be found in [47].

######  Lemma A.2 (Wendel’s double inequality).

For $0<s<1$ and $t>0$, the gamma function defined as in (38) satisfies

$\left(\frac{t}{t+s}\right)^{1-s}\leq\frac{\Gamma(t+s)}{t^{s}\Gamma(t)}\leq
1,$

or equivalently,

$t^{1-s}\leq\frac{\Gamma(t+1)}{\Gamma(t+s)}\leq(t+s)^{1-s}.$ (53)

######  Lemma A.3.

Assume $f$ is convex and $L_{\alpha}$-semi-smooth (i.e., satisfying (5), then
(7) holds for every $u,v\in\mathbb{R}^{d}$. Assume $f$ is convex and satisfies
(6), then (8) holds for every $u,v\in\mathbb{R}^{d}$.

Proof: We first consider the case when $f$ is convex and $L_{\alpha}$-semi-
smooth. It is easy to see that

$\displaystyle f(u)$ $\displaystyle=f(v)+\int_{0}^{1}\langle
f^{\prime}(v+\tau(v-u)),u-v\rangle\mathrm{d}\tau$ $\displaystyle=f(v)+\langle
f^{\prime}(v),u-v\rangle+\int_{0}^{1}\langle
f^{\prime}(v+\tau(v-u))-f^{\prime}(v),u-v\rangle\mathrm{d}\tau.$

Using the above identity, the Cauchy-Schwarz inequality, and (5), we have

$\displaystyle f(u)-f(v)-\langle f^{\prime}(v),u-v\rangle$
$\displaystyle=\int_{0}^{1}\langle
f^{\prime}(v+\tau(v-u))-f^{\prime}(v),u-v\rangle\mathrm{d}\tau$
$\displaystyle\leq\int_{0}^{1}\left\|f^{\prime}(v+\tau(v-u))-f^{\prime}(v)\right\|\|u-v\|\mathrm{d}\tau$
$\displaystyle\leq\int_{0}^{1}L_{\alpha}\tau^{\alpha}\|u-v\|^{\alpha+1}\mathrm{d}\tau=\frac{L_{\alpha}}{\alpha+1}\|u-v\|^{\alpha+1}.$

Hence, (7) holds. More generally, if $f$ satisfies (6), then (8) follows the
same argument.

##  Appendix B Missing proofs in Subsection 5.2

Proof of Lemma 5.6: It follows from the same argument as in the proof of Lemma
5.2 that (34) holds. Using (34), Lemma 3.1(e), and (8) with $(u,v)=(x,x^{*})$,
we conclude that

$f_{y}^{\eta}(x)-f_{y}^{\eta}(x^{*})\leq\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|x-x^{*}\|^{\alpha_{i}+1}+\frac{1}{2\eta}\|x-x^{*}\|^{2}.$

The lemma immediately follows from the above inequality and the definition of
$h_{2}$ in (47).

Proof of Proposition 5.8: If follows directly from the definition of $h_{2}$
in (47) that

$\displaystyle\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))\mathrm{d}x$
$\displaystyle=\exp(-f_{y}^{\eta}(x^{*}))\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x-x^{*}\|^{2}-\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|x-x^{*}\|^{\alpha_{i}+1}\right)\mathrm{d}x.$

It is easy to see that (52) implies that (48) holds with
$a_{i}=\frac{L_{\alpha_{i}}}{\alpha_{i}+1}$. Hence, by Proposition 5.7, we
have (5.7) holds with $a_{i}=\frac{L_{\alpha_{i}}}{\alpha_{i}+1}$, i.e.,

$\int_{\mathbb{R}^{d}}\exp\left(-\frac{1}{2\eta}\|x-x^{*}\|^{2}-\sum_{i=1}^{n}\frac{L_{\alpha_{i}}}{\alpha_{i}+1}\|x-x^{*}\|^{\alpha_{i}+1}\right)\mathrm{d}x\geq(2\pi\eta)^{\frac{d}{2}}\exp\left(-\frac{1}{2}+\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right).$

The above two inequalities and Lemma 5.6 imply that

$\int_{\mathbb{R}^{d}}\exp(-f_{y}^{\eta}(x))\mathrm{d}x\geq\int_{\mathbb{R}^{d}}\exp(-h_{2}(x))\mathrm{d}x\geq(2\pi\eta)^{\frac{d}{2}}\exp\left(-f_{y}^{\eta}(x^{*})-\frac{1}{2}+\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right).$

As in the proof of Proposition 5.4, (44) and (46) hold. Using (44), (46), and
the above inequality, we have

$\mathbb{P}\left(U\leq\frac{\exp(-f_{y}^{\eta}(X))}{\exp(-h_{1}(X))}\right)\geq\exp\left(f_{y}^{\eta}(\tilde{x}_{J})-f_{y}^{\eta}(x^{*})-\delta-\frac{1}{2}+\frac{\sum_{i=1}^{n}(\alpha_{i}-1)}{4}\right).$

The above inequality and the fact that $f_{y}^{\eta}(\tilde{x}_{J})\geq
f_{y}^{\eta}(x^{*})$ immediately imply that

$\frac{1}{\mathbb{P}\left(U\leq\frac{\exp(-f_{y}^{\eta}(X))}{\exp(-h_{1}(X))}\right)}\leq\exp\left(\delta+\frac{1}{2}+\frac{\sum_{i=1}^{n}(1-\alpha_{i})}{4}\right).$