Improved Convergence Rates for Distributed Resource Allocation

Angelia Nedić, Alex Olshevsky, Wei Shi

I Introduction

This paper deals with a decentralized resource allocation problem, which is defined over a connected network of nn agents, as follows:

I-B Our contributions

II Resource Allocation and Its Connection to Consensus Optimization

Some of the notation may not be standard but it enables us to present our algorithm and analysis in a compact form. Throughout the paper, we let agent ii hold a local variable xix_{i}, a function fif_{i}, and a constraint set Ωi\Omega_{i} of problem (1). We define

Our basic assumption is that problem (1) is convex, which is formalized as follows.

We define gig_{i} as the indicator function of the set Ωi\Omega_{i}, namely,

We also define a composite function hih_{i} for agent ii, as follows:

The equality in (2) holds when ri(dom{fi})⋂ri(dom{gi})≠∅\text{ri}(\text{dom}\{f_{i}\})\bigcap\text{ri}(\text{dom}\{g_{i}\})\neq\emptyset, where ri(⋅)\text{ri}(\cdot) denotes the relative interior of a set and dom{⋅}\text{dom}\{\cdot\} is the (effective) domain of a function (see Section 4 of for the definition of “(effective) domain”). See also Remark 16.46 and Corollary 16.48 of for more conditions and comments for (2) to hold. Note that we have not imposed any differentiability on fif_{i}’s. Moreover, since ∂gi(xi)\partial g_{i}(x_{i}) coincides with the normal cone of Ωi\Omega_{i} at xi∈Ωix_{i}\in\Omega_{i}, we have that

In addition to the network objective f(x)\mathbf{f}(\mathbf{x}) defined in (1a), we introduce two more network-wide aggregate functions,

Similarly, we define a matrix r\mathbf{r} by using the vectors rir_{i}, i=1,…,ni=1,\ldots,n.

Letting ∇~fi(xi)\widetilde{\nabla}f_{i}(x_{i}) be a subgradient of fif_{i} at xix_{i}, we construct a matrix ∇~f(x)\widetilde{\nabla}\mathbf{f}(\mathbf{x}) of subgradients ∇~fi(xi)\widetilde{\nabla}f_{i}(x_{i}), as follows:

and, similarly, the matrices ∇~g(x)\widetilde{\nabla}\mathbf{g}(\mathbf{x}) and ∇~h(x)\widetilde{\nabla}\mathbf{h}(\mathbf{x}) are defined using subgradients of gig_{i} and hi=fi+gih_{i}=f_{i}+g_{i} at xix_{i}, respectively. We drop the tilde in the notation ∇~\widetilde{\nabla} when the function under consideration is differentiable (i.e., a subdifferential set contains only a gradient). Each row ii of x\mathbf{x}, r\mathbf{r}, ∇~f(x)\widetilde{\nabla}\mathbf{f}(\mathbf{x}), ∇~g(x)\widetilde{\nabla}\mathbf{g}(\mathbf{x}), and ∇~h(x)\widetilde{\nabla}\mathbf{h}(\mathbf{x}) corresponds to the information available to agent ii only.

To model the underlying communication network for the agents, we use a simple (no self-loop) undirected graph, where [n]={1,2,…,n}[n]=\{1,2,\ldots,n\} is the vertex set and E\mathcal{E} is the edge set. We say that an n×nn\times n matrix AA is compatible with the graph G\mathcal{G} when the following property holds: ∀i,j∈[n]\forall i,j\in[n], the (i,j)(i,j)-th entry of AA is zero if neither {i,j}\{i,j\} is an element of E\mathcal{E} nor i≠ji\neq j. We use Ni{\cal N}_{i} to denote the set of neighbors of agent ii in the graph G\mathcal{G}, i.e., Ni={j∈[n]∣{i,j}∈E}.{\cal N}_{i}=\{j\in[n]\mid\{i,j\}\in\mathcal{E}\}.

Let \LG{\textbf{\L}}_{\mathcal{G}} denote the (standard) Laplacian matrix associated with the graph G\mathcal{G}, i.e., \LG=D−J{\textbf{\L}}_{\mathcal{G}}=D-J, where DD is the diagonal matrix with diagonal entries Dii=diD_{ii}=d_{i} and did_{i} being the number of edges incident to node ii, while JJ is the graph adjacency matrix (with Jij=1J_{ij}=1 when {i,j}∈E\{i,j\}\in\mathcal{E} and Jij=0J_{ij}=0 otherwise). A few facts about \LG{\textbf{\L}}_{\mathcal{G}} are that \LG{\textbf{\L}}_{\mathcal{G}} is compatible with G\mathcal{G}, symmetric and positive semidefinite.

In our algorithm, we will use a matrix \L=[\Lij]\textbf{{\L}}=[\text{\L}_{ij}] whose behavior is “close or the same” to \LG{\textbf{\L}}_{\mathcal{G}}, in the sense of the following assumption.

Since the graph Laplacian \LG{\textbf{\L}}_{\mathcal{G}} satisfies Assumption 2, we can choose \L=\LG\textbf{{\L}}={\textbf{\L}}_{\mathcal{G}}. In this case, each agent needs to know the number of its neighbors (its degree) and Ł can be constructed without any communication among the agents.

We can let \L=\LG/λmax⁡{\LG}\textbf{{\L}}={\textbf{\L}}_{\mathcal{G}}/\lambda_{\max}\{{\textbf{\L}}_{\mathcal{G}}\}. The network needs λmax⁡{\LG}\lambda_{\max}\{{\textbf{\L}}_{\mathcal{G}}\} to configure this matrix but a preprocessing to retrieve λmax⁡{\LG}\lambda_{\max}\{{\textbf{\L}}_{\mathcal{G}}\} is possible .

We can also choose \L=0.5(I−W)\textbf{{\L}}=0.5(I-W) where WW is a symmetric doubly stochastic matrix that is compatible with the graph G\mathcal{G} and λmax⁡{W−11⊤/n}\lambda_{\max}\{W-\mathbf{1}\mathbf{1}^{\top}/n\} is strictly less than 11. This matrix can be constructed in the network through a few rounds of local interactions between the agents since some local strategies for determining WW exist, such as the Metropolis-Hasting rule which requires only one round of local interactions .

We will discuss the specific choices of Ł in some of our results to simplify analysis or to point out to interesting results.

II-B The resource allocation and consensus optimization problems

In this subsection, we investigate the first-order optimality conditions for problem (1) and for consensus optimization. With the notation introduced in the preceding section, the resource allocation problem (1) can be compactly given by

where 1\mathbf{1} is a vector of appropriate dimension whose entries are all equal to 11. By using the Lagrangian function, we can write down the optimality conditions for problem (6) in a special form, as given in the following lemma.

where UU is the matrix defined in Assumption 2.

