Gradient-Free Multi-Agent Nonconvex Nonsmooth Optimization

TL;DR

This paper introduces a distributed gradient-free algorithm for multi-agent nonconvex, nonsmooth optimization, demonstrating convergence to stationary points with proven sublinear rates, suitable for scenarios with limited information access.

Contribution

It presents a novel primal-dual gradient-free method for distributed nonconvex nonsmooth optimization with theoretical convergence guarantees.

Findings

Algorithm converges to stationary solutions

Achieves global sublinear convergence rate

Effective in numerical experiments for complex problems

Abstract

In this paper, we consider the problem of minimizing the sum of nonconvex and possibly nonsmooth functions over a connected multi-agent network, where the agents have partial knowledge about the global cost function and can only access the zeroth-order information (i.e., the functional values) of their local cost functions. We propose and analyze a distributed primal-dual gradient-free algorithm for this challenging problem. We show that by appropriately choosing the parameters, the proposed algorithm converges to the set of first order stationary solutions with a provable global sublinear convergence rate. Numerical experiments demonstrate the effectiveness of our proposed method for optimizing nonconvex and nonsmooth problems over a network.

Equations

$$\underset{x\in {\mathbb{R}}^M}{\min }f(x):=\sum _{i=1}^N{f}_i(x).$$

$${\mathbb{E}}_{\xi }[{\mathcal{H}}_i(\hat{x},\xi )]=f_i(\hat{x}),$$

$$\underset{\{x_i\}}{\min }$$

$$s.t.$$

$$L^{-}$$

$$L^{+}$$

$$\underset{x\in {\mathbb{R}}^Q}{\min }f(x),s.t.Ax=0.$$

$$f_{\mu }(x)$$

$$\nabla f_{\mu }(x)=\frac{1}{(2\pi {)}^{\frac{Q}{2}}}\int \frac{f(x+\mu \varphi )-f(x)}{\mu }\varphi e^{-\frac{1}{2}\Vert \varphi {\Vert }^2}d\varphi .$$

$$\Vert \nabla f_{\mu }(x)-\nabla f_{\mu }(y)\Vert \le L_1\Vert x-y\Vert .$$

$$G_{\mu }(x,\varphi ,\xi )=\frac{\mathcal{H}(x+\mu \varphi ,\xi )-\mathcal{H}(x,\xi )}{\mu }\varphi ,$$

$${\mathbb{E}}_{\xi ,\varphi }[G_{\mu }(x,\varphi ,\xi )]$$

$$U_{\rho }(x,\lambda ):=f(x)+⟨\lambda ,Ax⟩+\frac{\rho }{2}{\Vert Ax\Vert }_2^2.$$

$${\bar{G}}_{\mu ,i}(x_i^r,{\zeta }_i^r)=\frac{1}{J}\sum _{j=1}^J\frac{{\mathcal{H}}_i(x_i^r+\mu {\varphi }_{i,j}^r,{\xi }_{i,j}^r)-{\mathcal{H}}_i(x_i^r,{\xi }_{i,j}^r)}{\mu }{\varphi }_{i,j}^r.$$

$${\mathbb{E}}_{\zeta }[\Vert G_{\mu }^{J,r}{\Vert }^2]\le \frac{L_0^2(Q+4{)}^2}{J^2}.$$

$$x^{r+1}$$

$${\lambda }^{r+1}$$

$$G_{\mu }^{J,r}+A^{\top }({\lambda }^r+\rho A{x}^r)+2\rho D(x^{r+1}-x^r)=0.$$

$$\displaystyle x^{r+1}=\frac{1}{2\rho}D^{-1}\bigg{[}\rho L^{+}x^{r}-G_{\mu}^{J,r}+A^{\top}\lambda^{r}\bigg{]}.$$

$${\lambda }_{ij}^{r+1}={\lambda }_{ij}^r+\rho (x_i^{r+1}-x_j^{r+1}),$$

$$⟨b-a,b⟩=\frac{1}{2}(\Vert b{\Vert }^2+\Vert a-b{\Vert }^2-\Vert a{\Vert }^2),$$

