11affiliationtext: Computer Science Dept.
Purdue University
West Lafayette, Indiana, USA 22affiliationtext: Computer Science Dept.
Indiana University
Bloomington, Indiana, USA 33affiliationtext: Institute for Artificial
Intelligence
Peking University
Beijing, China

# Provable Constrained Stochastic Convex Optimization
with XOR-Projected Gradient Descent

Fan Ding Yijie Wang  Jianzhu Ma  Yexiang Xue

###### Abstract

Provably solving stochastic convex optimization problems with constraints is
essential for various problems in science, business, and statistics. Recently
proposed XOR-Stochastic Gradient Descent (XOR-SGD) provides a convergence rate
guarantee solving the constraints-free version of the problem by leveraging
XOR-Sampling. However, the task becomes more difficult when additional
equality and inequality constraints are needed to be satisfied. Here we
propose XOR-PGD, a novel algorithm based on Projected Gradient Descent (PGD)
coupled with the XOR sampler, which is guaranteed to solve the constrained
stochastic convex optimization problem still in linear convergence rate by
choosing proper step size. We show on both synthetic stochastic inventory
management and real-world road network design problems that the rate of
constraints satisfaction of the solutions optimized by XOR-PGD is $10\%$ more
than the competing approaches in very large searching space. The improved XOR-
PGD algorithm is demonstrated to be more accurate and efficient than both XOR-
SGD and SGD coupled with MCMC based samplers. It is also shown to be more
scalable with respect to the number of samples and processor cores via
experiments with large dimensions.

##  1 INTRODUCTION

Stochastic Convex Optimization problem is of great importance given its wide
applicability in finance, control, robotics, management science, operational
research, ecology, and conservation [Sodomka et al., 2007, Ziukov, 2016, Gomes
et al., 2019]. Advancements made to address this problem have ramifications in
many domains. Mathematically, it minimizes an objective in expectation across
multiple probabilistic scenarios under uncertainty.

$\displaystyle\begin{split}&\min_{x}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathbb{E}_{\theta\sim Pr(\theta)}\leavevmode\nobreak\ \leavevmode\nobreak\
f(x,\theta),\\\ &s.t.\ \ x\in C:=\left\\{\forall i,\leavevmode\nobreak\
h_{i}(x)=0;\leavevmode\nobreak\ \forall j,\leavevmode\nobreak\ g_{j}(x)\leq
0\right\\}.\end{split}$ (1)

where each $f(x,\theta)$ is a a convex function with respect to $x$, $C$ is
the set of constraints, $g_{j}(x)$ are convex functions and $h_{i}(x)$ are
linear. Variable $\theta$ is sampled from distribution $Pr(\theta)$, which is
represented as a Markov random field (MRF) in this paper. We can see that this
problem is highly intractable ($\\#$P-hard) [Ding and Xue, 2021]. Despite the
intractability of computing the expectation over a general probability
distribution, a common operator in probabilistic inference, the problem is
harder to solve because of the existence of the constraints set $C$. A
recently proposed method, XOR-SGD [Ding and Xue, 2021] first harness XOR-
Sampler [Ermon et al., 2013a] to solve the constraints-free version of this
task with guaranteed convergence rate using Stochastic Gradient Descent (SGD).
However, they cannot provide such a bound once the solutions are in a
constrained space which is more common in the real world. In this paper we
consider the setting of the constrained stochastic convex optimization and
first provide an algorithm, XOR-Projected Gradient Descent (XOR-PGD), that has
a provably convergence rate towards the optimal value. The key of XOR-PGD is
to first draw a representative set of samples from $Pr(\theta)$ which yield an
accurate estimation of the gradient direction, and then adjust the estimation
by projecting the variable $x$ to the constraints space $C$. Notice that even
if PGD is a common method in traditional constrained optimization, it is not
trivial to extend it to a stochastic setting where unbiased estimation is
unavailable.

Like XOR-SGD, our XOR-PGD also leverages XOR-Sampling with a constant
approximation guarantee, which reduces the sampling problem into queries of NP
oracles via hashing and projection. Our key contribution is the extension of
classical convergence analysis of XOR-SGD on constrained stochastic convex
optimization problems, where we show that a constant multiplicative bound on
the expectation of the gradient direction is sufficient to bound the final
result against the true optimum. Our theoretic contribution also does not
depend on the unbiasness of the gradients, which was a necessary condition in
previous analysis.

On learning probabilistic graphical models, stochastic optimization is related
to the Marginal Maximum-a-posterior (MMAP) problem [Xue et al., 2016, Liu and
Ihler, 2013, Marinescu et al., 2014, Mauá and de Campos, 2012, Marinescu et
al., 2015, Domke, 2013]. These problems can be formulated as (albeit non-
convex) stochastic optimization problems. Convergence analysis of gradient
descent has been studied for objectives with or without constraints [Wang et
al., 2013, Dubey et al., 2016, Agarwal et al., 2017, Lee et al., 2015, Ruder,
2016, Jin et al., 2017, Ge et al., 2015, Duchi et al., 2011, Hinton et al.,
2012, Kingma and Ba, 2014, Duchi et al., 2018, Allen-Zhu, 2017, 2018]. The
constrained setting is related to a convex-concave saddle point minimax
optimization problem [Mokhtari et al., 2020, Xie et al., 2020, Wang and Li,
2020] where primal-dual methods are often used [Du and Hu, 2019, Hamedani and
Aybat, 2018]. A fruitful line of Stochastic alternating direction method of
multipliers (ADMM) [Ouyang et al., 2013, Zheng and Kwok, 2016, Liu et al.,
2017] has been proposed to deal with constrained optimization under
uncertainty, yet this uncertainty is represented by a uniform distribution
where the unbiased estimation is available. Similar idea of stochastic convex
optimization is also proposed by Ding et al. [2021] for machine learning,
where they leverage a XOR sampler to estimate the partition function in the
learning process of an energy-based model. However, none of them offer a
theoretical guarantee of convergence rate in the setting of constrained
stochastic convex optimization.

Experimental results reveal that XOR-PGD is effective in optimizing
constrained convex stochastic functions. XOR-PGD outperforms competing solvers
XOR-SGD and those running SGD with either MCMC, BP or BPChain samplers on both
the constrained stochastic inventory management and the constrained stochastic
network design problems on real-world data under various conditions in not
only accuracy and speed, but also the rate of constraints satisfaction. In
particular, 90% solutions obtained by XOR-PGD satisfy the constraints set $C$
even when the searching space is very large in the stochastic inventory
management problem, approximately 10% more that of competing methods. Besides,
the improved XOR-PGD algorithm converges faster than XOR-SGD by accessing 20
less XOR samples in each iteration and is able to find better solutions in the
stochastic network design problem. See the experiments section for more
details.

Notations For function $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$, we call it
$L$-smooth if for all $x,y$ in the convex domain $dom\leavevmode\nobreak\ f$,
$f(y)\leq f(x)+\nabla f(x)^{T}(y-x)+\frac{L}{2}||y-x||^{2}$. Denote $f^{+}(x)$
as the positive part of function $f(x)$. In other words,
$f^{+}(x)=\max\\{f(x),0\\}$. $f^{-}(x)$ is defined similarly. For a random
vector $x$, we define $\mathbb{E}[x]$ as the element-wise expectation and the
total variation $Var(x)=\mathbb{E}[||x||_{2}^{2}]-||\mathbb{E}[x]||_{2}^{2}$
where $||\cdot||_{2}^{2}$ is the square of $l_{2}$ norm.

##  2 PRELIMINARIES

###  2.1 Probabilistic Distribution

The probability distribution $Pr(\theta)$ can be defined in various different
forms. We consider $Pr(\theta)$ as a graphical model specified as a factor
graph with $N=|V|$ discrete random variables $\theta_{i}\in\Theta_{i},i\in V$
where $\Theta_{i}=\\{0,1\\}$. The global random vector
$\theta=[\theta_{1},\theta_{2},\ldots,\theta_{N}]$ takes value in the
cartesian product
$\Theta=\Theta_{1}\times\Theta_{2}\times\cdots\times\Theta_{N}$. We consider a
function over $\theta\in\Theta$ as follows:

$\displaystyle
f(\theta)=\prod_{\alpha\in\mathcal{I}}\phi_{\alpha}(\\{\theta\\}_{\alpha})$
(2)

which factors into potentials
$\phi_{\alpha}:\\{\theta\\}_{\alpha}\rightarrow\mathbb{R}^{+}$, where
$\mathcal{I}$ is the set of all the cliques of the graph,
$\\{\theta\\}_{\alpha}\subseteq V$ is a subset of variables that the factor
$\phi_{\alpha}$ depends on. We consider a normalized distribution
$p(\theta)=\frac{1}{Z}f(\theta)$ where $Z$, the normalization constant, also
known as the partition function, is defined as $Z=\sum_{\theta}f(\theta)$. The
structure of $Pr(\theta)$ or the set $\mathcal{I}$ can be built from domain
knowledge and potential functions can be learned from real-world data. The
focus of this paper is not on how to construct the MRF but is how we solve
problem in Equation 1 in general when $Pr(\theta)$ is given in the form shown
in Equation 2.

###  2.2 XOR-SGD

XOR-SGD outperforms SGD with MCMC based samplers in that it has a constant
bound on the probability of drawn samples and further guarantee a constant
bound between the expectation of the distribution formed by samples and the
expectation of true distribution, without requiring exponentially large number
of samples. These samples are drawn by XOR-Sampling [Ermon et al., 2013b], a
recently proposed sampling scheme with a constant approximation guarantee,
which reduces the sampling problem into queries of NP oracles via hashing and
projection. Indeed, XOR-SGD requires accessing NP-oracle queries at each
iteration. Specifically, XOR-SGD is guaranteed to converge to a solution that
is within a vanishing constant away from the true optimum in linear number of
iterations, which is shown as Theorem 1.

######  Theorem 1.

[Ding and Xue, 2021] Let $\rho,\kappa$ be the constant approximation factor as
in Ding and Xue [2021], function
$f(x,\theta):\mathbb{R}^{d}\times\\{0,1\\}^{n}\rightarrow\mathbb{R}$ be a
$L$-smooth convex function w.r.t. $x$. Denote
$OPT=\min_{x}\mathbb{E}_{\theta\sim Pr(\theta)}f(x,\theta)$ as the global
optimum. Let $\sigma^{2}=\max_{x}\\{Var(\nabla_{x}f(x,\theta))\\}$ and
$\varepsilon^{2}=\max_{x}\\{||\mathbb{E}[\nabla_{x}f(x,\theta)]||_{2}^{2}\\}$.
For any $1\leq\rho\kappa\leq\sqrt{2}$, step size
$t\leq\frac{2-\rho^{2}\kappa^{2}}{L\rho\kappa}$ and sample size $N\geq 1$,
$\overline{x_{K}}$ is the output of XOR-SGD and
$\mbox{obj}=\mathbb{E}_{\theta}[f(\overline{x_{K}},\theta)]$ is the objective
function value at $\overline{x_{K}}$. We have:

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[\mbox{obj}]-OPT$
$\displaystyle\leq\frac{\rho\kappa||x_{0}-x^{*}||_{2}^{2}}{2tK}+\frac{t(\sigma^{2}+\varepsilon^{2})}{N}.$
(3)

Theorem 1 states that in expectation, the difference in terms of the objective
function values between the output of XOR-SGD algorithm $\overline{x_{K}}$ and
the true optimum $OPT$ is bounded by a term that scales inversely proportional
to the number of SGD iterations $K$ and a tail term
$\frac{t(\sigma^{2}+\epsilon^{2})}{N}$. Although hard to compute, both
$\sigma^{2}$ and $\epsilon^{2}$ are from the input and do not depend on the
algorithm. However, although the authors provide a constrained version of XOR-
SGD in their paper, they are not able to give a theoretical guarantee of
convergence rate of this algorithm under the constrained setting, which is
known to be more intractable due to the extra constraints on the solution
space.

##  3 XOR-PROJECTED GRADIENT DESCENT

In this section we propose XOR-PGD, a new method to solve the stochastic
convex optimization problems with extra constraints. XOR-PGD converges to
solutions that are at most a constant away from the true optimum in linear
number of iterations. The detailed procedure of XOR-PGD for is shown in
Algorithm 1. To approximate the gradient
$\nabla_{x}\mathbb{E}_{\theta}f(x_{k},\theta)$ at step $k$, XOR-PGD draws $N$
samples $\theta_{1},\dots,\theta_{N}$ from $Pr(\theta)$ using XOR-Sampling.
Because XOR-Sampling has a failure rate, XOR-PGD repeatedly call XOR-Sampling
until all $N$ samples are obtained successfully (line 4 – 10). Once
$\theta_{1},\dots,\theta_{N}$ are obtained, XOR-PGD uses the empirical mean
$\overline{g_{k}}=\frac{1}{N}\sum_{i=1}^{N}\nabla_{x}f(x_{k},\theta_{i})$ as
an approximation for $\nabla_{x}\mathbb{E}_{\theta}f(x_{k},\theta)$.

Then, by setting the step size $t_{k}$ as $\frac{\rho\kappa}{\mu\kappa}$, we
update $x_{k}$ by first a minus of $t_{k}\overline{g_{k}}$ and then a
projection $\mathbb{P}_{C}(x)$ to project it to to constraints space $C$. Here
we define $\mathbb{P}_{C}(x)$ in Equation 4:

$\displaystyle\mathbb{P}_{C}(x):=\arg\min_{y\in C}\frac{1}{2}\|x-y\|_{2}^{2}$
(4)

Using this projection, we can prove that the output of XOR-PGD in expectation
converges to the true optimum within a small constant distance at a linear
speed w.r.t. the number of iterations $K$.

Input: $f(x,\theta),\mu,w(\theta),K,N,l,b,\delta,P,\alpha,C$

1 Initialize $x_{1}$ for function $f(x,\theta)$

2 for  _$k=1$ to $K-1$_ do

3 $i\leftarrow 1$

4 while  _$i\leq N$_ do

5 $s\leftarrow$ XOR-Sampling($w(\theta),l,b,\delta,P,\alpha$)

6 if  _$s\neq Failure$_ then

7 $\theta_{i}\leftarrow s$

8 $i\leftarrow i+1$

9 end if

10

11 end while

12 Compute
$\overline{g_{k}}\leftarrow\frac{1}{N}\sum_{i=1}^{N}\nabla_{x}f(x_{k},\theta_{i})$

13 Compute $t_{k}\leftarrow\frac{\rho\kappa}{\mu k}$

14 Update $x_{k+1}\leftarrow\arg\min_{y\in
C}\frac{1}{2}\|x_{k}-t_{k}\overline{g_{k}}-y\|_{2}^{2}$

15

16 end for

17${x_{K}}\leftarrow\overline{x_{K}}=\frac{1}{K}\sum_{k=1}^{K}x_{k}$

return ${x_{K}}$

Algorithm 1 XOR-PGD

1 … as in Algorithm 1, except

2 replacing line 12 by: $t_{k}\leftarrow\frac{2\rho\kappa}{\mu(k+1)}$

replacing line 15 by:
${x_{K}}\leftarrow\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}=\frac{2\sum_{k=1}^{K}kx_{k}}{K(K+1)}$