The Lagrangian function of problem (6) is

where ∂h(x∗)\partial\mathbf{h}(\mathbf{x}^{*}) is to be understood as a collection of all matrices whose every row ii is given by some subgradient (∇~hi(xi∗))⊤(\widetilde{\nabla}h_{i}(x_{i}^{*}))^{\top} of hi(xi)h_{i}(x_{i}) at xi=xi∗x_{i}=x_{i}^{*}. Noting that the condition 0∈∂h(x∗)+1(y∗)⊤\mathbf{0}\in\partial\mathbf{h}(\mathbf{x}^{*})+\mathbf{1}(\mathbf{y}^{*})^{\top} is equivalent to the requirement that there exists an x∗\mathbf{x}^{*} such that 0=∇~h(x∗)+1(y∗)⊤\mathbf{0}=\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{*})+\mathbf{1}(\mathbf{y}^{*})^{\top}, the optimality condition for the primal-dual pair (x∗,y∗)(\mathbf{x}^{*},\mathbf{y}^{*}) can be written as:

It turns out that the optimality conditions for the resource optimization problem, as given in Lemma 1, have an interesting connection with the optimality conditions for the consensus optimization problem. In order to expose this relation, we next discuss the consensus optimization problem, which is given as follows:

The local objective of each agent in (10) is the same as that in (6). Unlike the resource allocation problem, instead of having ∑i=1n(xi−ri)=0\sum_{i=1}^{n}(x_{i}-r_{i})=0 as constraints, here we have the consensus constraints, i.e., x1=x2=⋯=xnx_{1}=x_{2}=\cdots=x_{n}.

The first-order optimality condition of (10) is stated in the following lemma.

where UU is the matrix defined in Assumption 2.

The proof for this lemma is basically the same to that for Lemma 3.1 of reference only that we use the decomposition \L=U⊤U\textbf{{\L}}=U^{\top}U while reference uses U=\LU=\sqrt{\textbf{{\L}}}.

II-C The mirror relationship

It is known that the Lagrangian dual problem of the resource allocation problem is a consensus optimization problem (see the discussion around equations (4)∼\sim(6) in ). As having been pointed out in reference , a distributed optimization method that can solve the consensus optimization problem may also be used for the resource allocation problem through solving the dual of the resource allocation problem. Here, we will provide more special relations these two problems have which leads to a class of resource allocation algorithms following the design of a class of decentralized consensus optimization algorithms. Due to such special relations, it is possible that one can give a decentralized resource allocation algorithm without investigating the Lagrangian dual relationship between the above mentioned two problems.

The optimality conditions of the resource allocation problem (6) and the consensus optimization problem (10) are summed up in the following box:

These conditions share the same structure, i.e.,

Furthermore, to analyze a consensus convex optimization algorithm, the most crucial relation we need is the monotone inequality, namely, 0≤⟨∇~h(xk)−∇~h(x∗),xk−x∗⟩0\leq\langle\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{*}),\mathbf{x}^{k}-\mathbf{x}^{*}\rangle for all kk, which translates into verifying 0≤⟨A(xk)−A(x∗),B(xk)−B(x∗)⟩0\leq\langle\mathcal{A}(\mathbf{x}^{k})-\mathcal{A}(\mathbf{x}^{*}),\mathcal{B}(\mathbf{x}^{k})-\mathcal{B}(\mathbf{x}^{*})\rangle for all kk, when we substitute the key quantities by the images of those general maps we have discussed above. Apparently, this inequality still holds when we let A(x)=x−r\mathcal{A}(\mathbf{x})=\mathbf{x}-\mathbf{r} and B(x)=∇~h(x)\mathcal{B}(\mathbf{x})=\widetilde{\nabla}\mathbf{h}(\mathbf{x}) and assume h(x)\mathbf{h}(\mathbf{x}) is convex. It is possible that such substitutions of A\mathcal{A} and B\mathcal{B} do not affect the validity of some of the existing analyses for certain algorithms. For example, the subgradient form of the proximal method is xk+1=xk−α∇~h(xk+1)\mathbf{x}^{k+1}=\mathbf{x}^{k}-\alpha\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}) where α\alpha is a step size. This method can be proven to have ∇~h(xk+1)→0\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})\rightarrow\mathbf{0} if h(x)\mathbf{h}(\mathbf{x}) is convex. Its counterpart after the substitution is ∇~h(xk+1)=∇~h(xk)−α(xk+1−r)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})=\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})-\alpha(\mathbf{x}^{k+1}-\mathbf{r}), which can be resolved as

and updating according to the following rules:

It can be shown that the sequence {xk}\{\mathbf{x}^{k}\} of such an iterative algorithm will converge to r\mathbf{r}, corresponding to the fact that {∇~h(xk)}\{\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})\} will converge to 0\mathbf{0} in the proximal method.

These observations motivate the class of algorithms that we propose for solving resource allocation problems based on some existing consensus optimization algorithms. The consensus optimization algorithms that will be exploited in this paper are simple, efficient, and have recently been further accelerated akin to Nesterov’s fast methods , as well as enhanced to work over asynchronous and directed communication networks . Our algorithm design philosophy also implies possibilities of enhancing the resource allocation algorithms proposed in this paper by using techniques from those for advancing consensus optimization algorithms.

In Section III, we describe our resource allocation algorithms and provide their convergence analysis. Finally, we will illustrate some numerical experiments in Section VI and conclude the paper with remarks in Section VII.

III The Algorithms and their Convergence Analysis

Before we introduce our algorithms and conduct the analyses, let us introduce the solution set of the resource allocation problem (1), denoted by X∗\mathcal{X}^{*}. We make the following assumption for problem (1), which we use throughout the paper.

The set X∗\mathcal{X}^{*} is nonempty, for example, when the constraint set of the resource allocation problem (1) is compact, or the objective function satisfies some growth condition. Under the convexity conditions in Assumption 1, the optimal set X∗\mathcal{X}^{*} is closed and convex. Under Assumptions 1 and 3, the strong duality holds for problem (1) and its Lagrangian dual problem, and the dual optimal set is nonempty (see Proposition 6.4.2 of ).

This algorithm solves the original problem (1), i.e., the resource allocation with local constraints. This basic algorithm operates as follows (Algorithm 1). Each agent ii uses its local parameter βi>0\beta_{i}>0, which can be viewed as the stepsize.

The algorithm is motivated by the P-EXTRA algorithm from reference for a consensus optimization problem. The reason we refer to Algorithm 1 as Mirror-P-EXTRA will be clear from the following lemma. In the lemma and later on, we will use BB to denote the diagonal matrix that has βi\beta_{i} as its (i,i)(i,i)-th entry,