$$⟨a,b⟩\le \frac{1}{2ϵ}\Vert a{\Vert }^2+\frac{ϵ}{2}\Vert b{\Vert }^2;\forall ϵ>0.$$

$$\displaystyle\bigg{\|}\sum_{i=1}^{n}a_{i}\bigg{\|}^{2}\leq n\sum_{i=1}^{n}\big{\|}a_{i}\big{\|}^{2}.$$

$${\mathbb{E}}_{\zeta }\Vert G_{\mu }^{J,r}-\nabla f_{\mu }(x^r){\Vert }^2$$

$$\le {\mathbb{E}}_{\zeta }\Vert G_{\mu }^{J,r}{\Vert }^2$$

$$\le \frac{L_0^2(Q+4{)}^2}{J^2},$$

$$\frac{1}{\rho }\mathbb{E}\Vert {\lambda }^{r+1}-{\lambda }^r{\Vert }^2$$

$$V^{r+1}:=\frac{\rho}{2}\big{(}\|Ax^{r+1}\|^{2}+\|x^{r+1}-x^{r}\|^{2}_{B}\big{)},$$

$$P^{r+1}:=U_{\rho ,\mu }^{r+1}+c{V}^{r+1}.$$

$$\displaystyle\mathbb{E}\big{[}P^{r+1}-P^{r}\big{]}$$

$${\alpha }_1:={\rho }^2-(2c{L}_1+L_1/2+L_1^2/2+1/2)\rho -\frac{6L_1^2}{{\sigma }_{\min }},$$

$${\alpha }_2:=\frac{3\rho \Vert L^{+}\Vert }{{\sigma }_{\min }}-\frac{c\rho }{2},{\alpha }_3=\frac{9}{\rho {\sigma }_{\min }}+\frac{6c+1}{L_1}.$$

$$c>\frac{6\Vert L^{+}\Vert }{{\sigma }_{\min }},\rho >b+\sqrt{b^2+6L_1^2/{\sigma }_{\min }},$$

$$\mathbb{E}[P^{r+1}]\ge \underline{P}>-\infty ,\forall r\ge 1.$$

Full text

1 Introduction

Consider a network with $N$ distributed agents that collectively solve the following optimization problem

[TABLE]

Here, $f_{i}:\mathbb{R}^{M}\to\mathbb{R}$ is possibly a nonconvex nonsmooth function, that is only available to agent $i$. Such distributed optimization problems arise in many applications such as machine learning [1, 2], resource allocation [3], robotic networks [4], and signal processing [5]. See [6] for more applications.

Many distributed optimization methods have been proposed to solve problem (1). Many of them rely on consensus between the agents and assume that the cost functions are convex. One of the first such methods is the distributed subgradient (DSG) algorithm [7]. Subsequently, a number of similar consensus-based algorithms were also proposed to solve distributed convex optimization problems in the form of (1); see, e.g., [8, 9, 10, 11]. These methods only converge to a neighborhood of the solution set unless they use diminishing stepsizes, which often makes them slow. Faster algorithms using constant stepsizes include the incremental aggregated gradient (IAG) method [12], the exact first-order algorithm (EXTRA) [13], and Accelerated Distributed Augmented Lagrangian (ADAL) algorithm. See also [14, 15, 16] for optimization of convex problems.

Optimization of nonconvex functions is a much more challenging problem. Only recently, have there been developed a few nonconvex distributed optimization algorithms motivated by applications in resource allocation in ad-hoc network [17], sparse PCA [18], and flow control in communication networks [19]; see also [20, 21, 22, 19, 23, 24, 25, 26, 27] for additional algorithms developed for nonconvex optimization problems.