Algorithm 2 Improved XOR-PGD

######  Theorem 2.

Let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a $L$-smooth and $\mu$-strongly
convex function and $x^{*}=\arg\min_{x\in C}f(x)$, where $C$ is a convex
constraints set of $x$. In iteration $k$, $g_{k}$ is the estimated gradient,
i.e., $x_{k+1}=\mathbb{P}_{C}(x_{k}-t_{k}g_{k})$ where
$Var(g_{k})\leq\sigma^{2}$ and $\mathbb{P}_{C}(x):=\arg\min_{y\in
C}\frac{1}{2}\|x-y\|_{2}^{2}$ is the projection of $x$ onto $C$. If
$t_{k}=\frac{c}{\mu k}$ and there exists $1\leq c\leq\sqrt{2}$ s.t.
$\frac{1}{c}[\nabla f(x_{k})]^{+}\leq\mathbb{E}[g_{k}^{+}]\leq c[\nabla
f(x_{k})]^{+}$ and $c[\nabla
f(x_{k})]^{-}\leq\mathbb{E}[g_{k}^{-}]\leq\frac{1}{c}[\nabla f(x_{k})]^{-}$,
then the convergence rate of the XOR-PGD is $O(\frac{\log{K}}{K})$:

$\displaystyle\mathbb{E}[f(\overline{x_{K}})]-f(x^{*})\leq\frac{B}{2\mu
K}(1+\log(K)),$ (5)

where $\overline{x_{K}}=\frac{1}{K}\sum_{k=1}^{K}x_{k}$.

###### Proof.

(Theorem 2) Because orthogonal projections contract distances, we have

$\displaystyle\|x_{k+1}-x^{*}\|^{2}$
$\displaystyle\leq\|x_{k}-t_{k}g_{k}-x^{*}\|^{2}$
$\displaystyle=\|x_{k}-x^{*}\|^{2}+t_{k}^{2}\|g_{k}\|^{2}-2t_{k}\langle
g_{k},x_{k}-x^{*}\rangle$

Taking expectation on both side, we have

$\displaystyle\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}]$
$\displaystyle\leq\mathbb{E}[\|x_{k}-x^{*}\|^{2}]+t_{k}^{2}\mathbb{E}[\|g_{k}\|^{2}]$
$\displaystyle-2t_{k}\langle\mathbb{E}[g_{k}],x_{k}-x^{*}\rangle$

From Lemma 1 in Ding and Xue [2021], we know
$-\langle\mathbb{E}[g_{k}],x_{k}-x^{*}\rangle\leq-\frac{1}{c}\langle\nabla
f(x_{k}),x_{k}-x^{*}\rangle$. Then we have

$\displaystyle\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}]$
$\displaystyle\leq\mathbb{E}[\|x_{k}-x^{*}\|^{2}]+t_{k}^{2}\mathbb{E}[\|g_{k}\|^{2}]$
$\displaystyle-\frac{2t_{k}}{c}\langle\nabla f(x_{k}),x_{k}-x^{*}\rangle$
$\displaystyle\leq\mathbb{E}[\|x_{k}-x^{*}\|^{2}]+t_{k}^{2}\mathbb{E}[\|g_{k}\|^{2}]$
$\displaystyle-\frac{2t_{k}}{c}[f(x_{k})-f(x^{*})+\frac{\mu}{2}\|x_{k}-x^{*}\|^{2}]$

The last inequality is because $f$ is $\mu$-strongly convex. After some
rearranging and taking expectation on both sides, we have

$\displaystyle\mathbb{E}[f(x_{k})]-f(x^{*})\leq\frac{ct_{k}}{2}\mathbb{E}[\|g_{k}\|^{2}]\|]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{c-t_{k}\mu}{2t}\mathbb{E}[\|x_{k}-x^{*}\|^{2}]-\frac{c}{2t_{k}}\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}].$

We know
$\mathbb{E}[\|g_{k}\|^{2}]=\|\mathbb{E}[g_{k}]\|^{2}+\text{Var}(g_{k})$, and
from Lemma 1 in Ding and Xue [2021] we have

$\displaystyle\|\mathbb{E}[g_{k}]\|^{2}\leq$ $\displaystyle c\langle\nabla
f(x_{k}),\mathbb{E}[g_{k}]\rangle$ $\displaystyle\leq c\|\nabla
f(x_{k})\|\|\mathbb{E}[g_{k}]\|\leq cL\|\mathbb{E}[g_{k}]\|$

Therefore, $\|\mathbb{E}[g_{k}]\|^{2}\leq c^{2}L^{2}$ is bounded above and
$\mathbb{E}[\|g_{k}\|^{2}]\leq c^{2}L^{2}+\sigma^{2}$. Apply this bound, we
have