Let Assumptions 1 and 2 be satisfied, and let c>0c>0. Then, the sequence {xk,yk}\{\mathbf{x}^{k},\mathbf{y}^{k}\} generated by Algorithm 1 satisfies for k=0,1,…k=0,1,\ldots,

where x0\mathbf{x}^{0} is chosen the same as that in Algorithm 1, q−1=0\mathbf{q}^{-1}=\mathbf{0}, and q0=U∇~h(x0)\mathbf{q}^{0}=U\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{0}) where the matrix ∇~h(x0)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{0}) of subgradients is the same as those used in Algorithm 1.

By using the notation given in Subsection II-A, the updates of Algorithm 1 can be represented compactly as initializing with arbitrary x0∈Ω\mathbf{x}^{0}\in{\bm{\Omega}}, s0=∇~f(x0)\mathbf{s}^{0}=\widetilde{\nabla}\mathbf{f}(\mathbf{x}^{0}), and y−1=0\mathbf{y}^{-1}=\mathbf{0}, and then performing for k=0,1,2,…k=0,1,2,\ldots,

From (13b) and (13c), also considering the initialization Since x0∈Ω\mathbf{x}^{0}\in{\bm{\Omega}}, we can choose the subgradient ∇~g(x0)=0\widetilde{\nabla}\mathbf{g}(\mathbf{x}^{0})=\mathbf{0} and use the relation ∇~f(x0)+∇~g(x0)=∇~h(x0)\widetilde{\nabla}\mathbf{f}(\mathbf{x}^{0})+\widetilde{\nabla}\mathbf{g}(\mathbf{x}^{0})=\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{0}) (see (2)). s0=∇~f(x0)=∇~f(x0)+∇~g(x0)=∇~h(x0)\mathbf{s}^{0}=\widetilde{\nabla}\mathbf{f}(\mathbf{x}^{0})=\widetilde{\nabla}\mathbf{f}(\mathbf{x}^{0})+\widetilde{\nabla}\mathbf{g}(\mathbf{x}^{0})=\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{0}), we have sk=∇~h(xk)\mathbf{s}^{k}=\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k}) for k=0,1,…k=0,1,\ldots. Thus, by substituting sk=∇~h(xk)\mathbf{s}^{k}=\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k}) in (13a), we obtain that the given conditions are equivalent to:

Now, in relation (14b), we write 2cyk−cyk−1=cyk+c(yk−yk−1)2c\mathbf{y}^{k}-c\mathbf{y}^{k-1}=c\mathbf{y}^{k}+c(\mathbf{y}^{k}-\mathbf{y}^{k-1}) and use yk−yk−1=\L∇~h(xk)\mathbf{y}^{k}-\mathbf{y}^{k-1}=\textbf{{\L}}\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k}) (cf. (14a)), to obtain the following equivalent relations:

These relations are enough to generate the sequence {xk,yk}\{\mathbf{x}^{k},\mathbf{y}^{k}\}. By Assumption 2, we have that \L=U⊤U\textbf{{\L}}=U^{\top}U, so by introducing the notation yk=U⊤qk\mathbf{y}^{k}=U^{\top}\mathbf{q}^{k}, relation (15a) reduces to qk+1=qk+U∇~h(xk+1)\mathbf{q}^{k+1}=\mathbf{q}^{k}+U\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}). We note that a sequence {yk}\{\mathbf{y}^{k}\} generated by yk=U⊤qk\mathbf{y}^{k}=U^{\top}\mathbf{q}^{k} is the same to the sequence {yk}\{\mathbf{y}^{k}\} generated by (15a) with y−1=0\mathbf{y}^{-1}=\mathbf{0}. Using these relations and reorganizing (15), we obtain

which generates the same {xk,yk}\{\mathbf{x}^{k},\mathbf{y}^{k}\} sequence as Algorithm 1 does. Finally, relation (16b) is equivalent to qk=qk+1−U∇~h(xk+1)\mathbf{q}^{k}=\mathbf{q}^{k+1}-U\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}), which when substituted into (16a) gives relation (12a). Relations (16b) and (16c) coincide with (12b) and (12c), respectively.

To fulfill the optimal condition for the resource allocation problem, we would need ∇~h(xk+1)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}) to be consensual and the rows of xk+1−r\mathbf{x}^{k+1}-\mathbf{r} to sum up to 00. In view of the insight from Lemmas 1 and 2, the only thing we need to do is to replace ∇~h(xk+1)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}) by xk+1−r\mathbf{x}^{k+1}-\mathbf{r} and replace xk+1\mathbf{x}^{k+1} by ∇~h(xk+1)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}), as discussed in Section II-C. By doing so, we obtain

which is very similar to the relations (12a) and (12b) in Lemma 3. The only difference is in the term cW~(∇~h(xk+1)−∇~h(xk))c\widetilde{W}(\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})) of (18) whereas we have “replaced” this term by (B−c\L)(∇~h(xk+1)−∇~h(xk))(B-c\textbf{{\L}})(\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})) to obtain (12a). The key function of the term cW~(xk+1−xk)c\widetilde{W}(\mathbf{x}^{k+1}-\mathbf{x}^{k}) in (17) (corresponding to the term cW~(∇~h(xk+1)−∇~h(xk))c\widetilde{W}(\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})) in (18)) is to stabilize the iterative process and neutralize those terms that are not one-step decentralized implementable in P-EXTRA. Hypothetically, with a “large enough” positive (semi)definite matrix PP, any term P(xk+1−xk)P(\mathbf{x}^{k+1}-\mathbf{x}^{k}) in (17) (or P(∇~h(xk+1)−∇~h(xk))P(\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})) in (18)) will serve the purpose of stabilizing the iterative process. Here, we redesign this term as a more flexible one, (B−c\L)∇~h(xk+1)(B-c\textbf{{\L}})\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}), so that the recursive relations are resolvable and implementable in decentralized manner, while featuring per-agent-independent parameters.

Our analysis of Algorithm 1 will use the alternative description of the algorithm, as given in Lemma 3. To simplify our presentation, let us define the following quantities:

where c>0c>0 is the parameter of Algorithm 1. Using this particular cc in Lemma 1, with each solution x∗∈X∗\mathbf{x}^{*}\in\mathcal{X}^{*} we can identify q∗\mathbf{q}^{*} such that Lemma 1 holds, i.e., the optimality conditions in (7) are satisfied. This particular q∗\mathbf{q}^{*} and the matrix ∇~h(x∗)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{*}) constitute the matrix z∗\mathbf{z}^{*}.

Next, we will show that xk\mathbf{x}^{k} converges to a solution x∗\mathbf{x}^{*}.

Let Assumptions 1–3 hold. Let the parameters βi\beta_{i} and c>0c>0 be such that B−c\L≻0B-c\textbf{{\L}}\succ\mathbf{0}. Then, {\color[rgb]{0,0,0}\mathcal{M}\succ 0,} and the sequences {xk}\{\mathbf{x}^{k}\} and {qk}\{\mathbf{q}^{k}\} generated by Algorithm 1 satisfy the following relations:

where M\mathcal{M}, zk\mathbf{z}^{k} and z∗=z(x∗)\mathbf{z}^{*}=\mathbf{z}(\mathbf{x}^{*}), for x∗∈X∗\mathbf{x}^{*}\in\mathcal{X}^{*}, are defined by (19). Furthermore, the sequence {xk}\{\mathbf{x}^{k}\} converges to a point in the optimal set X∗\mathcal{X}^{*}.

The fact that M≻0\mathcal{M}\succ 0 follows directly from the assumptions of the theorem on the choice of the parameters βi\beta_{i} and c>0c>0. By the convexity of h\mathbf{h}, we have that for any arbitrary x∗∈X∗\mathbf{x}^{*}\in\mathcal{X}^{*},

By Lemma 1, where c>0c>0 is the chosen parameter in the algorithm, from (7a) we have

Using relation (12b) of Lemma 3 and U∇~h(x∗)=0U\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{*})=\mathbf{0} (see Lemma 1), it follows that

Recalling the definitions of M\mathcal{M}, zk\mathbf{z}^{k}, and z∗\mathbf{z}^{*} in (19), by applying the basic equality 2⟨zk+1−zk,M(z∗−zk+1)⟩=∥zk−z∗∥M2−∥zk+1−z∗∥M2−∥zk−zk+1∥M22\langle\mathbf{z}^{k+1}-\mathbf{z}^{k},\mathcal{M}(\mathbf{z}^{*}-\mathbf{z}^{k+1})\rangle=\|\mathbf{z}^{k}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}-\|\mathbf{z}^{k+1}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}-\|\mathbf{z}^{k}-\mathbf{z}^{k+1}\|_{\mathcal{M}}^{2}, from the preceding inequality we obtain

Recalling that {∇~h(x)}\{\widetilde{\nabla}\mathbf{h}(\mathbf{x})\} is the matrix with the iith row consisting of a subgradient ∇~hi(xi)\widetilde{\nabla}h_{i}(x_{i}) of hih_{i} at xix_{i}, we have for all ii,

To establish the rate of convergence, we use a convergence property of nonnegative monotonic scalar sequence, which has appeared in recent work . The result is stated in the following:

To simplify the notation, let us define Δxk+1≜xk−xk+1\Delta\mathbf{x}^{k+1}\triangleq\mathbf{x}^{k}-\mathbf{x}^{k+1}, Δqk+1≜qk−qk+1\Delta\mathbf{q}^{k+1}\triangleq\mathbf{q}^{k}-\mathbf{q}^{k+1}, Δzk+1≜zk−zk+1\Delta\mathbf{z}^{k+1}\triangleq\mathbf{z}^{k}-\mathbf{z}^{k+1}, and Δ∇~h(xk+1)≜∇~h(xk)−∇~h(xk+1)\Delta\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})\triangleq\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}). By the convexity of h\mathbf{h}, we have

Next, we use Lemma 3, where by taking the difference of (12a) at the kk-th and the (k+1)(k+1)-at iteration, we obtain

By combining (28) and (29), it follows that

Using the relation (12b) for qk\mathbf{q}^{k}, we see that Δqk+1=−U∇~h(xk+1)\Delta\mathbf{q}^{k+1}=-U\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1}). Thus, we have

which when substituted into (30) and re-arranging some terms yields

By applying the basic equality 2⟨MΔzk+1,Δzk−Δzk+1⟩=∥Δzk∥M2−∥Δzk+1∥M2−∥Δzk−Δzk+1∥M22\langle\mathcal{M}\Delta\mathbf{z}^{k+1},\Delta\mathbf{z}^{k}-\Delta\mathbf{z}^{k+1}\rangle=\|\Delta\mathbf{z}^{k}\|_{\mathcal{M}}^{2}-\|\Delta\mathbf{z}^{k+1}\|_{\mathcal{M}}^{2}-\|\Delta\mathbf{z}^{k}-\Delta\mathbf{z}^{k+1}\|_{\mathcal{M}}^{2} to (31), we finally have

which implies (27) and completes the proof.

Based on Lemma 4, we have the following rate result for the first-order optimality residual.

Under the assumptions of Theorem 1, the first-order optimality residual decays to 00 at an o(1k)o\left(\frac{1}{k}\right) rate, i.e.,

As an immediate consequence of Theorem 1 and Proposition 1, we can see that Lemma 4 holds for the sequence ak=∥zk−zk+1∥M2a_{k}=\|\mathbf{z}^{k}-\mathbf{z}^{k+1}\|_{\mathcal{M}}^{2}, which implies the stated result.

Each function fif_{i} is μi\mu_{i}-strongly convex.

Let us introduce μ=min⁡iμi\mu=\min_{i}\mu_{i} and L=max⁡iLiL=\max_{i}L_{i}.

Note that we have f=h\mathbf{f}=\mathbf{h}, so that by the strong convexity and the gradient Lipschitz property of h\mathbf{h}, it follows that (see Theorem 2.1.11 of )

Using arguments similar to those from (21) to (22), where (21) is replaced by (32), we obtain

To prove the linear convergence, it suffices to show that δ∥zk+1−z∗∥M2≤∥zk−z∗∥M2−∥zk+1−z∗∥M2\delta\|\mathbf{z}^{k+1}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}\leq\|\mathbf{z}^{k}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}-\|\mathbf{z}^{k+1}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2} for some δ>0\delta>0. Considering (33), this will hold as long as for some δ>0\delta>0 and all k≥0k\geq 0 we have

By the definition of zk\mathbf{z}^{k}, we have

where in the last equality we use the fact U∇h(x∗)=0U\nabla\mathbf{h}(\mathbf{x}^{*})=\mathbf{0} (see (7b)). Substituting (35) and (36) in (34), we conclude that for the linear convergence, it is sufficient to show the following relation:

Combining (37) and (41), we find that, in order to establish the rate result for ∥zk+1−z∗∥M2\|\mathbf{z}^{k+1}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}, it suffices to show that

Now, we utilize the basic properties of function h\mathbf{h} to verify the validness of (42). Considering the Lipschitz continuity of ∇h\nabla\mathbf{h}, we only need to determine a positive δ\delta such that

Let us use LHSi\text{LHS}i and RHSi\text{RHS}i (see Section II-A) to refer to key terms appeared in (42) and elaborate how (43) and (44) are obtained. The (matrix) inequality (43a) and (44a) are identical and due to requiring LHS1≤RHS1\text{LHS}1\leq\text{RHS}1; the matrix inequalities (43c) and (43e) are due to using the Lipschitz continuity of the gradient to generate a lower bound for RHS2\text{RHS}2 and requiring LHS2\text{LHS}2 bing bounded by it; the matrix inequalities (44c) and (44e) are due to using the strong convexity of h\mathbf{h} to generate an upper bound for −RHS2-\text{RHS}2, i.e.,