Regardless of convexity and/or smoothness, all the aforementioned methods require that either the first order (gradient/subgradient) information or the explicit form of the objective function is available to the agents. However, in many important practical situations such information can be expensive to obtain, or even impossible. Examples include, simulation-based optimization where the objective function can only be evaluated using repeated simulation [28], training deep neural networks where the relationship between the variables and the cost function is too complicated to derive an explicit form of the gradient [29], and bandit optimization where a player optimizes a sequence of cost functions having only knowledge of a single function value each time [30]. In these cases, the zeroth-order information about the objective function values is often readily available. Such zeroth-order information can be obtained through a stochastic zeroth-order oracle ($\mathcal{SZO}$). In particular, suppose that $\hat{x}\in\mbox{dom}(f_{i})$, then the $i$th $\mathcal{SZO}$ at agent $i$ returns a noisy version of $f_{i}(\hat{x})$ denoted by $\mathcal{H}_{i}(\hat{x},\xi)$ that satisfies

[TABLE]

where $\xi$ is a random variable representing the noise.

Recently, centralized zeroth-order optimization has received significant attention. In [31], Nesterov proposed a general framework for analyzing zeroth-order algorithms and provided the global convergence rates for both convex and nonconvex problems. In [32], the authors established a stochastic zeroth-order method, which again can deal with both convex and nonconvex (but smooth) optimization problems. In [33] a Mirror Descent based zeroth-oder algorithm was proposed for solving convex optimization problems. Both these zeroth-order algorithms are centralized and cannot be implemented over a multi-agent network. A few recent works [34, 35, 36] considered zeroth-order distributed convex (possibly nonsmooth) problems, but none of these works can address nonconvex problems. For distributed nonconvex but smooth optimization problem using zeroth-order information a primal-dual algorithm has been proposed in [37, 38].

In this paper, we propose a new algorithm for distributed nonconvex and nonsmooth optimization with zeroth-order information. Specifically, we first show that this problem can be reformulated as a linearly constrained optimization problem over a connected multi-agent network. Then, we propose a nonconvex primal-dual based algorithm, which requires only local communication among the agents, and utilizes local zeroth-order information. Theoretically, We show that the current solution converges approximately to a stationary solution of the problem. We also provide numerical results that corroborate the theoretical findings.

Notation. Given a vector $a$ and a matrix $A$, we use $\|a\|$ and $\|A\|$ to denote the Euclidean norm of vector $a$, and spectral norm of matrix $A$, respectively. $A^{\top}$ represents the transpose of matrix $A$. We define $\|a\|^{2}_{A}:=a^{T}Aa$. The notation $\langle a,b\rangle$ is used to denote the inner product of two vectors $a$, $b$. For matrices $A$ and $B$, $A\otimes B$ is the Kronecker product of $A$ and $B$. To denote an $M\times M$ identity matrix we use $I_{M}$. $\mathbb{E}[\cdot]$ denotes the expectation with respect to all random variables, and $\mathbb{E}_{v}[\cdot]$ denotes the expectation with respect to the random variable $v$.

2 Problem Definition and Proposed Algorithm

Let us define a graph $\mathcal{G}=\{\mathcal{V},\mathcal{E}\}$, where $\mathcal{V}$ is the node set with $|\mathcal{V}|=N$, and $\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$ is the edge set with $|\mathcal{E}|=E$. We assume that $\mathcal{G}$ is undirected, meaning that if $(i,j)\in\mathcal{E}$ then $(j,i)\in\mathcal{E}$. Moreover, every agent $i$ can only communicate with its direct neighbors in the set $\mathcal{N}_{i}=\{j\in\mathcal{V}\;|\;(i,j)\in\mathcal{E}\}$, and let $d_{i}=|\mathcal{N}_{i}|$ denote the degree of node $i$. We assume that the graph $\mathcal{G}$ is connected, meaning that there is a path, i.e., a sequence of nodes where consecutive nodes are neighbors, between any two nodes in $\mathcal{G}$.

In order to decompose problem (1) let us introduce $N$ new variables $x_{i}\in\mathbb{R}^{M}$ that are local to every agent $i$. Then, problem (1) can be reformulated as follows:

[TABLE]