$\displaystyle\mathbb{E}[f(x_{k})]-f(x^{*})$ (6)
$\displaystyle\leq\frac{ct_{k}(c^{2}L^{2}+\sigma^{2})}{2}+\frac{c-t_{k}\mu}{2t_{k}}\mathbb{E}[\|x_{k}-x^{*}\|^{2}]$
$\displaystyle-\frac{c}{2t_{k}}\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}]$

With $t_{k}=\frac{c}{\mu k}$, the above inequality becomes

$\displaystyle\mathbb{E}[f(x_{k})]-f(x^{*})\leq\frac{c^{2}(c^{2}L^{2}+\sigma^{2})}{2\mu
k}$ $\displaystyle+\frac{\mu(k-1)}{2}\mathbb{E}[\|x_{k}-x^{*}\|^{2}]-\frac{\mu
k}{2}\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}].$

Summing the above equations for $k=1,\ldots,K$, we get

$\displaystyle\mathbb{E}[f(\frac{1}{K}\sum_{k=1}^{K}x_{k})]-f(x^{*})\leq\frac{1}{K}\sum_{k=1}^{K}\mathbb{E}[f(x_{k})]-f(x^{*})$
$\displaystyle\leq\frac{B}{2\mu
K}\sum_{k=1}^{K}\frac{1}{k}+\frac{\mu}{2K}[0-K\mathbb{E}[\|x_{K+1}-x^{*}\|^{2}]$
$\displaystyle\leq\frac{B}{2\mu K}(1+\log(K))$

where $B=c^{2}(c^{2}L^{2}+\sigma^{2})$. The first inequality is obtained using
the convexity of $f$; the second inequality is obtained from a telescoping
sum. Let $\overline{x_{K}}=\frac{1}{K}\sum_{k=1}^{K}x_{k}$, the above
inequality can be written as

$\displaystyle\mathbb{E}[f(\overline{x_{K}})]-f(x^{*})\leq\frac{B}{2\mu
K}(1+\log(K))$

Therefore, the convergence rate is $O(\frac{\log(K)}{K})$. This completes the
proof. ∎

Theorem 2 states that by choosing the step size $t_{k}$ inverse proportional
to the step $k$, i.e., $t_{k}=\frac{c}{\mu k}$, our XOR-PGD can converge to
the optimal value in $O(\frac{\log K}{K})$. In practice, we can increase the
sample size from 1 to $N$ to further reduce the variance. By replacing the
objective $f(x)$ in Theorem 2 with $\mathbb{E}_{\theta\sim
p(\theta)}f(x,\theta)$, and noticing
$Var(\overline{g_{k}})=Var_{\theta}(\nabla_{x}f(x,\theta))/N$ due to the
sample size $N$, we have the following Theorem 5.

######  Theorem 3.

Let $\rho,\kappa$ be the constant approximation factor as in Ding and Xue
[2021], function
$f(x,\theta):\mathbb{R}^{d}\times\\{0,1\\}^{n}\rightarrow\mathbb{R}$ be a
$L$-smooth and $\mu$-strongly convex function w.r.t. $x$. Denote
$OPT=\min_{x\in C}\mathbb{E}_{\theta\sim Pr(\theta)}f(x,\theta)$ as the global
optimum in constraints set $C$. Let
$\sigma^{2}=\max_{x}\\{Var(\nabla_{x}f(x,\theta))\\}$ and
$\varepsilon^{2}=\max_{x}\\{||\mathbb{E}[\nabla_{x}f(x,\theta)]||_{2}^{2}\\}$.
For any $1\leq\rho\kappa\leq\sqrt{2}$, step size $t\leq\frac{\rho\kappa}{\mu
k}$ and sample size $N\geq 1$,
$\overline{x_{K}}=\frac{1}{K}\sum_{k=1}^{K}x_{k}$ is the output of XOR-PGD
(Algorithm 1) and $\mbox{obj}=\mathbb{E}_{\theta}[f(\overline{x_{K}},\theta)]$
is the objective function value at $\overline{x_{K}}$. We have:

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[\mbox{obj}]-OPT$
$\displaystyle\leq\left(\frac{\rho^{4}\kappa^{4}L^{2}}{2\mu}+\frac{\rho\kappa(\sigma^{2}+\varepsilon^{2})}{2\mu
N}\right)\frac{1+\log(K)}{K}.$

###### Proof.

(Theorem 5) Since we use $N$ samples at each iteration, we have
$\overline{g_{k}}=\frac{1}{N}\sum_{i=1}^{N}g_{k}^{i}$ and
$\mathbb{E}[\overline{g_{k}}]=\mathbb{E}[g_{k}^{i}]$. In each iteration $k$ we
can adjust the parameters in XOR-Sampling to make the tail
$\epsilon\eta_{\phi}$ zero, then for each sample $g_{k}^{i}$ we can obtain
from Theorem 2 that

$\displaystyle\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}$
$\displaystyle\leq\mathbb{E}[g_{k}^{i+}]\leq\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}.$ (7) $\displaystyle\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}$
$\displaystyle\leq\mathbb{E}[g_{k}^{i-}]\leq\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}.$ (8)

The variance of each sample $g_{k}^{i}$ can also be bounded by

$\displaystyle Var(g_{k}^{i})$ $\displaystyle=\mathbb{E}_{\theta^{\prime}\sim
p^{\prime}(\theta^{\prime})}[||\nabla
f(x_{k},\theta^{\prime})||_{2}^{2}]-||\mathbb{E}_{\theta^{\prime}\sim
p^{\prime}(\theta^{\prime})}[\nabla f(x_{k},\theta^{\prime})]||_{2}^{2},$
$\displaystyle\leq\rho\kappa\mathbb{E}_{\theta\sim p(\theta)}[||\nabla
f(x_{k},\theta)||_{2}^{2}],$ $\displaystyle=\rho\kappa(Var(\nabla
f(x_{k},\theta))+||\mathbb{E}_{\theta\sim p(\theta)}[\nabla
f(x_{k},\theta)]||_{2}^{2}),$
$\displaystyle\leq\rho\kappa(\sigma^{2}+\varepsilon^{2}).$

Denote $\overline{g_{k}}^{+}=\max\\{\overline{g_{k}},\textbf{0}\\}$ and
$\overline{g_{k}}^{-}=\min\\{\overline{g_{k}},\textbf{0}\\}$. Clearly,
$g_{k}^{i+}\geq 0$ and $g_{k}^{i-}\leq 0$. Moreover, for a given dimension,
either $g_{k}^{i+}=0$ for that dimension or $g_{k}^{i-}=0$. Evaluating
$\overline{g_{k}}$ dimension by dimension, we can see that
$\overline{g_{k}}^{+}=\frac{1}{N}\sum_{i=1}^{N}g_{k}^{i+}$ and
$\overline{g_{k}}^{-}=\frac{1}{N}\sum_{i=1}^{N}g_{k}^{i-}$. Combined with
Equation 14 and 15, we know

$\displaystyle\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}$
$\displaystyle\leq\mathbb{E}[\overline{g_{k}}^{+}]\leq\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}.$ $\displaystyle\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}$
$\displaystyle\leq\mathbb{E}[\overline{g_{k}}^{-}]\leq\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}.$

Because $\mathbb{E}[\overline{g_{k}}]=\mathbb{E}[g_{k}^{i}]$, we also have

$\displaystyle
Var(\overline{g_{k}})=\frac{1}{N^{2}}Var(\sum_{i=1}^{N}g_{k}^{i})=\frac{Var(g_{k}^{i})}{N}.$

Then the variance of $\overline{g_{k}}$ can be bounded as

$\displaystyle Var(\overline{g_{k}})$
$\displaystyle\leq\frac{\rho\kappa(\sigma^{2}+\varepsilon^{2})}{N}.$

Therefore, we can then apply Theorem 5 to get the result in equation 5.

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[\mathbb{E}_{\theta}[f(\overline{x_{K}},\theta)]]-\mathbb{E}_{\theta}[f(x^{*},\theta)]$
$\displaystyle\leq\frac{\rho^{4}\kappa^{4}L^{2}+\rho^{2}\kappa^{2}\max_{k}\\{Var(\overline{g_{k}})\\}}{2\mu
K}(1+\log(K))$
$\displaystyle\leq\frac{\rho^{4}\kappa^{4}L^{2}+\rho^{2}\kappa^{2}\frac{\sigma^{2}+\varepsilon^{2}}{N}}{2\mu
K}(1+\log(K))$

which can also be written as

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[obj]-OPT$
$\displaystyle\leq\left(\frac{\rho^{4}\kappa^{4}L^{2}}{2\mu}+\frac{\rho\kappa(\sigma^{2}+\varepsilon^{2})}{2\mu
N}\right)\frac{1+\log(K)}{K}.$

This completes the proof. ∎

In Theorem 5 we can see the difference to the optimum is inversely
proportional to both $K$ and $N$ to some extent. To tighten the bound with
fixed number of iterations $K$, we can either conduct more XOR-Sampling scheme
leading to smaller $\rho\kappa$, or generate more samples at each iteration to
reduce the variance term.

###  3.1 Improved XOR-PGD

Although we prove that the proposed XOR-PGD has a theoretical guarantee on the
convergence rate, this convergence rate of $O(\frac{\log K}{K})$ is not rather
surprising. Therefore, we further improve the algorithm by selecting a more
reasonable step size. By selecting $t_{k}=\frac{2\rho\kappa}{\mu(k+1)}$ as in
Algorithm 2, and compute the last output from a weigthed average
${x_{K}}\leftarrow\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}=\frac{2\sum_{k=1}^{K}kx_{k}}{K(K+1)}$,
we develop a improved algorithm of XOR-PGD named (I)XOR-PGD in Algorithm 2. We
prove that this (I)XOR-PGD is able to further accelerate the convergence rate
to $O(\frac{1}{K})$.

######  Theorem 4.