and requiring this upper bound being further bounded by −LHS2-\text{LHS}2.

By re-arranging terms and shrinking the feasible regions of (43) and (44), we re-write these conditions in more compact forms as follows:

For any γ>0\gamma>0, a positive small enough δ\delta exists that satisfies either \delta\leq\delta_{1}=\sup\{\delta>0\mid\delta\text{ satisfies \eqref{eq:geo_p7} }\} or \delta\leq\delta_{2}=\sup\{\delta>0\mid\delta\text{ satisfies \eqref{eq:geo_p8} }\}. Thus, we can find a positive δ≤max⁡{δ1,δ2}\delta\leq\max\{\delta_{1},\delta_{2}\} such that ∥zk+1−z∗∥M2≤11+δ∥zk−z∗∥M2\|\mathbf{z}^{k+1}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}\leq\frac{1}{1+\delta}\|\mathbf{z}^{k}-\mathbf{z}^{*}\|_{\mathcal{M}}^{2}.

In the above proof of Theorem 3, explicit bounds on δ\delta can be easily derived when we further choose B=βIB=\beta I. In this case, a sufficient condition on δ\delta would be either

where γ>0\gamma>0 is a parameter that we can choose. Consequently, we have the following two corollaries.

Since it appears in both (46) that the larger β\beta is, the worse the convergence rate is, let us choose β=c\beta=c.

From (49), we know that to reach ε\varepsilon-accuracy, the number of iterations needed is (omitted the factor ln⁡(1ε)\ln\left(\frac{1}{\varepsilon}\right) for briefness)

Under the same conditions as those in Corollary 1, by letting \L=0.5(I−W)\textbf{{\L}}=0.5(I-W) where WW is a lazy Metropolis matrix, that is,

and by setting the parameters in Mirror-P-EXTRA as c=Θ(nμ−0.5L−0.5)c=\Theta\left(n\mu^{-0.5}L^{-0.5}\right) and β=c\beta=c, to reach ε\varepsilon-accuracy, the number of iterations needed is of the order of O((n2+nκf0.5)ln⁡(ε−1))O\left((n^{2}+n\kappa_{\mathbf{f}}^{0.5})\ln\left(\varepsilon^{-1}\right)\right).

III-B The algorithm specialized for the unconstrained case: Mirror-EXTRA

In this section we concern the resource allocation without local constraints, i.e.,

Our proposed Algorithm 2 applies to (52) when the objective function has Lipschitz gradients, and it is given as follows.

Compared to Algorithm 1, Algorithm 2 has a lower per-iteration cost in the update of xikx_{i}^{k} which requires a gradient evaluation instead of “prox”-type (minimization) update. We next provide an alternative description Algorithm 2 that we use in the analysis of the method.

where x0\mathbf{x}^{0} is the same as in Algorithm 2, q−1=0\mathbf{q}^{-1}=\mathbf{0}, and q0=U∇f(x0)\mathbf{q}^{0}=U\nabla\mathbf{f}(\mathbf{x}^{0}).

Relations (53a), (53c), and (53d) are obtained using the analysis that is similar to that for deriving the corresponding relations in Lemma 3, so we omit their proofs. Relation (53b) is obtained from (53a) by using r=x∗+U⊤q∗\mathbf{r}=\mathbf{x}^{*}+U^{\top}\mathbf{q}^{*} (see Lemma 1).

Comparing the description of Algorithm 2 (Lemma 5) and the description of Algorithm 1 (Lemma 3), we can see that the only difference between these algorithms is in the scaling matrix of the “proximal term” (∇~h(xk+1)−∇~h(xk))(\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k+1})-\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{k})). In Algorithm 1 the coefficient matrix is βI−c\L\beta I-c\textbf{{\L}}, while in Algorithm 2 it is −c\L-c\textbf{{\L}}. We refer to the algorithm Mirror-EXTRA due to its relation to P-EXTRA and the fact that it features a gradient-based update, just as the EXTRA algorithm does for consensus optimization .

We start with a lemma that provides an important relation for the convergence analysis of Mirror-EXTRA.

We start by showing a relation valid for arbitrary vectors x,y,zx,y,z and for any ρ>0\rho>0:

Now, by multiplying the relation (58) with a(1−t)a(1-t) where t∈(0,1)t\in(0,1), we obtain

Letting x=vk\mathbf{x}=\mathbf{v}^{k}, y=vk+1\mathbf{y}=\mathbf{v}^{k+1}, and z=v∗\mathbf{z}=\mathbf{v}^{*}, we have

Next, by adding relations (59) and (54), we find that

Regarding the conditions of Lemma 6 in (55) and (56), we note that if a>3b>0a>3b>0, then it can be verified that the conditions are met by choosing ρ=a−b2b\rho=\frac{a-b}{2b} and any t∈[ba,a−b2a).t\in[\frac{b}{a},\frac{a-b}{2a}).

By the convexity and the gradient Lipschitz continuity of f\mathbf{f}, we have for any x∗∈X∗\mathbf{x}^{*}\in\mathcal{X}^{*},

where L′=Lλmax⁡{\L}L^{\prime}=L\lambda_{\max}\{\textbf{{\L}}\}. From Lemma 5 (cf. (53b)) we have an expression for xk+1−x∗\mathbf{x}^{k+1}-\mathbf{x}^{*}, which when substituted in (61) yields

where the first equality follows from relation (53c) of Lemma 5 and the optimality condition U∇f(x∗)=0U\nabla\mathbf{f}(\mathbf{x}^{*})=0 (cf. Lemma 1 with ∇~h(x∗)=∇f(x∗)\widetilde{\nabla}\mathbf{h}(\mathbf{x}^{*})=\nabla\mathbf{f}(\mathbf{x}^{*})). Therefore,

We now apply Lemma 6 with the following identification: in (63), we set {uk,vk}≜{qk,U∇f(xk)}\{\mathbf{u}^{k},\mathbf{v}^{k}\}\triangleq\{\mathbf{q}^{k},U\nabla\mathbf{f}(\mathbf{x}^{k})\}, the point {u∗,v∗}≜{q∗,U∇f(x∗)}\{\mathbf{u}^{*},\mathbf{v}^{*}\}\triangleq\{\mathbf{q}^{*},U\nabla\mathbf{f}(\mathbf{x}^{*})\} (note that \L=U⊤U\textbf{{\L}}=U^{\top}U), and the quantities b=cb=c, a=2L′−ca=\frac{2}{L^{\prime}}-c and ρ=a−b2b\rho=\frac{a-b}{2b}. The condition c∈(0,12L′)c\in\left(0,\frac{1}{2L^{\prime}}\right) is equivalent to a>3b>0a>3b>0 (see Remark 1). Thus, by applying Lemma 6, we obtain that for any c∈(0,12L′)c\in\left(0,\frac{1}{2L^{\prime}}\right) and t∈[cL′2−cL′,1−cL′2−cL′)t\in\left[\frac{cL^{\prime}}{2-cL^{\prime}},\frac{1-cL^{\prime}}{2-cL^{\prime}}\right),

