Decomposing Linearly Constrained Nonconvex Problems by a Proximal Primal Dual Approach: Algorithms, Convergence, and Applications

Introduction

Consider the following optimization problem

Distributed Optimization Over Networks. Consider a network consists of $N$ agents who collectively optimize the following problem

Introduce $N$ local variables $x=[x_{1},\cdots,x_{N}]^{T}$, and suppose the graph \{{\mbox{\mathcal{V}}},\mathcal{E}\} is connected. Then it is clear that the following formulation is equivalent to the global consensus problem, which is precisely problem (1)

Multi-Block Linearly Constrained Problem. Consider the following multi-block linearly constrained problem

Such problem, convex or nonconvex, has wide applications in practice, such as in distributed optimization and coordination (more specifically the sharing problem) , robust Principal Component Analysis and rate maximization problem in downlink broadcast communication channels .

2 Literature Review.

The Augmented Lagrangian (AL) methods, or the methods of multipliers, pioneered by Hestenes and Powell , is a classical algorithm for solving nonlinear nonconvex constrained optimization problems . Many existing packages such as LANCELOT are implemented based on this method. Recently, due to the need to solve very large scale nonlinear optimization problems, the AL and its variants regain their popularity, see recent developments in and the references therein. Also reference have developed an AL based algorithm for nonconvex nonsmooth optimization, where subgradients of the augmented Lagrangian are used in the primal update. When the problem is convex and the constraints are linear, Lan and Monterio have analyzed the iteration complexity for the AL method. More specifically, they have characterized the total number of Nesterov’s optimal iterations that are required to reach high quality primal-dual solutions. However, despite the generality of these methods, it appears that the AL methods does not decompose well over the optimization variables. Further, the AL method, at least in its classical forms, is difficult to be implemented in a distributed manner.

In this paper, we answer the following research question: Is it possible to develop augmented Lagrangian-like decomposition schemes for the linearly constrained nonconvex problem (1), with global convergence rate guarantee. Ideally, the resulting algorithm should be able to decompose the updates of different variable blocks so that each of its steps can be easily implemented in distributed manner and/or in parallel. Further, it is desirable that the resulting algorithm would have global convergence and rate of convergence guarantee. To this end, we study a primal-dual algorithm, where the primal step minimizes certain approximation of the augmented Lagrangian of problem (1), and the dual step performs an approximate dual ascent. The approximation used in the primal step is able to decompose the variables, making it possible to obtain simple subproblems by leveraging the problem structures. Theoretically, we show that whenever the penalty parameter in the augmented Lagrangian is larger than a given threshold, the Prox-PDA converges to the set of stationary solutions, globally and in a sublinear manner (i.e., certain measure of stationarity decreases in the rate of {\mbox{\mathcal{O}}}(1/r), where $r$ is the iteration counter). We also analyze various different extensions of the algorithm, and discuss their applications to the distributed nonconvex optimization problem (3).

The Proposed Algorithm

The proposed algorithm builds upon the classical augmented Lagrangian method (also known as the method of multipliers) . Let us introduce the augmented Lagrangian for problem (1) as

In Prox-PDA, the primal iteration (6a) minimizes the augmented Lagrangian plus a proximal term $\frac{\beta}{2}\|x-x^{r}\|^{2}_{B^{T}B}$. It is important to note that the proximal term is critical in both the algorithm implementation and the analysis. It is used to ensure the following key properties:

The primal problem is strongly convex, hence easily solvable;

The primal problem is decomposable over different variable blocks.

To see why the first point above is possible, suppose $B^{T}B$ is chosen such that $A^{T}A+B^{T}B\succeq I$, and that $f(x)$ has Lipschitz gradient. Then by a result in [45, Theorem 2.1], we know that for any $\beta>L$, the objective function of the $x$-problem (6a) is strongly convex.

The Convergence Analysis

In this section we provide convergence analysis for Algorithm 1. The key is the construction of a new potential function that decreases at every iteration of the algorithm. The constructed potential function is a conic combination the augmented Lagrangian, certain proximal term as well as the size of the violation of the equality constraint, thus it measures the progress of both the primal and dual updates.

The function $f(x)$ is differentiable and has Lipschitz continuous gradient, i.e.,

Further assume that $A^{T}A+B^{T}B\succeq I_{N}.$

There exists a constant $\delta>0$ such that

Without loss of generality and for the simplicity of notations, in the following we will assume that $\underline{f}=0$ We note that this is without loss of generality because we have assumed that $\underline{f}$ is finite, therefore we can consider an equivalent problem with the objective function $f(x)+\underline{f}$, which is always lower bounded by ..

Below we provide a few nonconvex smooth functions $f(x)$ that satisfy Assumption [A1] – [A3]. Note that the first three nonconvex functions are of particular interest in learning neural networks, as they are commonly used as activation functions.

The sigmoid function. The sigmoid function is given by

Clearly it satisfies [A2]. We have $\mbox{sigmoid}^{\prime}(x)=\frac{e^{-x}}{(1+e^{-x})^{2}}\in[0,\;1/4]$, and such boundedness of the first order derivative implies that [A1] is true (by first-order mean value theorem).

The $\arctan$ function. Note that $\arctan(x)\in$, so it clearly satisfies [A2]. $\arctan^{\prime}(x)=\frac{1}{x^{2}+1}\in[0,\;1]$ so it is bounded, which implies that [A1] is true. Finally note that

Therefore the function satisfies [A1]–[A2].

The logit function. Since the logistic function is related to the $\tanh$ function as follows

then Assumptions [A1]-[A2] are again satisfied.

The $\log(1+x^{2})$ function. This function has applications in structured matrix factorization . The function itself is obviously nonconvex and lower bounded. Its first order derivative is $\log^{\prime}(1+x^{2})=\frac{2x}{1+x^{2}}$ and it is also bounded.

The quadratic function $x^{T}Qx$. Suppose that $Q$ is a symmetric matrix but not necessarily positive semidefinite, and suppose that $x^{T}Qx$ is strongly convex in the null space of $A^{T}A$. Then it can be shown that there exists a $\delta$ large enough such that [A2] is true; see e.g., .

Other relevant functions include $\sin(x)$, $\mbox{sinc}(x)$, $\cos(x)$ and so on.

2 The Analysis Steps

Below we provide the analysis of Prox-PDA. Our analysis consists of a series of lemmas leading to a theorem characterizing the convergence and iteration complexity for prox-PDA.