The set of constraints $x_{i}=x_{j}$ enforce consensus on the local variables $x_{i}$ and $x_{j}$ for all neighbors $j\in\mathcal{N}_{i}$. We stack all the local variables $x_{i}$ in a vector $x:=\{x_{i}\}\in\mathbb{R}^{Q\times 1}$, where $Q=NM$. Moreover, we define the Degree matrix $\tilde{D}\in\mathbb{R}^{N\times N}$ to be a diagonal matrix where $\tilde{D}(i,i)=d_{i}$; let $D=\tilde{D}\otimes I_{M}\in\mathbb{R}^{Q\times Q}$. For a given graph $\mathcal{G}$, the incidence matrix $\tilde{A}\in\mathbb{R}^{E\times N}$ is a matrix where $\tilde{A}(k,i)=1$ and $\tilde{A}(k,j)=-1$, where $(i,j)$ is the $k$th edge of $\mathcal{G}$; the rest of the entries of $\tilde{A}$ are all zero. Let $A=\tilde{A}\otimes I_{M}\in\mathbb{R}^{EM\times Q}$. Finally, we define the Signed and Signless Laplacian matrix, denoted by $L^{-}$ and $L^{+}$, respectively as

[TABLE]

Using the above notations, problem (3) can be written in the following compact form:

[TABLE]

2.1 Preliminaries

In this section, we first introduce some standard techniques presented in [31] for approximating and smoothing the gradient of a given function. Suppose that $\phi\in\mathbb{R}^{Q}$ is a Gaussian random vector and let $\mu>0$ be some smoothing parameter. The smoothed version of function $f$ is defined as

[TABLE]

Then it can be shown that the function $f_{\mu}$ is differentiable and its gradient is given by Eq. (22) in [31]

[TABLE]

Further, assuming that the original function $f$ is Lipschitz continuous, i.e., there exists $L_{0}$ such that $|f(x)-f(y)|\leq L_{0}|x-y|$ for all $x,y\in\mbox{dom\,}(f)$, it can be shown that (see [31, Lemma 2]) $\nabla f_{\mu}$ is also Lipschitz continuous with constant $L_{1}=\frac{2L_{0}\sqrt{Q}}{\mu}$. In other words, for all $x,y\in\mbox{dom\,}(f_{\mu})$ we have

[TABLE]

Let $\mathcal{H}(x,\xi)$ denote the noisy functional value of the function $f$ obtained from an associated $\mathcal{SZO}$ as in equation (2). In view of (8), the gradient of $f_{\mu}(x)$ can be approximated as

[TABLE]

where the constant $\mu>0$ is the smoothing parameter. It can be easily checked that $G_{\mu}(x,\phi,\xi)$ is an unbiased estimator of $\nabla f_{\mu}(x)$, i.e.,

[TABLE]

For simplicity we define $\zeta:=(\xi,\phi)$. For a given number $J$ of independent samples of $\{\zeta_{j}\}_{j=1}^{J}$, we define the sample average $\bar{G}_{\mu}(x,\zeta):=\frac{1}{J}\sum_{j=1}^{J}G_{\mu}(x,\zeta_{j})$, where $\zeta:=\{\zeta_{j}\}_{j=1}^{J}$. It is easy to see that for any $J\geq 1$, $\bar{G}_{\mu}(x,\zeta)$ is also an unbiased estimator of $\nabla f_{\mu}(x)$.

2.2 The Proposed Algorithm

In this part we propose a primal-dual algorithm for the distributed optimization problem (6). Let $\lambda_{ij}\in\mathbb{R}^{M}$ be the multiplier associated with the consensus constraint $x_{i}-x_{j}=0$ for each $(i,j)\in\mathcal{E}$. Moreover, stack all $\lambda_{ij}$’s in a vector $\lambda=\{\lambda_{ij}\}_{(i,j)\in\mathcal{E}}\in\mathbb{R}^{EM}$. Then, the augmented Lagrangian (AL) function for problem (6) is given by

[TABLE]