Let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be a $L$-smooth and $\mu$-strongly
convex function and $x^{*}=\arg\min_{x\in C}f(x)$, where $C$ is a convex set
of $x$. In iteration $k$, $g_{k}$ is the estimated gradient, i.e.,
$x_{k+1}=\mathbb{P}_{C}(x_{k}-t_{k}g_{k})$ where $Var(g_{k})\leq\sigma^{2}$
and $\mathbb{P}_{C}(x):=\arg\min_{y\in C}\frac{1}{2}\|x-y\|_{2}^{2}$ is the
projection of $x$ onto $C$. If $t_{k}=\frac{2c}{\mu(k+1)}$ and there exists
$1\leq c\leq\sqrt{2}$ s.t. $\frac{1}{c}[\nabla
f(x_{k})]^{+}\leq\mathbb{E}[g_{k}^{+}]\leq c[\nabla f(x_{k})]^{+}$ and
$c[\nabla f(x_{k})]^{-}\leq\mathbb{E}[g_{k}^{-}]\leq\frac{1}{c}[\nabla
f(x_{k})]^{-}$, then the convergence rate of the Improved XOR-PGD algorithm is
$O(\frac{1}{K})$:

$\displaystyle\mathbb{E}\left[f\left(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}\right)\right]-f(x^{*})\leq\frac{2B}{\mu(K+1)}$

where
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}=\frac{2}{K(K+1)}\sum_{k=1}^{K}kx_{k}$.

###### Proof.

With $t_{k}=\frac{2c}{\mu(k+1)}$ and multiplying (6) by $k+1$, we have

$\displaystyle k(\mathbb{E}[f(x_{k})]-f(x^{*}))\leq\frac{kB}{\mu(k+1)}$
$\displaystyle+\frac{\mu}{4}\left(k(k-1)\mathbb{E}[\|x_{k}-x^{*}\|^{2}]-k(k+1)\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}]\right)$
$\displaystyle\leq\frac{B}{\mu}+\frac{\mu}{4}(k(k-1)\mathbb{E}[\|x_{k}-x^{*}\|^{2}]$
$\displaystyle\qquad\leavevmode\nobreak\
-k(k+1)\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}])$

Summing from $k=1$ to $k=K$, we have

$\displaystyle\sum_{k=1}^{K}k(\mathbb{E}[f(x_{k})]-f(x^{*}))$
$\displaystyle\leq\frac{KB}{\mu}+\frac{\mu}{4}(0-K(K+1)\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}]).$

Then, we have

$\displaystyle\mathbb{E}\left[f\left(\frac{2}{K(K+1)}\sum_{k=1}^{K}kx_{k}\right)\right]-f(x^{*})+\frac{\mu}{2}\mathbb{E}[\|x_{k+1}-x^{*}\|^{2}]$
$\displaystyle\leq\frac{2B}{\mu(K+1)},$

which implies

$\displaystyle\mathbb{E}\left[f\left(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}\right)\right]-f(x^{*})\leq\frac{2B}{\mu(K+1)}$

where
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}=\frac{2}{K(K+1)}\sum_{k=1}^{K}kx_{k}$.
Here we prove the $O(1/K)$ convergence rate. ∎

Similarly, we can get the final convergence rate of the Improved XOR-PGD as in
Algorithm 2 by increasing the sample size from 1 to $N$.

######  Corollary 1.

(Main) Let $\rho,\kappa$ be as before, function
$f(x,\theta):\mathbb{R}^{d}\times\\{0,1\\}^{n}\rightarrow\mathbb{R}$ be a
$L$-Lipschitz and $\mu$-strongly convex function w.r.t. $x$. Denote
$OPT=\min_{x\in C}\mathbb{E}_{\theta\sim Pr(\theta)}f(x,\theta)$ as the global
optimum. Let $\sigma^{2}=\max_{x}\\{Var(\nabla_{x}f(x,\theta))\\}$ and
$\varepsilon^{2}=\max_{x}\\{||\mathbb{E}[\nabla_{x}f(x,\theta)]||_{2}^{2}\\}$.
For any $1\leq\rho\kappa\leq\sqrt{2}$, step size
$t\leq\frac{2\rho\kappa}{\mu(k+1)}$ and sample size $N\geq 1$,
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}=\frac{2}{K(K+1)}\sum_{k=1}^{K}kx_{k}$
is the output of Improved XOR-PGD (algorithm 2) and
$\mbox{obj}=\mathbb{E}_{\theta}[f(\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox},\theta)]$
is the objective function value at
$\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{4.30554pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}$.
We have:

$\displaystyle\mathbb{E}_{\savestack{\tmpbox}{\stretchto{\scaleto{\scalerel*[width("x_{K}")]{\kern-0.6pt\bigwedge\kern-0.6pt}{\rule[-505.89pt]{3.01389pt}{505.89pt}}}{}}{0.5ex}}\stackon[1pt]{x_{K}}{\tmpbox}}[\mbox{obj}]-OP$
$\displaystyle\leq\left(\frac{\rho^{4}\kappa^{4}L^{2}}{\mu}+\frac{\rho\kappa(\sigma^{2}+\varepsilon^{2})}{\mu
N}\right)\frac{1}{K+1}.$

We skip the proof since it is very similar to the proof of Theorem 5. It
should be noticed that although hard to compute, $\sigma^{2}$ and
$\varepsilon^{2}$ are from the input which do not depend on the algorithm. In
addition, it should be noticed that the convergence rate of our algorithm is
determined by the approximation constant $\rho\kappa$ from XOR-sampling. The
smaller this constant is, the lower the success rate of XOR-Sampling has. By
setting proper parameter values, we can get $\rho\kappa=\sqrt{2}$. As a
consequence, we can collect $N$ samples successfully by running XOR-Sampling
around $40N$ times. The time complexity can be further reduced via parallel
sampling. In addition, the sampling procedure is independent of the
optimization step since $Pr(\theta)$ does not depend on $x$.

##  4 EXPERIMENTS

We evaluate our (Improved) XOR-PGD algorithm on the same benchmark of XOR-SGD
[Ding and Xue, 2021], the inventory management (Ziukov [2016], Shapiro and
Philpott [2007]) and the network design problems (Sheldon et al. [2012], Wu et
al. [2017, 2016]). For comparison, we consider XOR-SGD, and also those that
use SGD methods with Gibbs Sampling, Belief Propagation (BP) (Yedidia et al.
[2001], Murphy et al. [2013]), or Belief Propagation Chain (BPChain) (Fan and
Xue [2020]). Similar to the setting in Ding and Xue [2021], for each setting
of both applications, to produce a sample, we let Gibbs sampling have 100 burn
in samples, and then draws one sample every 30 steps. We fix the number of
iteration steps of both BP and BPChain as $20$, which is enough for belief
propagation to converge. We allow SGD with Gibbs sampling, BP and BPChain to
draw more samples than both XOR-PGD and XOR-SGD for a fair comparison. For
both applications, we use MRF as probabilistic models for $Pr(\theta)$. All
experiments were conducted using single core architectures on Intel Xeon Gold
6126 2.60GHz machines with 96GB RAM and a wall-time limit of $10$ hours. We
use IBM ILOG Cplex 12.71 as the solver of NP oracle to produce each XOR sample
for both XOR-PGD and XOR-SGD. Notice that the projection step in XOR-PGD is
also sovled by Cplex. Once a solution $x$ is generated by either algorithm, we
use an exact weighted counter ACE Barton et al. [2016] to evaluate
$\mathbb{E}_{\theta\sim Pr(\theta)}f(x,\theta)$ exactly. All objective values
reported are from ACE. Since we run each algorithm on one single core with a
wall-time limit of $10$ hours for a fair comparison thus not all algorithms
can complete all iterations, we report the best results found by each
algorithm within the time limit.

| Number of materials
---|---
| 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | 100
savings/Gibbs | 8.57% | 19.29% | 13.43% | 13.45% | 14.48% | 13.46% | 11.6% | 13.87% | 17.31% | 20.32%
savings/BP | 22.89% | 20.7% | 17.21% | 18.63% | 19.52% | 16.94% | 14.94% | 25.16% | 30.7% | 34.47%
savings/BPChain | 6.57% | 13.74% | 12.17% | 12.01% | 9.19% | 14.3% | 12.71% | 19.2% | 24.89% | 30.83%
savings/XOR-SGD | 1.54% | 1.74% | 5.31% | 2.86% | 1.21% | 2.64% | 2.24% | 2.87% | 5.38% | 10.27%
Table 1: The percentage of savings of the solutions found by Improved XOR-PGD ((I)XOR-SGD) in Algorithm 2 against other methods on 100% storage limit varying the number of materials. We can see the Improved XOR-PGD on average saves 10% cost against the first three methods and 2% against XOR-SGD, a naive method coupled with XOR-Sampler. | Number of materials
---|---
| 10 | 20 | 30 | 40 | 50 | 60 | 70 | 80 | 90 | 100
XOR-SGD | 98.3% | 97.3% | 95.3% | 94.1% | 93.4% | 91.6% | 89.6% | 87.3% | 83.4% | 80.9%
XOR-PGD | 99.4% | 99.3% | 98.5% | 98.1% | 97.4% | 96.6% | 95.2% | 93.4% | 92.9% | 91.3%
(I)XOR-PGD | 99.6% | 99.4% | 98.3% | 97.6% | 97.5% | 96.4% | 95.9% | 93.5% | 91.5% | 91.7%
Table 2: The average rate of constraints satisfaction of 1000 solutions found
by both XOR-PGD and (I)XOR-SGD against XOR-SGD on 100% storage limit and with
60 samples each iteration varying the number of materials. We can see
solutions found by both of our methods have more than $90\%$ rate of
constraints satisfaction even when the searching space is very large, which is
$10\%$ more than that from XOR-SGD.

Figure 1:  The upper row is experimental results on the inventory management
problem and the bottom row is the results on the network design problem. The
Improved XOR-PGD ((I)XOR-PGD) is better than the baselines in all cases.
(Upper left) The objective values found by all methods with 50 materials
varying storage limits. our method always finds better solutions than the
other approaches. (Upper right) The objective values found by all methods with
different number of samples for approximation in each iteration with 100%
storage limit and 50 materials. (Bottom left) The percentage of savings of
(I)XOR-PGD against other methods with 100% budget. (I)XOR-PGD saves on average
$>3\%$ commuting time on all 4 benchmarks. (Bottom middle) Commuting time of
the solutions found by (I)XOR-PGD and baselines varying budgets on network
“weak 20”. (I)XOR-PGD always leads to the least commuting time no matter how
much the budget is given. (Bottom right) Commuting time found by the three
methods with different number of samples (100% budget, “weak 20”). Both XOR-
PGD and (I) XOR-PGD outperforms XOR-SGD with less XOR samples.