In the first step we provide a bound of the size of the constraint violation using a quantity related to the primal iterates. Let $\sigma_{\min}(A^{T}A)$ denote the smallest non-zero eigenvalue for $A^{T}A$. We have the following result.

Suppose Assumptions [A1] and [A3] are satisfied. Then the following is true for Prox-PDA.

Proof. From the optimality condition of the $x$ problem (6a) we have

Note that by Assumption [A3], we have that $b$ lies in the column space of $A$. Therefore we must have

That is, the difference of the dual variable lies in the column space of $A$. Also from the fact that $\mu^{0}=0$, we have that the dual variable itself also lies in the column space of $A$

Applying these facts to (8), and let $\sigma_{\min}(A^{T}A)$ denote the smallest non-zero eigenvalue of $A^{T}A$, we have

This inequality combined with (8) implies that

Squaring both sides and dividing by $\beta$, we obtain the desired result. Q.E.D.

Our second step bounds the descent of the augmented Lagrangian.

Suppose Assumptions [A1] and [A3] are satisfied. Then the following is true for Algorithm 1

Proof. Since $f(x)$ has Lipschitz continuous gradient, and that $A^{T}A+B^{T}B\succeq I$ by Assumption [A1], it is known that if $\beta>L$, then the $x$-subproblem (6a) is strongly convex with modulus $\gamma:=\beta-L>0$; cf. [45, Theorem 2.1]. That is, we have

where in $(i)$ we have used (3.2) with the identification $z=x^{r+1}$ and $x=x^{r}$; in $(ii)$ we have used the optimality condition for the $x$-subproblem (6a). The claim is proved. Q.E.D.

A key observation from Lemma 3.2 is that no matter how large $\beta$ is, the rhs of (10) cannot be made negative, as the second term is increasing in $\beta$. This observation suggests that the augmented Lagrangian alone cannot serve as the potential function for Prox-PDA.

In search for an appropriate potential function, we need a new object that is decreasing in the order of $\beta\left\|(x^{r+1}-x^{r})-(x^{r}-x^{r-1})\right\|_{B^{T}B}^{2}$. The following lemma shows that the descent of the sum of the constraint violation $\|Ax^{r+1}-b\|^{2}$ and the proximal term $\|x^{r+1}-x^{r}\|_{B^{T}B}^{2}$ has the desired term.

Suppose Assumption [A1] is satisfied. Then the following is true

Proof. From the optimality condition of the $x$-subproblem (6a) we have

Plugging $x=x^{r}$ into the first inequality and $x=x^{r+1}$ into the second, adding the resulting inequalities and utilizing the $\mu$-update step (6b) we obtain

Let us bound the lhs and the rhs of (3.2) separately.

First the lhs of (3.2) can be expressed as

Second we have the following bound for the rhs of (3.2)

It is interesting to observe that the new object, ${\beta}/{2}\left(\|Ax^{r+1}-b\|^{2}+\|x^{r+1}-x^{r}\|^{2}_{B^{T}B}\right)$, increases in $L\|x^{r+1}-x^{r}\|^{2}$ and decreases in $\frac{\beta}{2}\left(\|(x^{r}-x^{r-1})-(x^{r+1}-x^{r})\|_{B^{T}B}^{2}\right)$, while the augmented Lagrangian behaves in an opposite manner (cf. Lemma 3.2). More importantly, in our new object, the constant in front of $\|x^{r+1}-x^{r}\|^{2}$ is independent of $\beta$. Although neither of these two objects decreases by itself, quite surprisingly, a proper conic combination of these two objects decreases at every iteration of Prox-PDA. To precisely state the claim, let us define the potential function for Algorithm 1 as

where $c>0$ is some constant to be determined later. We have the following result.

Suppose the assumptions made in Lemmas 3.1 – 3.3 are satisfied. Then we have the following estimate

Proof. Multiplying both sides of (3.3) by the constant $c$ and then add them to (10), we obtain

From the above analysis, it is easy to see that as long as $c$ and $\beta$ are chosen large enough, the potential function decreases at each iteration of Prox-PDA. Below we derive the precise bounds for $c$ and $\beta$.

First, it is clear that a sufficient condition for $c$ is

Note that the term“$\delta/L$” (Defined in Assumption [A2]) in the $\max$ operator is needed for later use. Importantly, such bound on $c$ is independent of $\beta$.

Second, for any given $c$, we need $\beta$ to satisfy

Clearly combining the bounds for $\beta$ and $c$ we see that $\beta>\delta$. We conclude that if both (19) and (20) are satisfied, then the potential function $P(x^{r+1},x^{r},\mu^{r+1})$ decreases at every iteration.

Our next step shows that by using the particular choices of $c$ and $\beta$ in (19) and (20), the constructed potential function is lower bounded.

Suppose [A1] - [A3] are satisfied, and $(c,\beta)$ are chosen according to (19) and (20). Then the following statement holds true

Proof. To prove this we need to utilize the boundedness assumption in [A2].

First, we can express the augmented Lagrangian function as following

Therefore, summing over $r=1\cdots,T$, we obtain

Suppose Assumption [A2] is satisfied and $\beta$ is chosen according to (20) and (19), then clearly the above sum is lower bounded since

This fact implies that the sum of the potential function is also lower bounded (note, the remaining terms in the potential function are all nonnegative), that is

Note that if $c$ and $\beta$ are chosen according to (19) and (20), then $P_{c,\beta}(x^{r+1},x^{r},\mu^{r+1})$ is nonincreasing. Combined with the lower boundedness of the sum of the potential function, we can conclude that the following is true

Now we are ready to present the main result of this section on the convergence and the rate of convergence of Prox-PDA. To this end, define $Q(x^{r+1},\mu^{r+1})$ as the optimality gap of problem (1), given by

It is easy to see that $Q(x^{r+1},\mu^{r})\to 0$ implies that the limit point $(x^{*},\mu^{*})$ is a KKT point of (1) that satisfies the following conditions

To see this, we can first observe that $Ax^{r+1}-b\to 0$, implying $\mu^{r+1}-\mu^{r}\to 0$, therefore the second condition in (24) hold. Second, using the fact that the first term in the optimality gap also goes to zero, we have

Therefore the first condition in (24) is true.

In the following result we show that the optimality gap $Q(\cdot)$ not only decreases to zero, but does so in a sublinear manner. This is the main result of this section.