where $\rho>0$ is a constant. Moreover, as in (7) define the smoothed version $f_{i,\mu}$ of the local function $f_{i}$. At iteration $r$ of the algorithm we obtain an unbiased estimation of the gradient of local function $f_{i,\mu}(x_{i}^{r})$ as follows. For every sample $j\in\{1,2,\cdots,J\}$ we generate a random vector $\phi^{r}_{i,j}\in\mathbb{R}^{M}$ from an i.i.d standard Gaussian distribution and calculate $\bar{G}_{\mu,i}(x_{i}^{r},\zeta^{r}_{i})\in\mathbb{R}^{M}$ similar to (10) by

[TABLE]

Define ${G^{J,r}_{\mu}}:=\{\bar{G}_{\mu,i}(z_{i}^{r},\zeta^{r}_{i})\}_{i=1}^{N}\in\mathbb{R}^{Q}$. The following theorem bounds the norm of $G^{J,r}_{\mu}$.

Theorem 1** **(Theorem 4 [31])

If $f$ is a Lipschitz continuous function with constant $L_{0}$, then

[TABLE]

Our proposed algorithm is summarized in Algorithm 1. In the primal step (15), an approximate gradient descent step is taken towards minimizing the augmented Lagrangian function with respect to $x$. In particular, in the first-order approximation of the AL function the true gradient of the function $f(x^{r})$ is approximated by the noisy zeroth-order estimate $G_{\mu}^{J,r}$ and then, a matrix-weighted quadratic penalty $\rho\|x-x^{r}\|_{D}$ is used. This term is critical for the algorithm itself, as well as for the analysis. The dual step (16) is then performed, which is a gradient ascent step over the dual variable $\lambda$.

To see how Algorithm 1 can be implemented in a distributed way, consider the optimality condition for (15) as

[TABLE]

Utilizing (4) and (5), we have

[TABLE]

To implement this primal iteration, each agent $i$ only requires local information as well as information from its neighbors $\mathcal{N}_{i}$. This is because $D$ is a diagonal matrix and the structure of the matrix $L^{+}$ ensures that the $i$th block vector of $L^{+}x^{r}$ is only related to $x^{r}_{j},\;j\in\mathcal{N}_{i}$. For the dual step w.l.o.g we assign the dual variable $\lambda_{ij}$ to node $i$ and therefore, from (16) we have

[TABLE]

which only requires the local information as well as information from the neighbors in $\mathcal{N}_{i}$.

3 The Convergence Analysis

In this section we study the convergence of Algorithm 1. We make the following assumptions.

Assumptions A. We assume that

• A1.

The function $f$ is Lipschitz continuous.

• A2.

The function $f$ is lower bounded.

The above assumptions on the objective $f$ are quite standard in the analysis of first order optimization Algorithms (1). To simplify notation let $\mathcal{F}^{r}:=\sigma(\zeta_{1},\zeta_{2},\cdots\zeta_{r})$ be the $\sigma$-field generated by the entire history of algorithm up to iteration $r$, $\sigma_{\min}$ be the smallest nonzero eigenvalue of $A^{\top}A$, and $w^{r}:=(x^{r+1}-x^{r})-(x^{r}-x^{r-1})$ be the successive difference of the differences of the primal iterates. In the analysis that follows we will make use of the following relations:

• •

For any given vectors $a$ and $b$ we have

[TABLE]

• •

For $n$ given vectors $a_{i}$ we have that

[TABLE]

Our proof consists of a series of lemmas leading to the main convergence rate result. In our presentation we try to provide insights about the analysis steps, while the proofs of the results are relegated to the appendix.

First we bound the difference between $\nabla f_{\mu}(x^{r})$ and its unbiased estimation $G^{J,r}_{\mu}$ as follows:

[TABLE]

where the last inequality follows from (14). The next lemma bounds the change of the dual variables by that of the primal variables. The proofs of the results that follow can be found in the appendix.

Lemma 1

Suppose Assumptions A hold true. Let $L_{1}$ denote the gradient Lipschitz constant for function $f_{\mu}$. Then we have the following inequity:

[TABLE]