###  4.1 Stochastic Inventory Management

We first investigate our algorithm on the stochastic inventory management
problem studied in Shapiro and Philpott [2007]. Assuming there are $n$
materials. The demand of material $i$ is $d_{i}$. Let
$d=(d_{1},\ldots,d_{n})^{T}$ be the demand vector. The manager stocks $x_{i}$
amount of material $i$ at the beginning of the season. Each unit of material
$i$ takes storage space $w_{i}$, and the total amount of pre-order is limited
by the available storage space $X$. At the end of the production season,
demand $d$ will be revealed to the the manager. We assume the cost of ordering
the $i$-th material is $c_{i}$ per unit. If the demand $d_{i}>x_{i}$, then a
back order is needed, of which one unit costs $b_{i}\geq c_{i}$. Overall, the
cost for back order is $b_{i}(d_{i}-x_{i})$ if $d_{i}>x_{i}$, and is zero
otherwise. On the other hand, if $d_{i}<x_{i}$, then a holding cost of $h_{i}$
per unit is incurred, leading to an additional total cost
$h_{i}(x_{i}-d_{i})$. Summing it up, the cost for material $i$ is
$G_{i}=c_{i}x_{i}+b_{i}[d_{i}-x_{i}]^{+}+h_{i}[x_{i}-d_{i}]^{+}$ where
$[a]^{+}$ denotes the maximum of $a$ and 0. Then, the total cost will be
$G(x,d)=\sum_{i=1}^{n}G_{i}$. The manager want to minimize his operational
cost, which translates to this problem:

$\displaystyle\min_{x\geq 0}\leavevmode\nobreak\ \leavevmode\nobreak\
\mathbb{E}_{d\sim Pr(d)}[G(x,d)],\quad s.t.\quad w^{T}x\leq X.$ (9)

where $G(x,d)$ is convex w.r.t. $x$. We run the experiments varying the number
of materials $n$, the storage limit, and the number of samples we use in XOR-
PGD and other methods. Parameters and experimental settings are the same as in
Ding and Xue [2021]. In terms of the parameters in XOR-Sampling we fix
$P=100,b=7,\epsilon=0.01$ and the others the same as in Ermon et al. [2013b]
to guarantee $\rho\kappa=\sqrt{2}$. Learning rate $t$ is $0.1$ at first and
divided by $10$ after $50$ iterations, then further divided by $10$ after
$100$ iterations. $\eta$ is 10 at first and divided by $10$ after $50$
iterations, then further divided by $10$ after $100$ iterations. The total
number of both $K$ and $M$ are set to be $200$. However, since we run each
algorithm on one single core with a wall-time limit of $10$ hours for a fair
comparison, not all algorithms can complete all iterations. The plots are
based on the best results found by each algorithm within the time limit.

Table 1 shows that our algorithm outperforms other methods varying number of
materials on the percentage reduction of the objective values of the
solutions. In math form, for example for Gibbs Sampling, the metric is
$(obj(\mbox{Gibbs})-obj(\mbox{(I)XOR-PGD}))/obj(\mbox{Gibbs})$ (metrics for
other approaches are analogous). The objective optimized by (I)XOR-PGD is on
average 2 percent better than XOR-SGD and $10\%$ more than that optimized by
other baselines. In addition, Table 2 shows that the rate of constraints
satisfaction of the solutions from either XOR-PGD or (I) XOR-PGD is $10\%$
more than that of XOR-SGD even in very large searching space. The first row in
Figure 1 shows the objective values of solutions varying the storage limit,
and the objective values varying the number of samples with $100\%$ storage
limit and 50 different materials. We can see from the left figure that with
the storage limit increasing, (I)XOR-PGD is always better than all baselines,
and from the right figure that (I)XOR-PGD found better solutions with less
samples than XOR-SGD. Since Ding and Xue [2021] already reported that XOR-SGD
is faster than the other competing approaches, here by comparing efficient
sample size with XOR-SGD we conclude that (I)XOR-PGD is faster than all the
other methods.

###  4.2 Stochastic Network Design

We then consider the expected commuting time of a random walk in network
optimization Ghosh et al. [2008], McClure et al. [2016], Inman et al. [2013],
which is also a benchmark in Ding and Xue [2021]. Given an undirected graph
$G=(V,E)$, where $|V|=m,|E|=n$. Let $g=(g_{1},\dots,g_{n})^{T}$, and edge
$e\in E$ is associated with a binary random variable $\theta_{e}$ that
describes the state of the edge during disasters. $\theta_{e}=0$ means that
the edge is destroyed, and 1 otherwise. Let
$\theta=(\theta_{1},\dots,\theta_{n})$. It will take money $c_{e}$ to increase
one unit of $g_{e}$ and we have a total budget of $B$. Then, From Ghosh et al.
[2008] we know the commuting time

$\displaystyle\overline{C}(g,\theta)=\frac{4(\mathbf{1}^{T}g)}{(m-1)}\Big{(}\mathbf{Tr}(L+\mathbf{1}\mathbf{1}^{T}/m)^{-1}-1\Big{)}$

which is convex w.r.t. $g$, where $L$ is the weighted Laplacian matrix. We
would like to find the best network improvement plan under the given budget,
which minimizes the expected commuting time averaged over all stochastic
events to maximizes the network connectivity. Mathematically, our problem can
be formulated as the following problem:

$\displaystyle\min_{\Delta g\geq 0}\leavevmode\nobreak\ \mathbb{E}_{\theta\sim
Pr(\theta)}[\overline{C}(g+\Delta g,\theta)],\leavevmode\nobreak\
s.t.\leavevmode\nobreak\ \sum_{e\in E}c_{e}\Delta g_{e}\leq B.$

We evaluate our algorithms on the Flood Preparation problem for the emergency
medical services (EMS) on road networks as studied in Wu et al. [2016] and
test our algorithm on four benchmarks involving the weak and the strong
network originally evaluated in Ding and Xue [2021] varying the percentage of
the largest budget size and samples size, which is shown in the bottom row of
Figure 1. The largest budget size $B$ is $1000$. Total number of SGD
iterations is 2000, while not all algorithms can complete all 2000 iterations
within the time limit of 10 hours. Parameters in XOR-Sampling are set to be
the same as in Ding and Xue [2021].

The results are similar to those for the inventory management problem, where
results clearly show that (I)XOR-PGD outperforms other methods both in
efficiency and in the quality of solutions. On all of the four different type
of networks (I)XOR-PGD saves on average $>3\%$ against XOR-SGD and more than
$9\%$ against the other competing methods. In addition, we would like to
emphasize that (I)XOR-PGD with 40 samples already outperforms XOR-SGD with 60
samples in the right figure while finds better solutions at the same time, and
even naive XOR-PGD finds better solutions than XOR-SGD given the same sample
size.

The left figure in Figure 3 shows the percentage of savings between SGD with
other sampling methods and XOR-SGD among all of the 4 different networks,
while the middle and the right figures show the averaged commuting time with
regard to different budget sizes and different number of samples,
respectively. For the left and the middle figures, we let XOR-SGD take 100
samples in each iteration while SGD with other methods take 10,000. We can see
from the left figure that objective optimized by XOR-SGD is at least 5% better
than that optimized by other methods for all the 4 different networks. In
addition, from the middle and the right figures we know that with the increase
of either budget size or the number of samples, our method can find
consistently better solutions than the compared methods. In particular, from
the right figure we can see even 40 samples in each iteration are enough for
XOR-SGD to compete with the result from Gibbs with 20,000 samples. Meanwhile,
XOR-SGD also runs faster than the compared method under this situation. In
this experiment, XOR-SGD with 40 samples take 1 minutes 40 seconds per SGD
iteration, while SGD with 20,000 Gibbs samples need 2.5 minutes per iteration.
Since sampling time of both BP and BPChain is no shorter than Gibbs Sampling,
we thus conclude that XOR-SGD outperforms other methods both in efficiency and
in the quality of solutions found.

##  5 CONCLUSION

We proposed XOR-PGD, a provable algorithm to attack constrained convex
stochastic optimization problems, which are crucial for many decision-making
applications with uncertainty. We showed theoretically that our algorithm has
a linear convergence rate to the global optimum by chossing proper step sizes.
Empirically, we demonstrated the superior performance of XOR-PGD on both the
stochastic inventory management and the stochastic network design problems. In
particular, 90% solutions obtained by XOR-PGD satisfy the constraints set $C$
even when the searching space is very large in the inventory management
problem, approximately 10% more that of competing methods. Besides, XOR-PGD
converges faster than XOR-SGD by accessing 20% less XOR samples in each
iteration and is able to find better solutions in the stochastic network
design problem. Overall, our paper demonstrates the power of integrating
cutting-edge computer science technology with real-world problems. Our paper
will also stimulate further academic progress in stochastic optimization,
constrained optimization, probabilistic inference with hashing and
randomization, and non-convex optimizations with insights from real-world
applications. Future work includes tightening the constant bound and
accelerating the convergence rate. We will keep active to investigate if our
approach can motivate new algorithms for non-convex stochastic optimization
problems.