Suppose Assumption A is satisfied. Further suppose that the conditions on $\beta$ and $c$ in (20) and (19) are satisfied. Then we have the following claims for the sequence generated by Prox-PDA.

(Eventual Feasibility). The constraint is satisfied in the limit, i.e.,

(Convergence to Stationary Points). Every limit point of the iterates $\{x^{r},\mu^{r}\}$ generated by Algorithm 1 converges to a stationary point of problem (1). Further, $Q(x^{r+1},\mu^{r})\to 0$.

(Sublinear Convergence Rate). For any given $\varphi>0$, let us define $T$ to be the first time that the optimality gap reaches below $\varphi$, i.e.,

Then there exists a constant $\nu>0$ such that the following is true

That is, the optimality gap $Q(x^{r+1},\mu^{r})$ converges sublinearly.

Proof. First we prove part (1). Combining Lemmas 3.4 and 3.5, we conclude that $\|x^{r+1}-x^{r}\|^{2}\to 0$. Then according to (3.1), in the limit we have $\mu^{r+1}\to\mu^{r}$, or equivalently $Ax^{r}\to b$. That is, the constraint violation will be satisfied in the limit.

Then we prove part (2). From the optimality condition of $x$-update step (6a) we have

Then we argue that $\{x^{r}\}$ is bounded if $f(x)+\frac{\beta}{2}\|Ax-b\|^{2}$ is coercive. Note that the potential function can be expressed as

and by our analysis in Lemma 3.5 we know that it is decreasing thus upper bounded. Suppose that $\{x^{r}\}$ is unbounded and let $\mathcal{K}$ denote an infinite subset of iteration index in which \lim_{r\in{\mbox{\mathcal{K}}}}x^{r}=\infty. Passing limit to $P_{c,\beta}(x^{r+1},x^{r},\mu^{r+1})$ over $\mathcal{K}$, and using the fact that $x^{r+1}\to x^{r}$, $\mu^{r+1}\to\mu^{r}$, we have

where the last equality comes from the coerciveness assumption. This is a contradiction to the fact that the potential function $P_{c,\beta}(x^{r+1},x^{r},\mu^{r+1})$ is upper bounded. This concludes the proof for the second part of the result.

Then we prove part (3). Let $\cal{K}$ denote any converging infinite iteration index such that $\{(\mu^{r},x^{r})\}_{r\in\cal{K}}$ converges to the limit point $(\mu^{*},x^{*})$. Passing limit in $\cal{K}$, and using the fact that $\|x^{r+1}-x^{r}\|\to 0$, we have

Combined with the fact that $Ax^{*}-b=0$, we conclude that $(\mu^{*},x^{*})$ is indeed a stationary point of the original problem (1), satisfying (24).

Additionally, even if the sequence $\{x^{r+1},\mu^{r+1}\}$ does not have a limit point, from part (1) we still have $\|\mu^{r+1}-\mu^{r}\|\to 0$ and $\|x^{r}-x^{r+1}\|\to 0$. Hence

where ${\rm(i)}$ is from the optimality condition of the $x$-subproblem (6a). Therefore we have that $Q(x^{r+1},\mu^{r})\to 0.$

Finally we prove part (4). Our first step is to bound the size of the gradient of the augmented Lagrangian. From the optimality condition of the $x$-problem (6a), we have

Therefore, by utilizing the estimate (3.1), we see that there must exist two constants $\xi_{1},\xi_{2}>0$ such that the following is true

From the descent estimate (10) we see that there must exist two constants $\nu_{1},\nu_{2}>0$ such that

Summing over $r$, and let $T$ denote the first time that $Q(x^{r+1},x^{r},\mu^{r+1})$ reaches below $\varphi$, we obtain

We conclude that the convergence in term of the optimality gap function $Q(x^{r+1},\mu^{r})$ is sublinear. Q.E.D.

Our result suggests that depending on the property of $f(x)$ and $\nabla f(x)$, the iterates $\{x^{r+1},\mu^{r+1}\}$ may or may not be bounded. However the optimality measure $Q(x,\mu)$ always converges to zero in a sublinear manner. Note that such sublinear complexity bound is in fact tight, even when applying first-order methods for nonconvex unconstrained problems; see the related discussions in .

Before leaving this section, we remark that a few recent works have analyzed the convergence of a family of splitting method for certain nonconvex problems (which does not cover our formulation (1)). All these works have used the augmented Lagrangian function as the potential function – a technique first developed in . Unfortunately this technique fails to apply to our algorithm because it appears difficult, if not impossible, to show that the augmented Lagrangian alone achieves the desired descent (cf. Lemma 3.2).

Discussion: The Convex Case

It is interesting to observe that the proof techniques used in the previous section apply to the convex case as well – only that for the convex case much milder conditions are required. That is, besides positivity, no additional requirement is needed for the penalty parameter $\beta$. This observation also suggests that the proof techniques used here are rather “tight”, in the sense that it would be difficult to further sharpen the bounds on $\beta$ for the nonconvex case. It is interesting to observe that for convex cases Prox-PDA is closely related to the Method of Multipliers , except that a proximal term is used in the $x$-step (6a). Our analysis below shows that this type of method converges sublinearly for convex problems, without taking any averaging on the iterates, and for arbitrary choice of the positive penalty parameter.

Below we briefly highlight the proof steps of the convex case. Throughout this section, we will assume that Assumption [A1]-[A3] hold true. Additionally assume that $f(x)$ is convex.

First it is easy to bound the difference of the dual variables as

The key difference compared with the proof in Lemma 3.1 is that we do not use the Lipschitz continuity of $\nabla f(x)$ to make the rhs explicitly dependent on $\|x^{r+1}-x^{r}\|$.

Similarly as in Lemma 3.2, the descent of the augmented Lagrangian can be bounded by

​​Note that in the first inequality we have replaced $\gamma$ with $\beta$ due to the convexity assumption of $f(x)$ as well as the assumption that $A^{T}A+B^{T}B\succeq I$.

where the last inequality is true due to the fact that the Lipschitz continuity of $\nabla f(x)$ and the convexity of $f(x)$ implies

For any fixed $\beta>0$ pick $c$ such that