The next step is the key in our analysis. We define the smoothed version of the AL function in a similar way as (7) and denote it by $U_{\rho,\mu}(x,\lambda)$. For notational simplicity let us define $U^{r+1}_{\rho,\mu}:=U_{\rho,\mu}(x^{r+1},\lambda^{r+1})$. From equation (9) we know that function $f_{\mu}$ is Lipschitz continuous with constant $L_{1}$. Now let $c>0$ be some positive constant and set k:=2\big{(}\frac{6{L_{1}}^{2}}{\rho\sigma_{\min}}+\frac{3c{L_{1}}}{2}\big{)}. Moreover, we define

[TABLE]

where $B:=L^{+}+\frac{k}{c\rho}I_{Q}$. Finally we define the following potential function:

[TABLE]

We study the behavior of the proposed potential function as the algorithm proceeds.

Lemma 2

Suppose Assumptions A hold true. We have that

[TABLE]

where

[TABLE]

Note that the constants $\alpha_{1}$ and $\alpha_{2}$ in (25) can be made positive as long as we choose constants $c$ and $\rho$ large enough. In particular, the following conditions are sufficient to have positive $\alpha_{1},\alpha_{2}$

[TABLE]

where $b=(cL_{1}+L_{1}/4+L_{1}^{2}/4+1/4)$.

The key insight obtained from Lemma 2 is that, a proper combination of a primal objective (i.e., the AL function), and the dual gap (i.e., the violation of the feasibility) can be served as the potential function that guides the progress of the algorithm.

In the next lemma we show that $P^{r+1}$ is lower bounded.

Lemma 3

Suppose Assumptions A hold true, and constant $c$ is selected as $c\geq\frac{2\|L^{+}\|}{\sigma_{\min}}.$ Then there exists a constant $\underline{P}$ that is independent of the total number of iterations $T$ so that

[TABLE]

To characterize the convergence rate of Algorithm 1, let us define the stationarity gap of the smoothed version of problem (6) as

[TABLE]

It can be easily checked that $\Phi_{\mu}(x^{*},\lambda^{*})=0$ if and only if $(x^{*},\lambda^{*})$ is a KKT point of the smoothed version of problem (6). For simplicity let us denote $\Phi_{\mu}^{r}:=\Phi_{\mu}(x^{r},\lambda^{r-1})$.

At this point we are ready to combine the previous results to obtain our main theorem.

Theorem 2

Suppose Assumptions A hold true, the penalty parameter $\rho$ satisfies the condition given in Lemma 2, and the constant $c$ satisfies $c\geq\frac{6\|L^{+}\|}{\sigma_{\min}}.$ Then, there exist constants $\gamma_{1},\gamma_{2}>0$ such that

[TABLE]

From Theorem 2 we can observe that there exists always a constant term in the right-hand-side of the stationarity gap. Therefore, no matter how many iterations we run the algorithm, we always converge to a neighborhood of a stationary point. However, if we choose the number of samples $J\in\mathcal{O}({\sqrt{T}})$, we have the following bound:

[TABLE]

which verifies the sublinear convergence rate for the algorithm.

4 Numerical Results

In this section we illustrate the proposed algorithm through numerical simulations. For our experiments we study two nonconvex distributed optimization problems.

First, we consider a simple nonconvex nonsmooth distributed optimization problem defined bellow:

[TABLE]

where for each agent we have

[TABLE]

In (32) the problem dimension $M=1$ and the number of nodes in the network is $N=10$. Therefore, $Q=MN=10$. The details of the underlying graph are discussed in [39]. We compare Algorithm 1 using a constant stepsize that satisfies (27) and the Randomize Gradient Free (RGF) algorithm proposed in [34] using diminishing stepsize $\frac{1}{\sqrt{r}}$. Note that in theory RGF only works for the convex problems. However, we include it here for the purpose of comparison only. We compare the two algorithms in terms of the stationarity gap defined in (29) and the constraint violation $\|Ax\|$. The stopping criterion is set to $T=1000$ iterations and the results are the average over $30$ independent trials. Figs 1(a) and 1(b) show our comparative results. We observe that the stationarity gap and the consensus error vanish faster for our proposed algorithm than for the RGF algorithm.

Next, we study a mini-batch binary classification problem using nonconvex nonsmooth regularizers, where each node stores $b$ (batch size) data points. For this problem the local function is given by