Choosing t=12(2−cL′)t=\frac{1}{2(2-cL^{\prime})} in (64), and using the preceding two relations, we obtain

By following a nearly identical line of analysis as in the proof of Theorem 1, from relations (23a) and (23b) onward, we can conclude that {xk}\{\mathbf{x}^{k}\} converges to a point in the optimal set X∗\mathcal{X}^{*}.

The bound c≤12Lλmax⁡{\L}c\leq\frac{1}{2L\lambda_{\max}\{\textbf{{\L}}\}} does not necessarily imply that the step size cc selection requires any knowledge of the graph G\mathcal{G} structure. For example, such a requirement can be avoided by using \L=0.5(I−W)\textbf{{\L}}=0.5(I-W), where WW is a symmetric stochastic matrix, in which case one may employ a step size c≤12Lc\leq\frac{1}{2L}.

We next investigate convergence rate properties of Mirror-EXTRA, where we make use of the following result which is based on Lemma 5.

where Δqk+1≜qk−qk+1\Delta\mathbf{q}^{k+1}\triangleq\mathbf{q}^{k}-\mathbf{q}^{k+1} and Δ∇f(xk+1)≜Δ∇f(xk)−Δ∇f(xk+1)\Delta\nabla\mathbf{f}(\mathbf{x}^{k+1})\triangleq\Delta\nabla\mathbf{f}(\mathbf{x}^{k})-\Delta\nabla\mathbf{f}(\mathbf{x}^{k+1}).

To simplify the notation, we also define Δxk+1≜xk−xk+1\Delta\mathbf{x}^{k+1}\triangleq\mathbf{x}^{k}-\mathbf{x}^{k+1}. By the convexity and the Lipschitz continuity of ∇f\nabla\mathbf{f}, we have

where L′=Lλmax⁡{\L}L^{\prime}=L\lambda_{\max}\{\textbf{{\L}}\}. Now we use relation (53a) of Lemma 5; specifically, by taking the difference between (53a) at the kk-th and (k+1)(k+1)-th iteration we obtain

From (66) we have an expression for Δxk+1\Delta\mathbf{x}^{k+1}, which when substituted in relation (65) yields

Relation (53c) of Lemma 5 implies Δqk+1=−U∇f(xk+1)\Delta\mathbf{q}^{k+1}=-U\nabla\mathbf{f}(\mathbf{x}^{k+1}), which in turn gives

By substituting (68) into (67), we can see that

Similarly, since \L=U⊤U\textbf{{\L}}=U^{\top}U, we have

By using the preceding two equalities in (69) and by reorganizing terms, we obtain

where the last inequality follows from ∥Δ∇f(xk)−Δ∇f(xk+1)∥\L2≤2∥Δ∇f(xk)∥\L2+2∥Δ∇f(xk+1)∥\L2.\|\Delta\nabla\mathbf{f}(\mathbf{x}^{k})-\Delta\nabla\mathbf{f}(\mathbf{x}^{k+1})\|_{\textbf{{\L}}}^{2}\leq 2\|\Delta\nabla\mathbf{f}(\mathbf{x}^{k})\|_{\textbf{{\L}}}^{2}+2\|\Delta\nabla\mathbf{f}(\mathbf{x}^{k+1})\|_{\textbf{{\L}}}^{2}. Therefore,

where in the last inequality we use c≤2L′−3cc\leq\frac{2}{L^{\prime}}-3c, which holds by the assumption that c∈(0,12L′)c\in\left(0,\frac{1}{2L^{\prime}}\right).

We have the following basic rate result for the iterates of Algorithm 2, when the objective function f\mathbf{f} is convex and has Lipschitz continuous gradients. The result follows directly from Theorem 4, Lemma 7, and Proposition 1.

Under the assumptions of Theorem 4, along the iterates of Algorithm 2, the first-order optimality residual decays to 00 at an o(1k)o\left(\frac{1}{k}\right) rate, i.e.,

We next show a linear convergence of Mirror-EXTRA under additional strong convexity assumption on the function f\mathbf{f}. To simplify the analysis, we will assume that λmax⁡{\L}≤1\lambda_{\max}\{\textbf{{\L}}\}\leq 1,which holds for example when \L=\LG/λmax⁡{\LG}\textbf{{\L}}={\textbf{\L}}_{\mathcal{G}}/\lambda_{\max}\{{\textbf{\L}}_{\mathcal{G}}\} or \L=0.5(I−W)\textbf{{\L}}=0.5(I-W) for some symmetric stochastic matrix WW compatible with the graph G\mathcal{G}.

where δ>0\delta>0 is such that for some γ>0\gamma>0,

By the strong convexity and the gradient Lipschitz continuity of the function f\mathbf{f}, and the assumption that λmax⁡{\L}≤1\lambda_{\max}\{\textbf{{\L}}\}\leq 1, we have

From relation (53b) of Lemma 5 it follows that

By the optimality condition U∇f(x∗)=0U\nabla\mathbf{f}(\mathbf{x}^{*})=0, it follows that

where the last equality follows from relation (53c) of Lemma 5. Therefore,

Upon substituting relations (73) and (76) in (72), after re-arranging the terms, we obtain

Now, we use (75) and we add (2L−2c)∥∇f(xk)−∇f(x∗)∥\L2−(2L−2c)∥∇f(xk+1)−∇f(x∗)∥\L2(\frac{2}{L}-2c)\|\nabla\mathbf{f}(\mathbf{x}^{k})-\nabla\mathbf{f}(\mathbf{x}^{*})\|_{\textbf{{\L}}}^{2}-(\frac{2}{L}-2c)\|\nabla\mathbf{f}(\mathbf{x}^{k+1})-\nabla\mathbf{f}(\mathbf{x}^{*})\|_{\textbf{{\L}}}^{2} to both sides of the preceding inequality, which gives

In view of the preceding relation, in order to prove the linear convergence, it suffices to show that for some δ>0\delta>0 the following relation holds for all k≥1k\geq 1,

To see this, note that assuming that relation (78) is valid, we will have

showing that the iterate sequence {xk}\{\mathbf{x}^{k}\} converges to the optimal solution x∗\mathbf{x}^{*} at an R-linear rate.