then the potential function will decrease at every iteration of the algorithm, where the descent quantity is composed of a conic combination of the following terms: $-\|x^{r+1}-x^{r}\|^{2}$, $-\|\nabla f(x^{r})-\nabla f(x^{r+1})\|$, $-\|(x^{r}-x^{r-1})-(x^{r+1}-x^{r})\|_{B^{T}B}^{2}$ and $-\|A(x^{r+1}-x^{r})\|^{2}$.

The rest of the proof of the (rate of) convergence follows the similar arguments as those leading to Theorem 3.1. To summarize this section, we provide the following corollary.

Suppose Assumption [A1] - [A3] are satisfied. Suppose that $f(x)$ is convex, and $\beta$ is any positive number. Then the same conclusions in Theorem 3.1 hold true for the Prox-PDA.

Extension: Inexactly Solving the Primal Problems

In this section, we discuss two important extensions of the Prox-PDA, in which the $x$-problem (6a) is solved inexactly. Our first extension replaces the $x$-step by a single gradient-type step, while our second algorithm solves the $x$ problem (6a) to some predefined error. The motivation is that for many practical applications, exactly minimizing the augmented Lagrangian may not be easy.

The proposed proximal gradient primal dual algorithm (Prox-GPDA) replaces the objective function $f(x)$ by the the surrogate function

The detailed algorithm is given in the following table.

Our second extension solves the $x$-step (6a) to certain $\epsilon$-optimality, where $\epsilon$ is some error term with small magnitude (the precise condition will be presented shortly). The quality of the solution is measured by the size of the gradient of the objective function for problem (6a).

Algorithm 3. The Inexact Proximal Primal Dual Algorithm (In-Prox-PDA) Initialize $\mu^{0}$ and $x^{0}$; At each iteration $r+1$, update variables by: $\displaystyle\mbox{Find}\;x^{r+1}\;\mbox{that satisfies the following}$ $\displaystyle\nabla f(x^{r+1})+A^{T}\mu^{r}+\beta A^{T}(Ax^{r+1}-b)+\beta B^{T}B(x^{r+1}-x^{r})=\epsilon^{r+1};$ (30a) $\displaystyle\mu^{r+1}=\mu^{r}+\beta(Ax^{r+1}-b).$ (30b)

The analysis of these two cases follows similar steps as that for Prox-PDA. For Prox-GPDA, the major difference is that there are several places in which we need to bound the term $\|\nabla f(x^{r-1})-\nabla f(x^{r})\|$ instead of $\|\nabla f(x^{r+1})-\nabla f(x^{r})\|$. Moreover, the potential function is no longer decreasing at each iteration. For the In-Prox-PDA case, an explicit condition on the size of the error sequence $\{\epsilon^{r+1}\}$ is needed. In the next subsection we provide an outline of the proof.

First, following the derivation leading to (3.1) we obtain

Note that the first term is now related to the difference squared of the previous two iterations.

Following the proof steps in Lemma 3.2, the descent of the augmented Lagrangian is given by

In the third step we have the following estimate

Note that the first two terms come from the following estimate

In the fourth step we have the following overall descent estimate

Note that there is a slight difference between this descent estimate and our previous estimate (3.4), because now there is a positive term in the rhs, which involves $\|x^{r}-x^{r-1}\|^{2}$. Therefore the potential function is difficult to decrease by itself. Fortunately, such extra term can be bounded by the descent of the previous iteration. We can take the summation over all the iterations and obtain

Clearly as long as the potential function is lower bounded, we have $x^{r+1}\to x^{r}$ and $x^{r+1}-x^{r}\to x^{r}-x^{r-1}$. The rest of the proof follows similar steps leading to Theorem 3.1, hence is omitted.

Suppose Assumption [A1]–[A3] are satisfied. Suppose $\beta$ and $c$ satisfy (20) and (19). Then all the conclusions in Theorem 3.1 hold true for the Prox-GPDA.

We remark that we can replace the approximation function $u(x,x^{r})$ by a larger family of upper-bound functions, following the BSUM (Block Successive Upper Bound Minimization) framework . The analysis follows similar lines of argument presented in this section.

2 The Analysis Outline for In-Prox-PDA

Similarly as in Lemma 3.1, we have the following bound

Then following the proof steps in Lemma 3.2, the descent of the augmented Lagrangian is given by

The inexactness update results in two additional terms in the descent of the augmented Lagrangian.

In the third step we have the following estimate

In the fourth step we have the following overall descent estimate

Therefore as long as the potential function is lower bounded, and that the error sequence satisfies

we have $x^{r+1}\to x^{r}$ and $x^{r+1}-x^{r}\to x^{r}-x^{r-1}$. The rest of the proof follows similar steps leading to Theorem 3.1, hence are omitted. We have the following convergence for In-Prox-GPDA.

Suppose Assumption [A1]–[A3] are satisfied. Suppose $\{\epsilon^{r}\}$ satisfies (35), and $\beta$ and $c$ satisfy the following conditions

Then all conclusions in Theorem 3.1 hold true for the In-Prox-GPDA.

Extension: Increasing the Penalty Sequence

In this section, we present an important variant of Prox-PDA in which there is no need to explicitly compute the bound for the penalty parameter $\beta$. Indeed, the bounds on $\beta$ derived in the previous sections are the worst case bounds, and algorithms that use stepsizes that strictly satisfy such bounds may be slow at the beginning. In practice, one may prefer to start with a small penalty parameter and gradually increase it. The following algorithm adopts such strategy.

We note that in the above algorithm, the primal proximal parameter, the primal penalty parameter, as well as the dual stepsize are all time-varying (in fact all of them increase unboundedly). This is the key feature of this variant. It would be challenging to achieve convergence if only a subset of these parameters grow unboundedly.

Throughout this section we will still assume that Assumption A holds true. Further, we will assume that ${\beta^{r}}$ satisfies the following conditions

Also without loss of generality we will assume that

Note that this is always possible, by adding an identity matrix to $B^{T}B$ if necessary.

The proof of convergence is long and technical, therefore we delegate it to Section 10. The key step is to construct a new potential function, given below

Note that the key to construct the potential function for Algorithm 4 is to make the coefficients in front of the terms $\|x^{r+1}-x^{r}\|^{2}_{B^{T}B}$ and $\|Ax^{r+1}-b\|^{2}$ proportional to $(\beta^{r})^{2}$. Our proof shows that after some finite number of iterations, the newly constructed potential function starts to descend, and the size of the descent is proportional to the following