## References

  * Agarwal et al. [2017] Naman Agarwal, Zeyuan Allen-Zhu, Brian Bullins, Elad Hazan, and Tengyu Ma.  Finding approximate local minima faster than gradient descent.  In _Proceedings of the 49th Annual ACM SIGACT Symposium on Theory of Computing_ , pages 1195–1199, 2017.
  * Allen-Zhu [2017] Zeyuan Allen-Zhu.  Katyusha: The first direct acceleration of stochastic gradient methods.  _The Journal of Machine Learning Research_ , 18(1):8194–8244, 2017.
  * Allen-Zhu [2018] Zeyuan Allen-Zhu.  Natasha 2: Faster non-convex optimization than sgd.  In _Advances in neural information processing systems_ , pages 2675–2686, 2018.
  * Barton et al. [2016] John P Barton, Eleonora De Leonardis, Alice Coucke, and Simona Cocco.  Ace: adaptive cluster expansion for maximum entropy graphical model inference.  _Bioinformatics_ , 32(20):3089–3097, 2016.
  * Ding and Xue [2021] Fan Ding and Yexiang Xue.  Xor-sgd: Provable convex stochastic optimization for decision-making.  In _Uncertainty in Artificial Intelligence_. UAI, 2021.
  * Ding et al. [2021] Fan Ding, Jianzhu Ma, Jinbo Xu, and Yexiang Xue.  Xor-cd: Linearly convergent constrained structure generation.  In _International conference on machine learning_. ICML, 2021.
  * Domke [2013] Justin Domke.  Learning graphical model parameters with approximate marginal inference.  _IEEE transactions on pattern analysis and machine intelligence_ , 35(10):2454–2467, 2013.
  * Du and Hu [2019] Simon S Du and Wei Hu.  Linear convergence of the primal-dual gradient method for convex-concave saddle point problems without strong convexity.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 196–205. PMLR, 2019.
  * Dubey et al. [2016] Kumar Avinava Dubey, Sashank J Reddi, Sinead A Williamson, Barnabas Poczos, Alexander J Smola, and Eric P Xing.  Variance reduction in stochastic gradient langevin dynamics.  In _NIPS_ , pages 1154–1162, 2016.
  * Duchi et al. [2011] John Duchi, Elad Hazan, and Yoram Singer.  Adaptive subgradient methods for online learning and stochastic optimization.  _Journal of machine learning research_ , 12(Jul):2121–2159, 2011.
  * Duchi et al. [2018] John Duchi, Feng Ruan, and Chulhee Yun.  Minimax bounds on stochastic batched convex optimization.  In _Conference On Learning Theory_ , pages 3065–3162, 2018.
  * Ermon et al. [2013a] Stefano Ermon, Carla P. Gomes, Ashish Sabharwal, and Bart Selman.  Taming the curse of dimensionality: Discrete integration by hashing and optimization.  In _Proceedings of the 30th ICML_ , 2013a.
  * Ermon et al. [2013b] Stefano Ermon, Carla P. Gomes, Ashish Sabharwal, and Bart Selman.  Embed and project: Discrete sampling with universal hashing.  In _Advances in Neural Information Processing Systems (NIPS)_ , 2013b.
  * Fan and Xue [2020] Ding Fan and Yexiang Xue.  Contrastive divergence learning with chained belief propagation.  In _International Conference on Probabilistic Graphical Models_ , 2020\.
  * Ge et al. [2015] Rong Ge, Furong Huang, Chi Jin, and Yang Yuan.  Escaping from saddle points—online stochastic gradient for tensor decomposition.  In _Conference on Learning Theory_ , pages 797–842, 2015.
  * Ghosh et al. [2008] Arpita Ghosh, Stephen Boyd, and Amin Saberi.  Minimizing effective resistance of a graph.  _SIAM review_ , 50(1):37–66, 2008.
  * Gomes et al. [2019] Carla Gomes, Thomas Dietterich, Christopher Barrett, Jon Conrad, Bistra Dilkina, Stefano Ermon, Fei Fang, Andrew Farnsworth, Alan Fern, Xiaoli Fern, et al.  Computational sustainability: Computing for a better world and a sustainable future.  _Communications of the ACM_ , 62(9):56–65, 2019\.
  * Hamedani and Aybat [2018] Erfan Yazdandoost Hamedani and Necdet Serhat Aybat.  A primal-dual algorithm with line search for general convex-concave saddle point problems.  _arXiv preprint arXiv:1803.01401_ , 2018.
  * Hinton et al. [2012] Geoffrey Hinton, Nitish Srivastava, and Kevin Swersky.  Neural networks for machine learning lecture 6a overview of mini-batch gradient descent.  2012\.
  * Inman et al. [2013] Robert M. Inman, Brent L. Brock, Kristine H. Inman, Shawn S. Sartorius, Bryan C. Aber, Brian Giddings, Steven L. Cain, Mark L. Orme, Jay A. Fredrick, Bob J. Oakleaf, Kurt Alt, Eric A. Odell, and Guillaume Chapron.  Developing priorities for metapopulation conservation at the landscape scale: Wolverines in the western united states.  2013\.
  * Jin et al. [2017] Chi Jin, Praneeth Netrapalli, and Michael I Jordan.  Accelerated gradient descent escapes saddle points faster than gradient descent.  _arXiv preprint arXiv:1711.10456_ , 2017.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ , 2014.
  * Lee et al. [2015] Jason D Lee, Qihang Lin, Tengyu Ma, and Tianbao Yang.  Distributed stochastic variance reduced gradient methods and a lower bound for communication complexity.  _arXiv preprint arXiv:1507.07595_ , 2015.
  * Liu and Ihler [2013] Qiang Liu and Alexander T. Ihler.  Variational algorithms for marginal MAP.  _Journal of Machine Learning Research_ , 14, 2013.
  * Liu et al. [2017] Yuanyuan Liu, Fanhua Shang, and James Cheng.  Accelerated variance reduced stochastic admm.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 31, 2017.
  * Marinescu et al. [2014] Radu Marinescu, Rina Dechter, and Alexander T. Ihler.  AND/OR search for marginal MAP.  In _Proceedings of the Thirtieth Conference on Uncertainty in Artificial Intelligence, UAI_ , 2014.
  * Marinescu et al. [2015] Radu Marinescu, Rina Dechter, and Alexander Ihler.  Pushing forward marginal map with best-first search.  In _Proceedings of the 24th International Conference on Artificial Intelligence (IJCAI)_ , 2015.
  * Mauá and de Campos [2012] Denis Deratani Mauá and Cassio Polpo de Campos.  Anytime marginal MAP inference.  In _Proceedings of the 29th ICML_ , 2012.
  * McClure et al. [2016] Meredith L McClure, Andrew J. Hansen, and Robert M. Inman.  Connecting models to movements: testing connectivity model predictions against empirical migration and dispersal data.  _Landscape Ecology_ , 31:1419–1432, 2016.
  * Mokhtari et al. [2020] Aryan Mokhtari, Asuman E Ozdaglar, and Sarath Pattathil.  Convergence rate of o(1/k) for optimistic gradient and extragradient methods in smooth convex-concave saddle point problems.  _SIAM Journal on Optimization_ , 30(4):3230–3251, 2020.
  * Murphy et al. [2013] Kevin Murphy, Yair Weiss, and Michael I Jordan.  Loopy belief propagation for approximate inference: An empirical study.  _arXiv preprint arXiv:1301.6725_ , 2013.
  * Ouyang et al. [2013] Hua Ouyang, Niao He, Long Tran, and Alexander Gray.  Stochastic alternating direction method of multipliers.  In _International Conference on Machine Learning_ , pages 80–88. PMLR, 2013.
  * Ruder [2016] Sebastian Ruder.  An overview of gradient descent optimization algorithms.  _arXiv preprint arXiv:1609.04747_ , 2016.
  * Shapiro and Philpott [2007] Alexander Shapiro and Andy Philpott.  A tutorial on stochastic programming.  2007\.
  * Sheldon et al. [2012] Daniel Sheldon, Bistra Dilkina, Adam N Elmachtoub, Ryan Finseth, Ashish Sabharwal, Jon Conrad, Carla P Gomes, David Shmoys, William Allen, Ole Amundsen, et al.  Maximizing the spread of cascades using network design.  _arXiv preprint arXiv:1203.3514_ , 2012.
  * Sodomka et al. [2007] Eric Sodomka, John Collins, and Maria Gini.  Efficient statistical methods for evaluating trading agent performance.  2007\.
  * Wang et al. [2013] Chong Wang, Xi Chen, Alexander J Smola, and Eric P Xing.  Variance reduction for stochastic gradient optimization.  In _Advances in Neural Information Processing Systems_ , pages 181–189, 2013.
  * Wang and Li [2020] Yuanhao Wang and Jian Li.  Improved algorithms for convex-concave minimax optimization.  _arXiv preprint arXiv:2006.06359_ , 2020.
  * Wu et al. [2016] Xiaojian Wu, Daniel R Sheldon, and Shlomo Zilberstein.  Optimizing resilience in large scale networks.  In _Proceedings of the 30th Conference of AAAI_ , 2016.
  * Wu et al. [2017] Xiaojian Wu, Yexiang Xue, Bart Selman, and Carla P. Gomes.  Xor-sampling for network design with correlated stochastic events.  In _Proceedings of the 26th IJCAI_ , pages 4640–4647, 2017.
  * Xie et al. [2020] Guangzeng Xie, Luo Luo, Yijiang Lian, and Zhihua Zhang.  Lower complexity bounds for finite-sum convex-concave minimax optimization problems.  In _International Conference on Machine Learning_ , pages 10504–10513. PMLR, 2020.
  * Xue et al. [2016] Yexiang Xue, Zhiyuan Li, Stefano Ermon, Carla P. Gomes, and Bart Selman.  Solving marginal map problems with np oracles and parity constraints.  In _Proceedings of the 29th Annual Conference on NIPS_ , 2016.
  * Yedidia et al. [2001] Jonathan S Yedidia, William T Freeman, and Yair Weiss.  Generalized belief propagation.  In _Advances in neural information processing systems_ , pages 689–695, 2001.
  * Zheng and Kwok [2016] Shuai Zheng and James T Kwok.  Fast-and-light stochastic admm.  In _IJCAI_ , pages 2407–2613, 2016.
  * Ziukov [2016] Serhii Ziukov.  A literature review on models of inventory management under uncertainty.  2016\.

## Appendix

##  Appendix A XOR-Sampling for the Weighted Case

The text here provides a synopsis for the approach in Ermon et al. [2013b]. We
still encourage the readers to read the original text for a better
explanation. Let $w(\theta),p(\theta)$ and $Z$ as defined before, the high-
level idea of XOR-Sampling is to first dicretize $w(\theta)$ to
$w^{\prime}(\theta)$ as in Definition 1, followed by embedding the weighted
$w^{\prime}(\theta)$ to the unweighted space $\Delta_{w}$. Finally, XOR-
sampling uses counting based on hashing and randomization to sample uniformly
from $\Delta_{w}$.