[TABLE]

where $v_{ij}\in\mathbb{R}^{M}$ and $y_{ij}\in\{1,-1\}$ are the feature vector and the label for the $j$th data point of $ith$ agent [40]. The nonconvex nonsmooth regularization term $\log(\epsilon+\|x_{i}\|_{1})$ imposes sparsity to vector $x_{i}$, the constant $\alpha$ controls the sparsity level, and $\epsilon>0$ is a small number. The network has $N=15$ nodes and each node contains randomly generated $b=100$ data point. Algorithm 1 and RGF run for $T=1000$ iterations and Figures 2(a) and 2(b) illustrate the stationarity gap and the constraint violation versus the iteration counter. From these plots we can observe that, as in the previous problem, Algorithm 1 is faster than the RGF. Note again that, in theory, RGF is designed for convex problems only. To the best of our knowledge, our algorithm is the first provable distributed zeroth-order method for nonconvex and nonsmooth problems.

5 Conclusion

In this work, we proposed a distributed gradient-free optimization algorithm to solve nonconvex and nonsmooth problems utilizing local zeroth-order information. We rigorously analyzed the convergence rate of the proposed algorithm and demonstrated its performance via simulation. To the best of our knowledge, this is the first distributed framework for the solution of nonconvex and nonsmooth distributed optimization problems that also has a provable sublinear convergence rate. The proposed framework can be used to solve a variety of problems where access to first or second order information is very expensive or even impossible.

APPENDIX

5.1 Proof of Lemma 1

From equation (16) we have

[TABLE]

Equation (33) implies that $\lambda^{r+1}-\lambda^{r}$ lies in the column space of $A$, therefore we have

[TABLE]

where $\sigma_{\min}$ denotes the smallest non-zero eigenvalue of $A^{\top}A$. Utilizing equation (33) and equation (17), we obtain

[TABLE]

Replacing $r$ with $r-1$ in equation (35) and then using the definition of $w^{r}:=(x^{r+1}-x^{r})-(x^{r}-x^{r-1})$ we obtain

[TABLE]

where the last inequality follows from (21). Adding and subtracting $\nabla f_{\mu}(x^{r-1})$ to the second term on the r.h.s of (5.1) and taking the expectation on both sides gives

[TABLE]

where the last inequity is true because of (3), the fact that $\nabla f_{\mu}(z)$ is gradient Lipschitz with constant $L_{1}$, and the inequality $\|L^{+}w^{r}\|^{2}\leq\|L^{+}\|\|w^{r}\|^{2}_{L^{+}}$. The proof is complete.

5.2 Proof of Lemma 2

First we prove that the function $g(x):=U_{\rho,\mu}(x,\lambda)+\frac{\rho}{2}\|x-x^{r}\|^{2}_{L^{+}}$ is strongly convex with respect to variable $x$ when $2\rho\geq L_{1}$. From Assumption A.1 and the fact that $D\succeq I$ we have

[TABLE]

This proves that $U_{\rho,\mu}(x,\lambda)+\frac{\rho}{2}\|x-x^{r}\|^{2}_{L^{+}}$ is strongly convex with modulus $2\rho-L_{1}$. Using this fact, we can bound $U_{\rho,\mu}^{r+1}-U_{\rho,\mu}^{r}$ as follows:

[TABLE]

where the last inequity holds true due to the strong convexity of $U_{\rho,\mu}(x,\lambda)+\frac{\rho}{2}\|x-x^{r}\|^{2}_{L^{+}}$ and (33). Now using (35) we have

[TABLE]

here the last inequality follows from (20) with $\epsilon=L_{1}^{2}$. Taking expectation on both sides we have

[TABLE]

where the first inequality follows from (5.1) and the second inequality follows from (3). Now we bound $V^{r+1}-V^{r}$. From the optimality condition for problem (15) and the dual update (16) we have

[TABLE]

Similarly, for the $(r-1)$th iteration, we have

[TABLE]

Setting $x=x^{r}$ in first equation, $x=x^{r+1}$ in second equation, and adding them, we obtain