Therefore if the potential function is lower bounded, then we can conclude that

Using these two inequalities, we can show the desired convergence to stationary solution of problem (1). We refer the readers to to Section 10 for proof details.

We have the following theorem regarding to the convergence of Prox-PDA-IP.

Suppose Assumption A is satisfied. Further suppose that the sequence of penalty parameters $\{\beta^{r}\}$ satisfies (37), and that $B$ is selected such that (38) holds true. Then we have the following claims for Prox-PDA-IP.

(Eventual Feasibility). The constraint is satisfied in the limit, i.e.,

(Convergence to Stationary Points). Every limit point of the iterates $\{x^{r},\mu^{r}\}$ generated by Algorithm 4 converges to a stationary point of problem (1). Further, $Q(x^{r+1},\mu^{r})\to 0$.

We remark the same analysis technique also applies to the following schemes, where $\beta^{r+1}$ also satisfies the condition in (20), and the function $u(x,z)$ is given in (28).

Application: Distributed Nonconvex Optimization

In this section, we discuss the applications of Algorithms 1, 2 and 4 to the nonconvex distributed optimization problem. Our focus will be given to providing explicit update rules for each distributed node, as well as to relating the resulting algorithms to some well-known algorithms in the literature for distributed convex optimization.

We assume throughout this section that each component function $f_{i}$ has Lipschitz continuous gradient with constant $L_{i}$. Then clearly Assumption [A1] is still satisfied, with $L=\max_{i}L_{i}$.

First, we present a direct application of Prox-GPDA (Algorithm 2) to the distributed optimization setting. Quite interestingly, in this case we recover the so-called EXTRA algorithm for distributed convex optimization.

The derivation is in fact rather straightforward. The optimality condition of the $x$-update step (29a) is given by

Utilizing the fact that $A^{T}A=L_{-}$, $B^{T}B=L_{+}$ and $L_{+}+L_{-}=2D$ (where $L_{-}$, $L_{+}$ and $D$ denote respectively the signed Laplacian, the signless Laplacian and the degree matrix), we have

Subtracting the same equation evaluated at the previous iteration, we obtain

where we have used the fact that $A^{T}(\mu^{r}-\mu^{r-1})=\beta A^{T}Ax^{r}=\beta L_{-}x^{r}$. Rearranging terms, we have

where in the last equality we have defined the weight matrix $W:=\frac{1}{2}D^{-1}(L_{+}-L_{-})$, which is a row stochastic matrix.

Iteration (7) has exactly the same form as the EXTRA algorithm given in , therefore we can conclude that the EXTRA is a special case of Prox-GPDA. Moreover, by appealing to our analysis in Section 5, it readily follows that iteration (7) works for the nonconvex distributed optimization problem as well, as long as the parameter $\beta$ is selected appropriately. In particular, in this setting, we need $\beta$ to satisfy

Note that for EXTRA, a sufficient condition for the penalty parameter is that $\beta>L$. Clearly the above requirement for $\beta$ is at least $\frac{L}{2}\left(4+\sqrt{25+\frac{16L^{2}}{\sigma_{\min}(L_{-})}}\right)$ larger. However, this is reasonable since Prox-GPDA is capable of dealing with nonconvex problems as well.

Note that our analysis is by no means any extension to that of EXTRA. Fundamentally, the analysis of EXTRA relies upon showing the descent of the conventional potential function $\|x^{*}-x^{r}\|_{G}$ (i.e., the distance to the optimal solution set), where $G\succeq 0$ being some problem dependent matrix and $x^{*}$ is a solution in the global optimal solution set. In our nonconvex setting, such measure is not useful anymore.

We remark that each node $i$ can distributedly implement iteration (7) by performing the following

where $\mathcal{N}(i)$ denotes the set of neighbors for node $i$

Clearly, at iteration $r+1$, besides the local gradient information, node $i$ only needs the aggregated information from its neighbors, $\sum_{j\in\mathcal{N}(i)}x^{r}_{j}$. Also such aggregated sum is required to be stored in local memory for at least one more iteration, because in order to carried out the $x^{r+1}_{i}$ update, node $i$ also needs $\sum_{j\in\mathcal{N}(i)}x^{r-1}_{j}$.

We also remark that one can apply Prox-PDA (Algorithm 1) to the distributed setting as well. The resulting iteration is given below

This algorithm can be implemented as follows. At iteration $1$, assuming that $x^{0}$ and $x^{-1}$ have been properly initialized. We generate $x^{1}$ according to the following

Equivalently, each node generates $x^{1}_{i}$ according to the following

To obtain such $x^{1}_{i}$, each node solves the following optimization problem

Also, we remark that applying either Algorithm 1 or Algorithm 2, the penalty parameter $\beta$ will still be selected according to (19) and (20). Intuitively, $\beta$ is decreasing with respect to the smallest nonzero eigenvalue of $L_{-}$, increasing with the maximum eigenvalue of $L_{+}$, and finally proportional to $L=\max_{i}L_{i}$.

Finally, in applications where explicitly selecting the stepsizes is difficult, an alternative is to use the Prox-GPDA-IP. To derive the iterates, let us select the $B$ matrix such that $B^{T}B=L_{+}+I_{N}$ (in order to satisfy (38)). Using this choice of $B$, the optimality condition for the $x$-subproblem (40a) is given by

Subtracting the same equation evaluated at the previous iteration, we obtain

where we have used the fact that $A^{T}(\mu^{r}-\mu^{r-1})=\beta^{r}A^{T}Ax^{r}=\beta^{r}L_{-}x^{r}$. Rearranging terms, we have

where in the last equality we have defined the weight matrix $W:=\frac{1}{2}(D+\frac{1}{2}I_{N})^{-1}(L_{+}-L_{-}+I_{N})$, which is a row stochastic matrix.

The discussion in this section is summarized in the following corollary .

Consider the distributed optimization problem (3). Suppose that the graph ({\mbox{\mathcal{V}}},\mathcal{E}) is connected. Then we have the following claims.

Suppose Assumption A is satisfied. Suppose $\beta$ and $c$ satisfy (20) and (19). Then all the conclusions in Theorem 3.1 hold true for the distributed iterations (7) and (43), respectively.

Suppose Assumption A is satisfied, and $\{\beta^{r}\}$ satisfies (37). Then all the conclusions in Theorem 6.1 hold true for the distributed iteration (7).