######  Definition 1.

Assume $w(\theta)$ has both upper and lower bound, namely,
$M=\max_{\theta}w(\theta)$ and $m=\min_{\theta}w(\theta)$. Let $b\geq
1,\epsilon>0,r=2^{b}/(2^{b}-1)$ and $l=\lceil\log_{r}(2^{n}/\epsilon)\rceil$.
Partition the configurations into the following weight based disjoint buckets:
$\mathcal{B}_{i}=\\{\theta|w(\theta)\in(\frac{M}{r^{i+1}},\frac{M}{r^{i}}]\\},i=0,\ldots,l-1$
and $\mathcal{B}_{l}=\\{\theta|w(\theta)\in(0,\frac{M}{r^{l}}]\\}$. The
discretized weight function
$w^{\prime}:\\{0,1\\}^{n}\rightarrow\mathbb{R}^{+}$ is defined as follows:
$w^{\prime}(\theta)=\frac{M}{r^{i+1}}$ if
$\theta\in\mathcal{B}_{i},i=0,\ldots,l-1$ and $w^{\prime}(\theta)=0$ if
$\theta\in\mathcal{B}_{l}$. This leads to the corresponding discretized
probability distribution $p^{\prime}(\theta)=w^{\prime}(\theta)/Z^{\prime}$
where $Z^{\prime}$ is the normalization constant of $w^{\prime}(\theta)$.

For the weighted case, the goal of XOR-sampling is to guarantee that the
probability of sampling one $\theta$ is proportional to the unnormalized
density (up to a multiplicative constant). By Definition 1, we obtain a
distribution $p^{\prime}(x)$ which satisfying $\frac{1}{\rho}p(x)\leq
p^{\prime}(x)\leq\rho p(x)$ where $\rho=\frac{r^{2}}{1-\epsilon}$. Then, XOR-
sampling implements a horizontal slice technique to transform a weighted
problem into an unweighted one. For the easiness of illustration, we denote
$M^{\prime}=\max_{\theta}w^{\prime}(\theta)$ and $m^{\prime}$ as the smallest
non-zero value of $w^{\prime}(\theta)$. Then consider the simple case where
$b=1$ and $r=2$, where we have $M^{\prime}=2^{l-1}m^{\prime}$. Let
$\delta=(\delta_{0},\ldots,\delta_{l-2})^{T}\in\\{0,1\\}^{l-1}$ be a binary
vector of length $l-1$, XOR-sampling samples $(\theta,\delta)$ uniformly at
random from the following set $\Delta_{w}$ using the unweighted version of
sampling based on hashing and randomization:

$\displaystyle\Delta_{w}=\\{(\theta,\delta):w^{\prime}(\theta)\leq
2^{i+1}m^{\prime}\Rightarrow\delta_{i}=0\\}.$ (10)

If we sample $(\theta,\delta)$ uniformly at random from $\Delta_{w}$ and then
only return $\theta$, it can be proved that the probability of sampling
$\theta$ from $w^{\prime}(\theta)$ is proportional to $m^{\prime}2^{i-1}$ when
$w(\theta)$ is sandwiched between $m^{\prime}2^{i-1}$ and $m^{\prime}2^{i}$.
Therefore, this technique leads to the constant approximation guarantee of
XOR-Sampling. The precise statement of the guarantee is in Theorem 2. For
general case of $b$ and $r$, please refer to Ermon et al. [2013b].

Setting $\epsilon\eta_{\phi}$ to Zero In Definition 1 we can make $b$ larger
and $\epsilon$ smaller enough, then there will be a possibly large but finite
value of $l$ such that $\frac{M}{r^{l}}$ is smaller than $m$ , which leads
$\mathcal{B}_{l}$ to be empty and $\epsilon\eta_{\phi}$ to be zero.

##  Appendix B Proofs

###  B.1 Proof of Lemma 1

We define two functions $g_{k}^{+}=\max\\{g_{k},\mathbf{0}\\}$ and
$g_{k}^{-}=\min\\{g_{k},\mathbf{0}\\}$ where $\mathbf{0}$ is a vector of all 0
which has the same dimension as $g_{k}$. We have $g_{k}=g_{k}^{+}+g_{k}^{-}$.
We define both $\nabla f(x_{k})^{+}$ and $\nabla f(x_{k})^{-}$ in the similar
way. Then Lemma 1 gives the new bounds of two terms assuming the constant
bound on the gradient, which are essential to the proof of convergence rate.
The proof of Lemma 1 is as follows:

###### Proof.

(Lemma 1) Since we have the constant bound that

$\displaystyle\frac{1}{c}\nabla f(x_{k})^{+}$
$\displaystyle\leq\mathbb{E}[g_{k}^{+}]\leq c\nabla f(x_{k})^{+}.$ (11)
$\displaystyle c\nabla f(x_{k})^{-}$
$\displaystyle\leq\mathbb{E}[g_{k}^{-}]\leq\frac{1}{c}\nabla f(x_{k})^{-}.$
(12)

and because of $g_{k}^{+}\geq\mathbf{0}$ and $g_{k}^{-}\leq\mathbf{0}$ we can
obtain

$\displaystyle\frac{1}{c}||\mathbb{E}[g_{k}^{+}]||_{2}^{2}$
$\displaystyle=\frac{1}{c}\langle\mathbb{E}[g_{k}^{+}],\mathbb{E}[g_{k}^{+}]\rangle\leq\langle\nabla
f(x_{k})^{+},\mathbb{E}[g_{k}^{+}]\rangle$ $\displaystyle\leq
c\langle\mathbb{E}[g_{k}^{+}],\mathbb{E}[g_{k}^{+}]\rangle=c||\mathbb{E}[g_{k}^{+}]||_{2}^{2}.$
$\displaystyle\frac{1}{c}||\mathbb{E}[g_{k}^{-}]||_{2}^{2}$
$\displaystyle=\frac{1}{c}\langle\mathbb{E}[g_{k}^{-}],\mathbb{E}[g_{k}^{-}]\rangle\leq\langle\nabla
f(x_{k})^{-},\mathbb{E}[g_{k}^{-}]\rangle$ $\displaystyle\leq
c\langle\mathbb{E}[g_{k}^{-}],\mathbb{E}[g_{k}^{-}]\rangle=c||\mathbb{E}[g_{k}^{-}]||_{2}^{2}.$

which exactly means

$\displaystyle\frac{1}{c}||\mathbb{E}[g_{k}]||_{2}^{2}$
$\displaystyle\leq\langle\nabla f(x_{k}),\mathbb{E}[g_{k}]\rangle\leq
c||\mathbb{E}[g_{k}]||_{2}^{2}.$

To prove the second inequality, we need to take advantage of the convexity of
$f$. Denote $[x_{k}-x^{*}]^{+}=\max\\{x_{k}-x^{*},\textbf{0}\\}$ and
$[x_{k}-x^{*}]^{-}=\min\\{x_{k}-x^{*},\textbf{0}\\}$, we know
$x_{k}-x^{*}=[x_{k}-x^{*}]^{+}+[x_{k}-x^{*}]^{-}$. In addition, because $f$ is
convex, the index set of non-zero entries of $[x_{k}-x^{*}]^{+}$ and $\nabla
f(x_{k})^{+}$ is the same. The index set of non-zero entries of
$[x_{k}-x^{*}]^{-}$ and $\nabla f(x_{k})^{-}$ is also the same. In addition,
because of Equation 11 and 12, the index set of non-zero entries of
$\mathbb{E}[g_{k}^{+}]$ ($\mathbb{E}[g_{k}^{-}]$) is the same with $\nabla
f(x_{k})^{+}$ ($\nabla f(x_{k})^{-}$). Combining these facts with Equations 11
and 12, we have

$\displaystyle\frac{1}{c}\langle\mathbb{E}[g_{k}^{+}],[x_{k}-x^{*}]^{+}\rangle$
$\displaystyle\leq\langle\nabla f(x_{k})^{+},[x_{k}-x^{*}]^{+}\rangle$
$\displaystyle\leq c\langle\mathbb{E}[g_{k}^{+}],[x_{k}-x^{*}]^{+}\rangle.$
$\displaystyle\frac{1}{c}\langle\mathbb{E}[g_{k}^{-}],[x_{k}-x^{*}]^{-}\rangle$
$\displaystyle\leq\langle\nabla f(x_{k})^{-},[x_{k}-x^{*}]^{-}\rangle$
$\displaystyle\leq c\langle\mathbb{E}[g_{k}^{-}],[x_{k}-x^{*}]^{-}\rangle.$

Combining these two equations, we have

$\displaystyle\frac{1}{c}\langle\mathbb{E}[g_{k}],x_{k}-x^{*}\rangle$
$\displaystyle\leq\langle\nabla f(x_{k}),x_{k}-x^{*}\rangle\leq
c\langle\mathbb{E}[g_{k}],x_{k}-x^{*}\rangle.$

This completes the proof. ∎

###  B.2 Proof of Theorem 4

######  Theorem 5.

(Main) Let $b,\epsilon,l,\delta,P,\alpha,\rho,\kappa$ and $\mathcal{B}_{l}$ be
as in section A in appendix, function
$f(x,\theta):\mathbb{R}^{d}\times\\{0,1\\}^{n}\rightarrow\mathbb{R}$ be a
$L$-smooth convex function w.r.t. $x$. Denote
$OPT=\min_{x}\mathbb{E}_{\theta\sim Pr(\theta)}f(x,\theta)$ as the global
optimum. Let $\sigma^{2}=\max_{x}\\{Var(\nabla_{x}f(x,\theta))\\}$ and
$\varepsilon^{2}=\max_{x}\\{||\mathbb{E}[\nabla_{x}f(x,\theta)]||_{2}^{2}\\}$.
For any $1\leq\rho\kappa\leq\sqrt{2}$, step size
$t\leq\frac{2-\rho^{2}\kappa^{2}}{L\rho\kappa}$ and sample size $N\geq 1$,
$\overline{x_{K}}$ is the output of XOR-SGD and
$\mbox{obj}=\mathbb{E}_{\theta}[f(\overline{x_{K}},\theta)]$ is the objective
function value at $\overline{x_{K}}$. We have:

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[\mbox{obj}]-OPT$
$\displaystyle\leq\frac{\rho\kappa||x_{0}-x^{*}||_{2}^{2}}{2tK}+\frac{t(\sigma^{2}+\varepsilon^{2})}{N}.$
(13)