furthermore, it always holds that ∥∇f(xk)−∇f(xk+1)∥\L2≤2∥∇f(xk)−∇f(x∗)∥\L2+2∥∇f(xk+1)−∇f(x∗)∥\L2\|\nabla\mathbf{f}(\mathbf{x}^{k})-\nabla\mathbf{f}(\mathbf{x}^{k+1})\|_{\textbf{{\L}}}^{2}\leq 2\|\nabla\mathbf{f}(\mathbf{x}^{k})-\nabla\mathbf{f}(\mathbf{x}^{*})\|_{\textbf{{\L}}}^{2}+2\|\nabla\mathbf{f}(\mathbf{x}^{k+1})-\nabla\mathbf{f}(\mathbf{x}^{*})\|_{\textbf{{\L}}}^{2}. From the preceding two inequalities it follows that relation (78) is valid as long as the following relation holds

We next further examine some sufficient relations for (80) to be valid. In particular, by the assumption that λmax⁡{\L}≤1\lambda_{\max}\{\textbf{{\L}}\}\leq 1 and by the Lipschitz continuity of ∇f\nabla\mathbf{f}, we have

The preceding relation will hold, as long as δ>0\delta>0 is small enough so that, for some γ>0\gamma>0, we have

Given a γ>0\gamma>0, to satisfy the conditions in (81), one can choose δ>0\delta>0 so that

As a consequence of Theorem 6, we have the following corollary regarding the scalability of the Mirror-EXTRA.

From the bound on δ\delta in Theorem 6 we see that, to reach ε\varepsilon-accuracy, the number of iterations needed is (the factor ln⁡(1ε)\ln\left(\frac{1}{\varepsilon}\right) is omitted for simplicity)

The complexity of the Mirror-EXTRA is slightly worse than the complexity of Mirror-P-EXTRA, which is O((κ\L+κ\Lκf)ln⁡(1ε))O\left((\kappa_{\textbf{{\L}}}+\sqrt{\kappa_{\textbf{{\L}}}\kappa_{\mathbf{f}}})\ln\left(\frac{1}{\varepsilon}\right)\right). The Mirror-EXTRA has a lower per-iteration cost at the expense of less favorable scalability, and O(κ\Lκfln⁡(1ε))O\left(\kappa_{\textbf{{\L}}}\kappa_{\mathbf{f}}\ln\left(\frac{1}{\varepsilon}\right)\right) scalability of Mirror-EXTRA also coincides with that of the algorithm proposed in reference .

In the following subsections, we will first provide a gradient projection algorithm which can be understood as a hybrid of Algorithm 1 and Algorithm 2. Then we will conduct two case studies to show how the methodology of using the “mirror relation” can help to design more powerful distributed resource allocation algorithms based on existing consensus optimization algorithms.

III-C A projection gradient-based algorithm: Mirror-PG-EXTRA

Since Algorithm 1 has to solve a constrained optimization problem per iteration over each agent (which can be costly) while Algorithm 2 cannot be applied to constrained problems, we consider another algorithm that uses a gradient-projection update which can be understood as a hybrid of Algorithms 1 and 2. We name the newly constructed algorithm Mirror-PG-EXTRA (see Algorithm 3) after a previous algorithm for consensus optimization, PG-EXTRA .

In Algorithm 3, the matrix BB is a diagonal matrix with entries Bii=βiB_{ii}=\beta_{i} on its diagonal, while PΩi[⋅]\mathcal{P}_{\Omega_{i}}[\cdot] is the projection operator on the set Ωi\Omega_{i} with respect to the Euclidean norm. In the absence of the per-agent-constraints, Algorithm 3 degenerates to Algorithm 2 (this can be verified by writing out the recursive relations of {xk}\{\mathbf{x}^{k}\}, {yk}\{\mathbf{y}^{k}\}, and {sk}\{\mathbf{s}^{k}\} in Algorithm 3 and eliminating the sequence {sk}\{\mathbf{s}^{k}\} in the relations). We do not formally establish the convergence or the convergence rates for Algorithm 3, though we believe it has 1/k1/k convergence rate under convexity and smoothness assumptions. An immediate guess on the convergence property of Algorithm 3 is that it will be slower than Algorithm 1 in terms of the number of iterations needed to reach a given accuracy. We will evaluate it numerically in our simulations (see Section VI).

IV Improving the rates by Nesterov’s acceleration

To make the discussion concise, without loss of generality, let us choose Ł such that ∥U∥2=σmax⁡{U}=1\|U\|_{2}=\sigma_{\max}\left\{U\right\}=1. This can be done by choosing, for example, \L=0.5(I−W)\textbf{{\L}}=0.5(I-W) where WW is a symmetric doubly stochastic matrix that is compatible with the graph (see Subsection II-A and Corollary 2). Also suppose that Assumption 5 (Lipscthiz continuity of ∇f(⋅)\nabla\mathbf{f}(\cdot)) holds and L=max⁡iLiL=\max_{i}L_{i} is the Lipschitz constant of ∇f(⋅)\nabla\mathbf{f}(\cdot), then the gradient mapping of J\mathbf{J} is given by

and is ∇J(⋅)\nabla\mathbf{J}(\cdot) is LL-Lipschitz as well.

Next we will show that an equivalent form of the Nesterov’s accelerated gradient method on finding the minimizer of J(z)\mathbf{J}(\mathbf{z}) can be implemented in decentralized way and clearly this can be utilized to obtain x∗\mathbf{x}^{*}. The recursion of the Nesterov accelerated gradient method (see for the introduction of the algorithm and its analysis) for minimizing J(z)\mathbf{J}(\mathbf{z}) is

For any zk\mathbf{z}^{k} and yk\mathbf{y}^{k}, let us define

Then from (82), we can obtain the following recursions of xk\mathbf{x}^{k} and vk\mathbf{v}^{k} which inherit the convergence properties of (82):

The algorithm given in (83) can be implemented in a fully decentralized fashion. We thus can apply Nesterov’s accelerated gradient method for the resource allocation problem over a network and obtain improved convergence rates. Specifically, we have the following two propositions for the claims of rates. Proofs are omitted due to the fact they are direct applications of Nesterov’s method.

Suppose that Assumption 4 holds and Ł is chosen such that σmax⁡{\L}=1\sigma_{\max}\left\{\textbf{{\L}}\right\}=1. If we choose α=1L\alpha=\frac{1}{L} and βk=kk+3\beta_{k}=\frac{k}{k+3}, with the algorithm given in (83), it is guaranteed that

To explore the geometric convergence, it is expected that one assumes strong convexity. Suppose that Assumption 5 holds, then for any a\mathbf{a} and b\mathbf{b} with proper dimensions,