Moreover, in either case, a stationary solution of the original unconstrained problem (2) is achieved.

The first and the second statements directly follow the results in Theorem 3.1, Corollary 5.1 and Theorem 6.1. The third statement is easy to see because the stationary solution of problem (3) is given by

where the first equality is true because $1$ is in the null space of the incidence matrix $A$. The fact that $Ax^{*}=0$ implies that $x_{i}=x_{j}$, for all $i\neq j$. Therefore we conclude that

This corollary suggests that iterations (7) and (43) also achieve a global sublinear convergence. Note that there has been a few recent works on distributed nonconvex optimization, for example and the references therein. These works design algorithms under different assumptions on the problem as well as on the network structure. However, central to these works is the use of certain diminishing stepsize for th local updates, which results in no global convergence rate guarantees. To the best of our knowledge, the Prox-PDA based distributed algorithms are the first ones that provably achieve global sublinear convergence rate for nonconvex distributed optimization.

Generalization: Distributed Nonconvex Matrix Factorization

In this section we study a variant of the Prox-PDA algorithm for a distributed matrix factorization problem.

Consider the following matrix factorization problem

Consider a distributed scenario where $N$ agents form a graph \{{\mbox{\mathcal{V}}},\mathcal{E}\}, each having a column of $Y$ (note, this can be easily generalized to the case where a subset of columns are available for each agent), we reformulate problem (49) as

Using this condition, we formulate the distributed matrix factorization problem as

Clearly the above problem does not fall into the form of (1), because there are two block variables $\{X_{i}\}$ and $\{y_{i}\}$ in the objective, but the linear constraint only has to do with the $X$-block. Moreover, the objective function couples among the variable blocks $\{X_{i}\}$ and $\{y_{i}\}$ in a nonconvex manner, and neither $\{X_{i}\}$ nor $\{y_{i}\}$ has Lipschitz continuous gradient. The latter fact poses significant difficulty in algorithm development and analysis.

Define the block-signed and the block-signless Laplacians as follows

The augmented Lagrangian for the above problem is given by

Let us consider the following generalization of Algorithm 1 for distributed matrix factorization.

In the above algorithm we have introduced a new sequence $\theta^{r}_{i}\geq 0$, which is some iteration-dependent coefficient representing the size of the local factorization error. We note that including the proximal term $\frac{\theta^{r}_{i}}{2}\|y_{i}-y^{r}_{i}\|^{2}$ is the key to achieving convergence for Algorithm 6. Ideally, the effect of such proximal term will disappear as the algorithm approaches convergence, since we would expect that the local factorization error becomes small. Again one should note that $\frac{\beta}{2}\langle\mathbf{A}\mathbf{X},\mathbf{A}\mathbf{X}\rangle+\frac{\beta}{2}\langle\mathbf{B}(\mathbf{X}-\mathbf{X}^{r}),\mathbf{B}(\mathbf{X}-\mathbf{X}^{r})\rangle$ is strongly convex in $\mathbf{X}$.

We briefly comment on how the algorithm can be implemented in a distributed manner. First note that the $y$ subproblem (52b) is naturally distributed to each node, that is, only local information is needed to perform the update. Second, the $\mathbf{X}$ subproblem (52c) can also be decomposed into $N$ subproblems, one for each node. To be more specific, let us examine the terms in (52c) one by one. First, the term $f(\mathbf{X},Y^{r+1})=\sum_{i=1}^{N}\left(\frac{1}{2}\|X_{i}y^{r+1}_{i}-z_{i}\|^{2}+h_{i}(y_{i})+\gamma\|X_{i}\|^{2}_{F}\right)$, hence it is decomposable. Second, the term \langle{\mbox{\boldmath\Omega}}^{r},\mathbf{A}\mathbf{X}\rangle can be expressed as

where the sets $U(i)$ and $H(i)$ are defined as

Finally, one can verify that the $\Omega$ update step (52d) can be implemented by each edge $e\in\mathcal{E}$ as follows

To show convergence rate of the algorithm, we need the following definition.

In the above derivation, we have defined the proximity operator for a given convex lower semi-continuous function $p(\cdot)$ as

We have used $\mathcal{Y}:=\bigcup_{i}\left\{\|y_{i}\|^{2}\leq\tau\right\}$ to denote the feasible set of $Y$, and used $\iota(\mathcal{Y})$ to denote the indicator function of such set. Similarly as in Section 3, we can show that Q(\mathbf{X}^{r+1},Y^{r+1},{\mbox{\boldmath\Omega}}^{r+1})\to 0 implies that every limit point of (\mathbf{X}^{r+1},Y^{r+1},{\mbox{\boldmath\Omega}}^{r+1}) is a stationary point of problem (50).

Next we present the main convergence analysis for Algorithm 6. The proof is long therefore we delegate it to Section 11.

Consider using Algorithm 6 to solve the distributed matrix factorization problem (50). Suppose that $h(Y)$ is lower bounded over $\mbox{dom}\;h(x)$, and that the penalty parameter $\beta$, together with two positive constants $c$ and $d$, satisfies the following conditions

Then in the limit, consensus will be achieved, i.e.,

Further, the sequences $\{\mathbf{X}^{r+1}\}$ and \{{\mbox{\boldmath\Omega}}^{r+1}\} are both bounded, and every limit point generated by Algorithm 6 converges to a stationary point for problem (49).

Additionally, Algorithm 6 converges sublinearly, i.e., the measure Q(\mathbf{X}^{r+1},Y^{r+1},{\mbox{\boldmath\Omega}}^{r+1}) decreases to in the same manner as in Theorem 3.1. Specifically, for any given $\varphi>0$, define $T$ to be the first time that the optimality gap reaches below $\varphi$, i.e.,

Then there exists a constant $\nu>0$ such that the following is true

We can see that it is always possible to find the tuple $\{\beta,c,d>0\}$ that satisfies (54). For example $c$ can be solely determined by the last inequality; $d$ needs to be chosen large enough such that $1/2-\frac{c}{d}>0$ and $1/2-\frac{c\tau}{d}>0$. After $c$ and $d$ are fixed, one can always choose $\beta$ large enough to satisfy the first three conditions. In practice, we typically prefer to choose $\beta$ as small as possible to improve the convergence speed. Therefore empirically one can start with (for some small $\nu>0$)