###### Proof.

(Theorem 4) Since we use $N$ samples at each iteration, we have
$\overline{g_{k}}=\frac{1}{N}\sum_{i=1}^{N}g_{k}^{i}$ and
$\mathbb{E}[\overline{g_{k}}]=\mathbb{E}[g_{k}^{i}]$. In each iteration $k$ we
can adjust the parameters in XOR-Sampling to make the tail
$\epsilon\eta_{\phi}$ zero, then for each sample $g_{k}^{i}$ we can obtain
from Theorem 2 that

$\displaystyle\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}$
$\displaystyle\leq\mathbb{E}[g_{k}^{i+}]\leq\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}.$ (14) $\displaystyle\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}$
$\displaystyle\leq\mathbb{E}[g_{k}^{i-}]\leq\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}.$ (15)

The variance of each sample $g_{k}^{i}$ can also be bounded by

$\displaystyle Var(g_{k}^{i})$ $\displaystyle=\mathbb{E}_{\theta^{\prime}\sim
p^{\prime}(\theta^{\prime})}[||\nabla
f(x_{k},\theta^{\prime})||_{2}^{2}]-||\mathbb{E}_{\theta^{\prime}\sim
p^{\prime}(\theta^{\prime})}[\nabla f(x_{k},\theta^{\prime})]||_{2}^{2},$
$\displaystyle\leq\rho\kappa\mathbb{E}_{\theta\sim p(\theta)}[||\nabla
f(x_{k},\theta)||_{2}^{2}],$ $\displaystyle=\rho\kappa(Var(\nabla
f(x_{k},\theta))+||\mathbb{E}_{\theta\sim p(\theta)}[\nabla
f(x_{k},\theta)]||_{2}^{2}),$
$\displaystyle\leq\rho\kappa(\sigma^{2}+\varepsilon^{2}).$

Denote $\overline{g_{k}}^{+}=\max\\{\overline{g_{k}},\textbf{0}\\}$ and
$\overline{g_{k}}^{-}=\min\\{\overline{g_{k}},\textbf{0}\\}$. Clearly,
$g_{k}^{i+}\geq 0$ and $g_{k}^{i-}\leq 0$. Moreover, for a given dimension,
either $g_{k}^{i+}=0$ for that dimension or $g_{k}^{i-}=0$. Evaluating
$\overline{g_{k}}$ dimension by dimension, we can see that
$\overline{g_{k}}^{+}=\frac{1}{N}\sum_{i=1}^{N}g_{k}^{i+}$ and
$\overline{g_{k}}^{-}=\frac{1}{N}\sum_{i=1}^{N}g_{k}^{i-}$. Combined with
Equation 14 and 15, we know

$\displaystyle\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}$
$\displaystyle\leq\mathbb{E}[\overline{g_{k}}^{+}]\leq\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{+}.$ $\displaystyle\rho\kappa\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}$
$\displaystyle\leq\mathbb{E}[\overline{g_{k}}^{-}]\leq\frac{1}{\rho\kappa}\mathbb{E}_{\theta}[\nabla
f(x_{k},\theta)]^{-}.$

Because $\mathbb{E}[\overline{g_{k}}]=\mathbb{E}[g_{k}^{i}]$, we also have

$\displaystyle
Var(\overline{g_{k}})=\frac{1}{N^{2}}Var(\sum_{i=1}^{N}g_{k}^{i})=\frac{Var(g_{k}^{i})}{N}.$

Then the variance of $\overline{g_{k}}$ can be bounded as

$\displaystyle Var(\overline{g_{k}})$
$\displaystyle\leq\frac{\rho\kappa(\sigma^{2}+\varepsilon^{2})}{N}.$

Therefore, we can then apply Theorem 3 to get the result in equation 5.

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[\mathbb{E}_{\theta}[f(\overline{x_{K}},\theta)]]-\mathbb{E}_{\theta}[f(x^{*},\theta)]$
$\displaystyle\leq\frac{\rho\kappa||x_{0}-x^{*}||_{2}^{2}}{2tK}+\frac{t\max_{k}\\{Var(\overline{g_{k}})\\}}{\rho\kappa},$
$\displaystyle\leq\frac{\rho\kappa||x_{0}-x^{*}||_{2}^{2}}{2tK}+\frac{t(\sigma^{2}+\varepsilon^{2})}{N}.$

which can also be written as

$\displaystyle\mathbb{E}_{\overline{x_{K}}}[\mbox{obj}]-OPT$
$\displaystyle\leq\frac{\rho\kappa||x_{0}-x^{*}||_{2}^{2}}{2tK}+\frac{t(\sigma^{2}+\varepsilon^{2})}{N}.$
(16)

This completes the proof. ∎

##  Appendix C Experiments

We evaluate our XOR-SGD algorithm on the inventory management Ziukov [2016],
Shapiro and Philpott [2007] and the network design problems Sheldon et al.
[2012], Wu et al. [2017, 2016]. For each setting of both applications, to
produce a sample, Gibbs sampling first takes 100 steps to burn in, and then
draws samples every 30 steps. We fix the iteration step of both BP and BPChain
as $20$, which is enough for BP to converge. We allow SGD with Gibbs sampling,
BP and BPChain to draw more samples than XOR-SGD for a fair comparison. All
experiments were conducted using single core architectures on Intel Xeon Gold
6126 2.60GHz machines with 96GB RAM and a wall-time limit of $10$ hours. For
both applications, we use MRF as probabilistic models for $Pr(\theta)$, which
can be seen in the next section. For a fair comparison, once a solution $x$ is
generated by either algorithm, we use an exact weighted counter ACE Barton et
al. [2016] to evaluate $\mathbb{E}_{\theta\sim Pr(\theta)}f(x,\theta)$
exactly. All objective values reported here are from ACE.

###  C.1 Settings of Stochastic Inventory Management

Taking into account of the storage constraint, the original problem is
equivalent to the following problem:

$\displaystyle\min_{x\geq 0}\max_{\mu\geq 0}\mathbb{E}_{d\sim
Pr(d)}[G(x,d)]+\mu(w^{T}x-X).$ (17)

For inventory management problem, we assume each $d_{i}$ can take two
different values, one corresponding to the high demand one corresponding to
the low demand. Then, we introduce a new vector $\theta$ where $\theta_{i}=1$
means $d_{i}$ is the high value while $\theta_{i}=0$ otherwise. In the
experiment we range $n$ from 10 to 100 increased by a step size of 10 and draw
10 instances for each setting. Under each setting, we draw every $c_{i}$
uniformly from $(0,5]$, $h_{i}$ uniformly from $(0,10]$, sample $s_{i}$
uniformly drawn from $(0,10]$ and let $b_{i}=c_{i}+s_{i}$. The two values of
each $d_{i}$ are also uniformly drawn from $(0,10]$. We model $Pr(\theta)$ as
a MRF with several cliques. The variables in each clique are highly correlated
with each other. For a problem with $n$ products, we draw the number of
cliques uniformly from $[n,2n]$. The domain size of each clique
$\phi_{\alpha}$ is chosen from the range of $[1,6]$ at random. The potential
function of a clique involving $l$ variables is in the form of a table of size
$2^{l}$. The $i$-th entry of this table, denoted as $v_{i}$, is modeled as
$v_{i}=v_{i1}+v_{i2}v_{i3}$, where $v_{i1}$ is uniformly drawn from $(0,1)$,
$v_{3}$ uniformly from $(10,1000)$ and binary variable $v_{i2}$ uniformly
randomly drawn from $\\{0,1\\}$. Each storage requirement $w_{i}$ is drawn
from $(0,10]$ uniformly at random. The largest storage limit $X$ is set to be
$5n$. We also evaluate our method given different percentages of the largest
storage limit, which is shown in Figure 2 (middle). In the SGD algorithm, $x$
is initialized with the absolute value of a Gaussian random variable from
$\mathcal{N}(5,3)$ to ensure it is non-negative.

###  C.2 Settings of Stochastic Network Design

The task in equation 8 is equivalent to solving the following problem:

$\displaystyle\min_{\Delta g\geq 0}\max_{\mu\geq 0}\leavevmode\nobreak\
\mathbb{E}_{\theta\sim Pr(\theta)}[\overline{C}(g+\Delta
g,\theta)]+\mu(\sum_{e\in E}c_{e}\Delta g_{e}-B).$ (18)

Because of the convexity of $\overline{C}(g+\Delta g,\theta)$ and strong
duality, both problems have the same optimal solution.

We test our algorithm on a real-world problem, the so-called Flood Preparation
problem for the emergency medical services (EMS) on road networks Wu et al.
[2016]. The problem setup, including the graph structure and the definition of
$Pr(\theta)$, are the same as that in Wu et al. [2016]. The original network
is unweighted, hence we set the initial conductance value for each edge as 1.
$c_{e}$ is initialized uniformly from the range $(0,10)$. The largest budget
size $B$ is $1000$. We evaluate our method varying the percentage allowed of
the largest budget size, which is shown in Figure 3 (middle). In the
experiment, each entry of $\Delta g$ is initialized with the absolute value of
a Gaussian random variable from $\mathcal{N}(0,1)$. Total number of SGD
iterations is 2000, while not all algorithms can complete all 2000 iterations
within the time limit of 10 hours. The experimental results reported in the
plots are based on the best solutions found by each algorithm within the time
limit. Learning rate $t$ is $1$ at first and divided by $10$ after 20
iterations, further by $10$ after 100 iterations. Parameters in XOR-Sampling
are set to be the same as in the inventory management problem.