Suppose that Assumptions 4 and 5 hold and \L=0.5(I−W)\textbf{{\L}}=0.5(I-W) where WW is the lazy Metropolis matrix. If we choose α=1L\alpha=\frac{1}{L} and βk=β=Q−1Q+1\beta_{k}=\beta=\frac{\sqrt{Q}-1}{\sqrt{Q}+1} where Q=71n^2L/μQ=71\hat{n}^{2}L/\mu, then in O(n^κf0.5ln⁡(ε−1))O\left(\hat{n}\kappa_{\mathbf{f}}^{0.5}\ln(\varepsilon^{-1})\right) iterations, the algorithm will give f(vk)−f(x∗)≤ε\mathbf{f}(\mathbf{v}^{k})-\mathbf{f}(\mathbf{x}^{*})\leq\varepsilon.

The above discussion only applies to problems with smooth objectives assuming no local constraints. Adapting Nesterov’s accelerated proximal gradient method for the resource allocation problem with local constraints will be one of our future directions.

V Extensions: Using the Mirror to Conquer More Complicated Scenarios

In this section, we will conduct two case studies to show how our methodology can help to design more powerful distributed resource allocation algorithms based on existing consensus optimization algorithms.

Before introducing the so called Mirror-Push-DIGing algorithm for solving the decentralized resource allocation problem, let us review the Push-DIGing algorithm which is proposed in reference for solving the decentralized consensus optimization over time-varying directed graphs with geometric convergence guarantees. The procedure of Push-DIGing is as follows:

Since the discussion is under the time-varying set-up, instead of using the right upper corner mark k as previous, we now use (k)(k) to indicate the kk-th iteration as well as the time index kk. In the Push-DIGing algorithm, the aggregated symbols u\mathbf{u}, v\mathbf{v}, x\mathbf{x}, y\mathbf{y}, ∇f(x)\nabla\mathbf{f}(\mathbf{x}) are all defined in a similarly way as we have explained for the quantities x\mathbf{x} and ∇~f(x)\widetilde{\nabla}\mathbf{f}(\mathbf{x}) in Section II-A (each agent ii maintains the ii-th row). The matrix C\mathbf{C} is a column stochastic matrix (namely, each row sums to 11) that admits the topology of the network (directed graph). A popular choice of initialization sets d(0)=y(0)=∇f(x(0))\mathbf{d}(0)=\mathbf{y}(0)=\nabla\mathbf{f}(\mathbf{x}(0)). Using uncoordinated step sizes is possible but we will keep using the same step size α\alpha across agents for the sake of simplicity.

To develop an algorithm for resource allocation from Push-DIGing (also considered as a gradient-based method, or termed as explicit method or forward method in different research contexts), we need to first tweak the recursion to produce a proximal method (also termed as implicit method or backward method in different research contexts). The recursive relations of such proximal variant is given by A proximal-gradient variant can also be derived following a similar idea.

The original Push-DIGing algorithm has adopted the so-termed “Adapt-then-Combine” (ATC) strategy in its u\mathbf{u}-update to accelerate the convergence (see Remark 3 of reference ). But in the proximal variant (89), such strategy cannot be applied due to the implementability (see reference for a similar story of not being able to take advantage of the ATC strategy to accelerate the proximal update). Also, we have replaced ∇f\nabla\mathbf{f} by ∇~h\widetilde{\nabla}\mathbf{h} to allow a non-smooth objective. Finally the recursive relations (89) with a simple initialization can be resolved as follows.

Having the above intuitive introduction on Push-DIGing and its proximal variation, now we are ready to give the Mirror-Push-DIGing algorithm for decentralized resource allocation over time-varying directed graphs. The design of Mirror-Push-DIGing will start from keeping the recursive relation

V-B Dealing with local couplings

A generalization of the resource allocation problem which explicitly considers the local linear coupling of multiple resources is as follows

Note that, unlike previous sections, we are no longer using the aggregated/compact notation in matrix forms since it becomes inconvenient when xix_{i}’s have different dimensions. But still as that in previous sections, the function hi≜fi+gih_{i}\triangleq f_{i}+g_{i} is only assumed to be convex and already include the indicator function of the local constraint xi∈Ωix_{i}\in\Omega_{i} in itself.

To eventually converge to a point that meets (94), we can construct sequences that satisfy the following recursive relations:

Algorithm 5: Mirror-P-EXTRA handling local couplings

It can be seen that this algorithm degenerates to Mirror-P-EXTRA (Algorithm 1) when Ai=IA_{i}=I for all ii. We have not analyzed the convergence properties of this algorithm. We believe that its convergence behavior is similar to that of Mirror-P-EXTRA, since it is a generalization of Mirror-P-EXTRA with an intuitive modification introduced above.

VI Numerical Experiments

where the interval boundaries ω‾i,j\overline{\omega}_{i,j} are randomly generated following the uniform distribution over the interval $forallfor alli\in[n]andandj=1,2$.

To compare the competitiveness of the algorithms of Section III with those in the existing literature, we implement the DPDA-S algorithm of . DPDA-S has one system-level parameter γ\gamma and two per-agent parameters, τi\tau_{i} and κi\kappa_{i}. Based on the recommended parametric structure as given in Remark II.1 of , which sets τi\tau_{i} and κi\kappa_{i} automatically to produce another parameter cic_{i}, we tuned the parameter cic_{i} and γ\gamma to obtain one plot for this algorithm. In another plot for this algorithm, we hand-optimized all the parameters to achieve a better performance. In Mirror-P-EXTRA, there is a system-level parameter cc which can be set based on the per-agent parameters βi\beta_{i}. Compared to Mirror-EXTRA, DPDA-S requires a finer tune in its parameters to achieve a competitive performance. The convergence curves are shown in Fig. 1.

In the above numerical test (see Fig. 1), the outcome of Mirror-P-EXTRA outperforming the other algorithms in the number of iterations is in expectation. Mirror-P-EXTRA has to pay extra computational effort (solving a constrained convex optimization problem) at each iteration compared to other competitive algorithms. However, in decentralized computing, it may worth it to perform more computations before exchanging information in the following round of communication because communication costs and delays are usually considered more significant compared to that caused by local computations. There is actually a trade-off between the total computational cost/time and the total communication cost/delay. Which algorithm is more preferable depends on the hardness of the specific optimization problem and the performance of the underlying cyber system.

VII Conclusion

In this paper, we have presented an interesting relationship between resource allocation problem and the consensus optimization problem. Based on this relation, we have proposed two algorithms, namely Mirror-P-EXTRA and Mirror-EXTRA, for distributed resource allocation in a static connected undirected graph. We have established the convergence and convergence rate properties of the algorithms. In particular, we have shown that both of the algorithms enjoy an RR-linear convergence rate when the resource allocation problem has strongly convex objective function with Lipschitz continuous gradients and does not have additional set constraints. We also have illustrated the convergence behavior of the Mirror-P-EXTRA and its computationally less expensive projection-based variant (Mirror-PG-EXTRA) by some numerical experiments.

References