and then gradually increase $d$ to find an appropriate $\beta$ that satisfies the first three conditions. Of course, one also has the option of utilizing increasing penalty parameters, just as what we have done in Section 10.

Concluding Remarks

In this paper, we have proposed a decomposition approach for certain linearly constrained nonconvex smooth optimization problem. Our developed algorithms, mainly based upon a novel proximal primal-dual augmented Lagrangian method, are able to decompose the optimization variables, resulting in simple subproblems that can often be solved in closed-form. By constructing a new potential function, which is a conic combination of the augmented Lagrangian, the size of the constraint violation and q certain proximal term, we have shown that the proposed Prox-PDA and its various extensions converge globally sublinearly to the set of stationary solutions. Surprisingly, when specializing a variant of Prox-PDA to the nonconvex distributed optimization problem, the proposed algorithm recovers the popular EXTRA algorithm, indicating that such algorithm converges globally sublinearly even for nonconvex problems.

The proposed proximal primal-dual based algorithm can have many extensions or generalizations. In the paper we have discussed one such generalization to a (distributed) matrix factorization problem. Can we deal with nonsmooth nonconvex terms in the objective, such as indicator functions of convex/nonconvex sets? Can we randomized the algorithm so that each time a randomly selected subproblem is solved instead of the full subproblems? Can we apply our approach to stochastic optimization problems where the objective function involves the expectation of certain nonconvex function? Can we show that some variant of the Prox-PDA converges to local optimal solutions instead of stationary solutions? These are all very interesting research questions that require further investigation.

Proof of Convergence for Algorithm 4

Our analysis consists of a series of steps.

Step 1. Our first step is again to bound the size of the successive difference of $\{\mu^{r}\}$. To this end, write down the optimality condition for the $x$-update (36a) as

Subtracting the previous iteration, we obtain

Also from the optimality condition we have the following relation

where we have defined the primal update direction $v^{r+1}$ as

Step 2. In the second step we analyze the descent of the augmented Lagrangian. We have the following estimate

where in ${\rm(i)}$ we have used the optimality of the $x$-subproblem (cf. the derivation in (3.2)); in ${\rm(ii)}$ we have applied (57).

Step 3. In the third step, we construct the remaining part of the potential function. We have the following two inequalities from the optimality condition of the $x$-update (36a)

Plugging $x=x^{r}$ and $x=x^{r+1}$ to these two equations and adding them together, we obtain

The lhs of the above inequality can be expressed as

Therefore, combining the above three inequalities we obtain

Multiplying both sides by $\beta^{r}$, we obtain

where in the last equality we have merged the terms $\frac{\beta^{r+1}-\beta^{r}}{2\beta^{r}}\|\mu^{r-1}-\mu^{r}\|^{2}$ and $\frac{\beta^{r}(\beta^{r}-\beta^{r-1})}{2}\|Ax^{r}-b\|^{2}$.

Step 4. In this step we construct and estimate the descent of the potential function. For some given $c>0$, let us define the potential function as

Note that this potential function has some major differences compared with the one we used before; cf. (17). In particular, the second and the third terms are now quadratic, rather than linear, in the penalty parameters. This new construction is the key to our following analysis.

Then combining the estimate in (10) and (10), we obtain

where in the inequality we have also used the fact that $\beta^{r}\geq\beta^{r-1}$.

Taking the sum of $r$ from $t$ to $T+1$ (for some $T>t>1$) and utilize again the estimate in (57), we have

First, note that for any $c\in(0,1)$, the coefficient in front of $\left\|(x^{r+1}-x^{r})-(x^{r}-x^{r-1})\right\|_{B^{T}B}^{2}$ becomes negative for sufficiently large (but finite) $t$. This is because $\{\beta^{r}\}\to\infty$, and that the first term in the parenthesis scales in $\mathcal{O}((\beta^{r})^{2})$ while the second term scales in $\mathcal{O}(1)$ . For the first term to be negative, we need $c>0$ to be small enough such that the following is true for large enough $r$

Suppose that $r$ is large enough so that $(\beta^{r+1}-L)/2>{\beta^{r+1}}/{3}$, or equivalently $\beta^{r+1}>3L$. Also choose $c=\min\{1/(4L),1/(12\omega\|B^{T}B\|)\}$, where $\omega$ is given in (37). Then we have

For this given $c$, we can also show that the following is true for sufficiently large $r$

In conclusion we have that for sufficiently large but finite $t_{0}$, we have

Therefore we conclude that if $\{\beta^{r+1}\}$ satisfies (37), and for $c$ sufficiently small, there exits a finite $t_{0}>0$ such that for all $T>t_{0}$, the first two terms of the rhs of (10) is negative.

Step 5. Next we show that the potential function must be lower bounded. Observe that the augmented Lagrangian is given by

where we have used the fact that $\beta^{r+1}\geq\beta^{r}$. Note that $t_{0}$ in (10) is a finite number hence $\frac{1}{2\beta^{t_{0}}}\|\mu^{t_{0}}\|^{2}$ is finite, and utilize Assumption [A2], we conclude that

Combining the above with the fact that the remaining terms of the potential function are all nonnegative, we conclude

Combining (66) and the bound (10) (which is true for a finite $t_{0}>0$), we conclude that the potential function $P_{\beta^{r+1},c}(x^{r+1},x^{r},\mu^{r+1})$ is lower bounded for all $r$.

Step 6. In this step we show that the successive differences of various quantities converge.

The lower boundedness of the potential function combined with the bound (10) (which is true for a finite $t_{0}>0$) implies that

These two facts applied to (10), combined with $\mu^{r+1}-\mu^{r}\in\mbox{col}(A)$, indicate that the following is true

Also (10) implies that the potential function is upper bounded as well, and this indicates that

The second of the above inequality implies that $\beta^{r+1}B^{T}B(x^{r+1}-x^{r})$ is bounded. If we further assume that $\nabla f(x)$ is bounded, and use (55), we can conclude that $\{\mu^{r}\}$ is bounded.

Step 7. Next we show that every limit point of $(x^{r},\mu^{r})$ converges to a stationary solution of problem (1). Let us pass a subsequence $\mathcal{K}$ to $(x^{r},\mu^{r})$ and denote $(x^{*},\mu^{*})$ as its limit point. For notational simplicity, in the following the index $r$ all belongs to the set $\mathcal{K}$.

From relation (67a) we have that any given $\epsilon>0$, there exists $t$ large enough the following is true