[TABLE]

The l.h.s can be expressed as follows:

[TABLE]

where the first equality follows from (16) and the second equality follows from (19). For the r.h.s of (40) we have

[TABLE]

where the first inequality follows from (20). To get the second inequality we add and subtract $\nabla f_{\mu}(x^{r})+\nabla f_{\mu}(x^{r-1})$ to $G_{\mu}^{J,r}-G_{\mu}^{J,r-1}$ and use (21). Taking expectation on both sides, we have

[TABLE]

where the inequality follows from (3) and the last equality follows from (19). Combining (5.2) and (5.2), we obtain

[TABLE]

Rearranging terms in (5.2), and using the definition of $V^{r}$ and $B$, we have

[TABLE]

Recall the definition of $P^{r+1}:=U_{\rho,\mu}^{r+1}+cV^{r+1}$. Utilizing (39) and (5.2) eventually we obtain

[TABLE]

where we have

[TABLE]

The lemma is proved.

5.3 Proof of Lemma 3

From (35) we have

[TABLE]

Recall that $\sigma_{\min}$ is the smallest nonzero eigenvalue of $A^{T}A$. Also from the fact that $\lambda^{0}=0$ we have that the dual variable lies in the column space of $A$, thus we get that

[TABLE]

From the definition of the potential function we have

[TABLE]

where $B:=L^{+}+\frac{k}{c\rho}I$. From Assumption A.2 we have that $f$ is lower bounded, therefore, $f_{\mu}$ is lower bounded too, i.e. there exists $\underline{f}$ such that $\underline{f}\leq f_{\mu}(x)$ for all $x\in\mbox{dom\,}(f_{\mu})$. Plugging (46) in (5.3), and utilizing the fact that $\|Ax^{r+1}+\frac{1}{\rho}\lambda^{r+1}\|^{2}\geq 0$ we obtain

[TABLE]

where the last inequality is due to the fact that $\frac{k}{2}\|x^{r+1}-x^{r}\|^{2}\geq 0$. Notice that $L^{+}$ is a symmetric PSD matrix. Therefore, picking constant $c$ large enough such that $c\geq\frac{2\|L^{+}\|}{\sigma_{\min}}$, we have $\frac{c\sigma_{\min}}{2}L^{+}-(L^{+})^{2}\succeq 0$. Hence, with this choice of $c$ we get the following bound for the potential function

[TABLE]

Taking expectation on both sides we have

[TABLE]

where the last inequality follows from (14). To complete the proof we only need to set $\underline{P}=-\frac{L_{0}(Q+4)^{2}}{\sigma_{\min}J^{2}}+\underline{f}$.

5.4 Proof of Theorem 2

First we bound the stationarity gap given in (29). We have

[TABLE]

where the first equality follows from (16), the second equality follows from (35), and the first inequality follows from (21). Taking expectation on both sides gives

[TABLE]

where the inequality follows from (3). Next, we bound the expected value of the constraint violation. Utilizing the equation (16) we have

[TABLE]

Taking expectation and utilizing (23) we reach

[TABLE]

Summing (50) and (5.4), we have the following bound for the stationarity gap

[TABLE]

where $\beta_{1},\beta_{2},\beta_{3},\beta_{4}$ are positive constants given by

[TABLE]

Summing both sides of (52) over $T$ iterations, we get

[TABLE]

Applying Lemma 2 and summing both sides of (25) over $T$ iterations, we obtain

[TABLE]

Let us set $\tau=\frac{\max(\beta_{1}+\beta_{2},\beta_{3})}{\min(\alpha_{1},\alpha_{2})}$. Combining the two inequalities (53) and (54) and utilizing the fact that $\mathbb{E}[P^{T+1}]$ is lower bounded by $\underline{P}$, we arrive at the following inequality

[TABLE]

Because $u$ is a uniform random number from $\{0,1,\cdots,T-1\}$ we have

[TABLE]

Combining (55) and (56) implies the following

[TABLE]

Setting

[TABLE]

we conclude the proof.