Utilizing (58), we have that the following is true

The first relation implies that $\lim\inf_{r\to\infty}\|v^{r+1}\|=0$. Applying these relations to (57), we have

This implies that for any given $\epsilon>0$, $c>0$, there exists an index $t$ sufficiently large such that

Applying this inequality and (71) to (10), we have that for large enough $t$ and for any $T>t$ the following is true

Next we modify a classical argument in [59, Proposition 3.5] to show that

We already know from the first relation in (LABEL:eq:v:summable) that $\lim\inf_{r\to\infty}\|v^{r+1}\|=0$. Suppose that $\|v^{r+1}\|$ does not converge to , then we must have $\lim\sup_{r\to\infty}\|v^{r+1}\|>0$. Hence there exists an $\epsilon>0$ such that $\|v^{r+1}\|<\epsilon/2$ for infinitely many $r$, and $\|v^{r+1}\|>\epsilon$ for infinitely many $\epsilon$. Then there exists an infinite subset of iteration indices $\mathcal{R}$ such that for each $r\in\mathcal{R}$, there exits a $t(r)$ such that

Using the fact that \lim_{r\in{\mbox{\mathcal{K}}}}\mu^{r}=\mu^{*}, we have that for $r$ large enough, the following is true for all $t\geq 0$

Without loss of generality we can assume that this relation holds for all $r\in\mathcal{R}$. Note that the following is true

where in the last inequality we have used (75) and the fact that for all $t\in(r+1,t(r))$, we have $\|v^{t}\|<{\epsilon}$. This implies that

Using the descent of the potential function (74) we have, for $r\in\mathcal{R}$ and $r$ large enough

where in ${\rm(i)}$ we have used the fact that for all $r\in\mathcal{R}$, $\|v^{r+i}\|\geq\frac{\epsilon}{2}$ for $i=1,\cdots,t(r)$; in ${\rm(ii)}$ we have used (77). However we know that the potential function is converging, i.e.,

which contradicts to (10). Therefore we conclude that $\|v^{r+1}\|\to 0$.

Finally, combining $\|v^{r+1}\|\to 0$ with the convergence of $\mu^{r+1}-\mu^{r}$ (cf. (69)), we conclude that every limit point of $\{x^{r},\mu^{r}\}$ satisfies

Therefore it is a stationary solution for problem (1). This completes the proof.

Proof of Convergence for Algorithm 6

To make the derivation compact, define the following matrix

Step 1. First we note that the optimality condition for the $\mathbf{X}$-subproblem (52c) is given by

By utilizing the fact that that {\mbox{\boldmath\Omega}}^{r+1}-{\mbox{\boldmath\Omega}}^{r} lies in the column space of $\mathbf{A}$, and the eigenvalues of $A^{T}A$ equal to the eigenvalue of $\mathbf{A}^{T}\mathbf{A}$, we have the following bound

Next let us analyze the first term in the rhs of the above inequality. The following identity holds true

where in the last inequality we have defined the constant $\theta^{r}_{i}$ as

Therefore, combining the above two inequalities, we obtain

Step 2. Next let us analyze the descent of the augmented Lagrangian. First we have

where in the second to the last equality we have used the convexity of $h_{i}$, and $\zeta^{r+1}_{i}\in\partial h_{i}(y_{i}^{r+1})$; the last inequality uses the optimality condition of the $y$-step (52b). Similarly, we can show that

where we have utilized the fact that $\mathbf{A}^{T}\mathbf{A}+\mathbf{B}^{T}\mathbf{B}=2\mathbf{D}\succeq{\bf I}_{NM}.$ Therefore, combining the estimate (11), we obtain

Step 3. This step follows Lemma 3.3 in the analysis of Algorithm 1. In particular, after writing down the optimality condition of the $\mathbf{X}^{r+1}$ and $\mathbf{X}^{r}$ step, we can obtain

Then it is easy to show that the above inequality implies the following (which utilizes the convexity of $h$)

where in ${\rm(i)}$ we utilize the convexity of $f(\mathbf{X},y)$ wrt $\mathbf{X}$ for any fixed $y$; in ${\rm(ii)}$ we use the Cauchy-Swarch inequality, where $d>0$ is a constant (to be determined later); ${\rm(iii)}$ is true due to a similar calculation as in (11).

Step 4. Let us define the potential function as

Then utilize the bound in (11) in Step 2 and bounds (11) in Step 3, we obtain

Therefore the following are the condition that guarantees the descent of the potential function

To see that it is always possible to find the tuple $(\beta,c,d)$, first let us set $c$ such that the last inequality is satisfied

Second, let us pick any $d$ such that the following is true

Then clearly it is possible to make $\beta$ large enough such that all the four conditions in (90) are satisfied.

Step 5. We need to prove that the potential function is lower bounded. We lower bound the augmented Lagrangian as follows

Then by the same argument leading to (22), we conclude that as long as $h_{i}$ is lower bounded over its domain, then the potential function will be lower bounded.

Step 6. Combining the result in Step 5 and Step 4, we conclude the following

That is, in the limit the network-wide consensus is achieved. Next we show that the primal and dual iterates are bounded.

Note that the potential function is both lower and upper bounded. Combined with (93) we must have that the augmented Lagrangian is both upper and lower bounded. Using the expression (11), the assumption that $h_{i}(y_{i})$ is lower bounded, and the fact that $y_{i}$ is bounded, we have that in the limit, the following term is bounded

This implies that the primal variable sequence $\{X^{r+1}_{i}\}$ are bounded for all $i$. To show the boundedness of the dual sequence, note that {\mbox{\boldmath\Omega}}^{r+1}\in\mbox{col}(\mathbf{A}) (due to the initialization that {\mbox{\boldmath\Omega}}^{0}=\mathbf{0}). Therefore using (80) we have

Note that from the expression of $\mathbf{M}$ in (11), we see that $\{\mathbf{M}^{r+1}\}$ is bounded because both $\mathbf{X}^{r+1}$ and $Y^{r+1}$ are bounded. Similarly, the second term on the rhs of the above inequality is bounded because $\mathbf{X}^{r+1}\to\mathbf{X}^{r}$. These two facts imply that \{{\mbox{\boldmath\Omega}}^{r+1}\} is bounded as well.

Arguing the convergence to stationary point as well as the convergence rate follows exact the same steps as in the proof of Theorem 3.1.

References