Solving Multiple-Block Separable Convex Minimization Problems Using Two-Block Alternating Direction Method of Multipliers

Introduction

In this paper, we consider the following convex optimization problem with $m$ block variables and the objective being the sum of $m\ (m\geq 2)$ separable convex functions:

while $\|\cdot\|_{*}$, $\|\cdot\|_{1}$ and $\|\cdot\|_{F}$ denote respectively the matrix nuclear norm (i.e., the sum of the matrix singular eigenvalues), the $L_{1}$ and Frobenius norm of a matrix. Clearly problem (1.2) corresponds to the case of $m=3$ in problem (1), with $x_{1}=L$, $x_{2}=S$, $x_{3}=Z$ and

where $\mathcal{I}_{\mathcal{H}}(\cdot)$ denotes the indicator function for the set ${\mathcal{H}}$.

Similarly, the so-called latent variable Gaussian graphical model selection (LVGGMS) problem , which is closely related to the inverse covariance matrix estimation problem, is also in the form of (1). In particular, suppose $(X,Y)$ is a pair of $(p+r)$-dimensional joint multivariate Gaussian random variables, with covariance matrix denoted by $\Sigma_{(X,Y)}:=[\Sigma_{X},\Sigma_{XY};\Sigma_{YX},\Sigma_{Y}]$ and its inverse $\Theta_{(X,Y)}:=[\Theta_{X},\Theta_{XY};\Theta_{YX},\Theta_{Y}]$ respectively. The random variable $X:=(X_{1},X_{2},\cdots,X_{p})^{T}$ is observable while $Y:=(Y_{1},Y_{2},\cdots,Y_{r})^{T}$ is the latent (or hidden) random variable. In many applications, we typically have $r\ll p$. Moreover, the marginal distribution of the observed variables $X$ usually follows a sparse graphical model and hence its concentration matrix $\Theta_{X}$ is sparse. Notice that the inverse of the covariance matrix for $X$ can be expressed as

which is the difference between the sparse term $\Theta_{X}$ and the low-rank term $\Theta_{XY}\Theta_{Y}^{-1}\Theta_{YX}$ (since $r$ is much less than $p$). Thus, the task of estimating the sparse marginal concentration matrix $\Theta_{X}$ can be accomplished by solving the following regularized maximum likelihood problem

where the constraint $R\succ 0$ is implicitly imposed by having the term $-\hbox{logdet}(R)$ in the objective function. It is easily seen that LVGGMS corresponds to the three block case $m=3$ in problem (1) with $x=(R,S,L)$ and

where the linear constraint is $R-S+L=0$.

Problem (1) is a structured convex problem with a separable objective function and a single linear equality constraint. A popular algorithm for solving this class of problem is the so-called Alternating Direction Method of Multipliers (ADMM). To outline the basic steps of the ADMM, we first introduce the augmented Lagrangian function for problem (1)

where $\beta$ is the penalty parameter for the violation of the linear constraint and $\langle\cdot,\cdot\rangle$ denotes the standard trace inner product. The ADMM method is a Gauss-Seidel iteration scheme in which the primal block variables $\{x_{i}\}$ and the Lagrangian multiplier $\lambda$ for the equality constraint are updated sequentially. Specifically, for a fixed stepsize $\beta>0$, the ADMM for solving problem (1) can be described as follows:

The history of ADMM dates back to 1970s in where the method was first developed for solving 2-block separable convex problems. In , it is shown that ADMM can be interpreted as a special case of an operator splitting method called Douglas-Rachford Splitting Method (DRSM) for finding a zero of the sum of two monotone operators $\mathcal{A}$ and $\mathcal{B}$ . Moreover, ADMM can also be related to Spingarn’s method called Partial Inverse . Recently, ADMM has found its application in solving a wide range of large-scale problems from statistics, compressive sensing, image and video restoration, and machine learning, see e.g., and the references therein.

The ADMM convergence has long been established in the literature for the case of two block variables (i.e. $m=2$ in problem (1)). References show that the algorithm converges globally when each subproblem is solved exactly. Such convergence results have also been obtained in the context of DRSM. In , the authors show that DRSM is a special case of the so-called Proximal Point Algorithm (PPA), for which the global convergence and the rate of convergence have been established (see ). Accordingly, under certain regularity conditions, the global convergence and the rate of convergence of DRSM (and hence ADMM) follow directly. On the other hand, for large scale problems such as those arising from compressive sensing, the global optimal solution for subproblems related to certain block variables may not be easily computable . In these cases the classical ADMM needs to be modified accordingly so that those difficult subproblems are solved inexactly. In , the authors show that by performing a simple proximal gradient step for each subproblem, global convergence results similar to those for the classical ADMM can also be obtained. Beyond global convergence, there are a few results characterizing the iteration complexity and the convergence rate of the ADMM. The authors of have shown that to obtain an $\epsilon$-optimal solution, the worst-case iteration complexity of both exact and inexact ADMM is $\mathcal{O}(1/\epsilon)$, where the $\epsilon$ optimality is defined using both the constraint and objective violation. In , Rockafellar has shown that if the inverse of the considered operator is Lipschitz continuous at origin point, the PPA converges linearly when the resolvent operator is solved either exactly or inexactly. Therefore the linear convergence of DRSM and ADMM follow directly under some assumptions on $\mathcal{A}$ and $\mathcal{B}$ in DRSM or $\{f_{i}\}$ and $\{A_{i}\}$ in ADMM. Further in , Lions and Mercier have proved that when operator $\mathcal{B}$ is both coercive and Lipschitz, then DRSM converges linearly. Further in , Eckstein and Bertsekas have shown the linear convergence rate of ADMM for linear programming. More recently the authors of show that for both exact version and inexact version involving a proximal term, ADMM converges linearly if the objective function is strongly convex and its gradient is Lipschitz continuous in at least one block variable, and that the matrices $\{A_{i}\}$ satisfy certain rank assumptions.

However, when the number of block variables is greater than two, the convergence of ADMM has not been well understood. Recently, the authors of prove the global (linear) convergence of the ADMM for the multiple-block problem (1), under the following assumptions: a) for each $i$, $A_{i}$ is full column rank, b) $\beta$ is sufficiently small and c) certain error bounds hold for the problem (1). The full column rank condition a) can be dropped if each subproblem is solved inexactly. However when conditions b) and c) are not satisfied, even the global convergence of the algorithm is still open in the literature. As a result, a number of variants of the classical ADMM have been proposed; see . For example, in , by adding an additional Gaussian back substitution correction step in each iteration after all the block variables are updated, the authors establish the global convergence, iteration complexity and linear convergence rate (under some technical assumption on the iterates) of the modified algorithm. However such correction step is not always computationally efficient, especially when $A_{i}$’s are not efficiently invertible. In , an alternating proximal gradient method is proposed, in which the proximal version of the ADMM is applied to problem (1), after grouping the $m$ block variables into two blocks. However, the way that the block variables should be grouped is highly problem dependent. There is no general criterion guiding how such step should be done. Also recently, some multiple block splitting algorithms have been proposed for solving some models similar to (1) in .

In this work, we systematically study ADMM algorithms for solving the multi-block problem (1). We first propose two novel algorithms that apply the two-block ADMM to certain reformulation of the original multi-block problem. We show in detail how these algorithms are derived and analyze their convergence properties. We then report numerical results comparing the original multi-block ADMM with the proposed approaches on problems with multiple block structures such as the basis pursuit problem, robust PCA and latent variable Gaussian graphical model selection. Our numerical experiments show that the multi-block ADMM performs much better than the two-block ADMM algorithms as well as many other existing algorithms.

where $\tau>0$ and $c$ are given. For instance when $g(x)=\|x\|_{1}$,

Let ${\mathcal{T}}$ be a set-valued operator, $c$ is any positive scalar, and denote the following operator

The rest of this paper is organized as follows: the primal and dual splitting ADMM algorithms are presented and analyzed in Section 2.1 and Section 2.2, respectively. In Section 3, numerical results for both synthetic and real problems are reported. Finally, some concluding remarks are given in Section 4.

The Proposed Algorithms

Various modifications of the ADMM algorithm have been proposed in the literature to deal with the multi-block problem (1). Instead of proposing yet another ADMM variant, we propose to transform any multi-block problem into an equivalent two-block problem to which the classical two-block ADMM can be readily applied. The main idea is to appropriately introduce some auxiliary variables so that the original variables are completely decoupled from each other. In this section, this technique will be explained in details for both the primal and dual versions of problem (1).

Problems (1) and (2.1) are equivalent in the sense that they share the same primal optimal solution set for the variables $\{x_{i}\}_{i=1}^{m}$, and achieve the same global optimal objective value. To apply the ADMM algorithm to the above reformulated problem, we write its (partial) augmented Lagrangian function, by keeping the constraint $\sum\limits_{i=1}^{m}y_{i}=0$ and penalizing the rest of the constraints:

We note that the subproblem for $y$ is a projection onto the hyperplane $\sum_{i}y_{i}=0$. As such, it admits the following closed-form solution

Further, it is easy to see that the subproblem for $x_{i}$ can be solved efficiently by using (1.7), provided that $A_{i}$ is an identity matrix (or any constant multiple of it). Else, $x_{i}$ can be updated by simply using a proximal gradient step:

where $\tau_{i}$ denotes the penalty parameter for the distance between $x_{i}^{k+1}$ and $x_{i}^{k}$.

Next we discuss the convergence of the above primal splitting ADMM algorithm.

For any $\beta>0$, suppose the subproblem for $x_{i}$ is either exactly solved, or is solved inexactly using (2.4) with $\tau_{i}>\beta\cdot\rho(A_{i}^{T}A_{i})$. Let (\mbox{\boldmathx}^{k},\mbox{\boldmathy}^{k},\mbox{\boldmath\lambda}^{k}) be any sequence generated by Algorithm 2. Then starting with any initial point (\mbox{\boldmathx}^{0},\mbox{\boldmathy}^{0},\mbox{\boldmath\lambda}^{0})\in{\mathcal{X}}\times{\mathcal{Y}}\times\Lambda, we have

The sequence \{\mbox{\boldmath\lambda}^{k}\} converges to \mbox{\boldmath\lambda}^{*}, where \mbox{\boldmath\lambda}^{*} is the dual optimal solution for problem (2.1);

The sequence $\{\sum\limits_{i=1}^{m}f_{i}(x_{i}^{k})\}$ converges to $p^{*}$, where $p^{*}$ is the primal optimal value for problem (2.1);

The residual sequence $\{A_{i}x_{i}^{k}-y^{k}_{i}-\frac{b}{m}\}$ converges to 0 for each $i=1,\cdots,m$;

If the subproblem for $x_{i}$ is exactly solved for each $i$, then the sequence $\{A_{i}x_{i}^{k}\}$ and $\{y^{k}_{i}\}$ converge. Moreover, if $A_{i}$ has full column rank, then $\{x_{i}^{k}\}$ converges. When the subproblem for $x_{i}$ is solved inexactly using (2.4) with $\tau_{i}>\beta\cdot\rho(A_{i}^{T}A_{i})$, the sequence \{({\mbox{\boldmathx}}^{k},{\mbox{\boldmathy}}^{k})\} converges to an optimal solution to problem (2.1).

Proof. When the subproblems are solved exactly, we are actually using the classical two-block ADMM to solve the equivalent formulation (2.1). As a result, the first three conclusions as well as the convergence of $\{A_{i}x_{i}^{k}\}$ and $\{y^{k}_{i}\}$ follow directly from the classical analysis of the two-block ADMM (see, e.g., [2, Section 3.2]). The convergence of $\{x^{k}_{i}\}$ is a straightforward consequence of the convergence of $\{A_{i}x_{i}^{k}\}$ and the assumption that $A_{i}$’s are all full column rank. When the subproblems are solved inexactly via (2.4), because $\tau_{i}>\beta\cdot\rho(A_{i}^{T}A_{i})$, all the conclusions are implied by the result in [23, Theorem 1]. $\square$

Besides global convergence, we can elaborate on other convergence properties of Algorithm 2. First, for both the exact case and the inexact proximal case with $\tau_{i}>\beta\cdot\rho(A_{i}^{T}A_{i})$, we can obtain the following iteration complexity result by adopting the variational inequality framework developed in [27, Theorem 4.1]. To illustrate, let \mbox{\boldmathw}:=(\mbox{\boldmathy}^{T},\mbox{\boldmathx}^{T},\mbox{\boldmath\lambda}^{T})^{T}\in{\mathcal{Y}}\times{\mathcal{X}}\times\Lambda, and let \{\mbox{\boldmathw}^{k}\} denote the sequence generated by Algorithm 2. Further we define

where C_{p}^{1}=\max_{{\mbox{\boldmathw}}\in{\mathcal{W}}}\|\mbox{\boldmathw}-\mbox{\boldmathw}^{0}\|_{H}^{2} and $H$ is a positive semi-definite matrix which is associated with $\{A_{i}\}$, $\beta$ and $\{\tau_{i}\}$. Note that at optimality, the left hand side of (2.11) is no greater than zero, therefore the above inequality is indeed a possible measure of the optimality gap, although it is implicit. In the following, we show explicitly that the objective values decrease at the rate ${\mathcal{O}}(1/K)$.

Let \{{\mbox{\boldmathx}}^{k},{\mbox{\boldmathy}}^{k},{\mbox{\boldmath\lambda}}^{k}\} be the sequence generated by Algorithm 2, and ({\mbox{\boldmathx}}^{*},{\mbox{\boldmathy}}^{*},{\mbox{\boldmath\lambda}}^{*}) be any optimal solution, then we have

It is worth noting that the complexity results presented in Theorem 2.2 and the one presented in (2.11) do not imply each other. Moreover, in the proof of Theorem 2.2, we have explored certain structure of Algorithm 2, therefore this result does not carry over to the general ADMM algorithm.

Next we show that the linear rate of convergence for Algorithm 2 can also be established using existing results for the Dauglas-Rachford Splitting Method (DRSM). To this end, we first derive the relationship between Algorithm 2 and the DRSM. Recall that DRSM solves the following problem

by generating two sequences $\{u^{k}\}$ and $\{v^{k}\}$ according to:

To see the exact form of the operators ${\mathcal{A}}$ and ${\mathcal{B}}$ for Algorithm 2, let us consider the dual formulation of (2.1), stated below

where ${\mathcal{I}}(\cdot)$ denotes the indicator function. By setting

we can rewrite the dual form of (2.1) (i.e., eq. (2.13)) equivalently as finding a {\mbox{\boldmath\lambda}}^{*} that satisfies

Applying DRSM to solve (2.16), we obtain the ($k+1$)th iterate as follows

Further, applying the Fenchel-Rockafellar Duality [1, Definition 15.19], we can write the dual problem for (2.17) and (2.18) (with dual variables $y$ and $x$) as

Then obviously when substituting $u_{i}^{k}=v_{i}^{k}-\beta Ax_{i}^{k}$ into the subproblem about $y$, we obtain

Similarly, when substituting $v_{i}^{k+1}=u_{i}^{k}-\beta\left(-y_{i}^{k+1}-\frac{b}{m}\right)$ into the subproblem about $x$, we can get the following equivalent problem for $x$

Combining $u_{i}^{k+1}=v_{i}^{k+1}-\beta Ax_{i}^{k+1}$ and $v_{i}^{k+1}=u_{i}^{k}-\beta\left(-y_{i}^{k+1}-\frac{b}{m}\right)$, we obtain the update of $u$,

The above analysis indicates that the sequence \{{\mbox{\boldmathu}}^{k}\} is the same as the multiplier sequence \{{\mbox{\boldmath\lambda}}^{k}\} in Algorithm 2. As a result, Algorithm 2 (or in general the two-block ADMM) can be considered as a special case of DRSM.

The linear convergence of DRSM has been well studied in [29, Proposition 4] with an assumption that operator ${\mathcal{B}}$ is both strongly monotone and Lipschitz, which means there exists $\alpha>0$ and $M$ such that

By using the results in , we can show that if $f_{i}$’s are all strongly convex with Lipschitz continuous gradients, and when $A_{i}$’s are all full row rank, then the operator ${\mathcal{B}}$ is strongly monotone and Lipschitz. As a result, the sequences \{{\mbox{\boldmathx}}^{k}\}, \{{\mbox{\boldmathy}}^{k}\} and \{{\mbox{\boldmath\lambda}}^{k}\} generated by Algorithm 2 converge linearly. Similarly, if each subproblem cannot be solved exactly, then the linear convergence of the inexact version of Algorithm 2 (cf. (2.4), with $\tau_{i}>\beta\cdot\rho(A_{i}^{T}A_{i})$) can be established by following [10, Theorem 4]. Again we require that $f_{i}$’s are all strongly convex with Lipschitz continuous gradients, and $A_{i}$’s all have full row rank.

To this point, all the convergence results characterize the behavior of Algorithm 2 for solving problem (2.1). As problem (1) is an equivalent reformulation of (2.1), we can readily conclude that the sequence \{\mbox{\boldmathx}^{k}\} generated by Algorithm 2 converges to the primal optimal solution of (1), if either each subproblem is exactly solved and $A_{i}$’s are all full rank, or the subproblems are solved using the proximal step (2.4) with $\tau_{i}>\beta\cdot\rho(A_{i}^{T}A_{i})$. Further, if (\mbox{\boldmathx}^{k},{\mbox{\boldmathy}}^{k}) is linearly convergent to an optimal solution of problem (2.1), then \mbox{\boldmathx}^{k} converges linearly to the primal optimal solution of (1).

2. Dual Splitting ADMM

We can also apply the splitting technique to the dual formulation (1) to derive a dual splitting ADMM algorithm. In particular, let us first write (1) in its saddle point form

By exchanging the order of $\max$ and $\min$, and using the definition of the conjugate function of $f_{i}$, we can rewrite (2.1) equivalently as

where $\lambda$ denotes the dual variable of (1). We then split the dual variable $\lambda$ by introducing a set of auxiliary variables $\{\lambda_{i}\}_{i=1}^{m}$, and rewrite (2.2) as

It is obvious that each primal optimal solution $\{\lambda_{i}^{*}\}_{i=1}^{m}$, $\lambda^{*}$ of (2.3) corresponds to a dual optimal solution of (1) ($\lambda_{i}^{*}=\lambda^{*}$). The augmented Lagrangian function for this dual problem can be expressed as follows:

By the classical Fenchel-Rockafellar duality , the dual problem of (2.5) can be expressed as

where $\{x_{i}\}$ is precisely the set of primal variables of (1). The relationship between $\lambda_{i}^{k+1}$ and $x_{i}^{k+1}$ is as follows

The dual splitting ADMM is stated formally in the following table.

Similar to the case of primal splitting, the subproblem of $\lambda$ can be solved easily in closed-form

Furthermore, the subproblem for the block variable $x_{i}$ can be solved efficiently if $A_{i}$ is an identity matrix (or any of its constant multiples), because (2.6) can be efficiently solved by computing the proximity operator (1.7). Else, a proximal gradient step can be performed, i.e.,

where $\tau_{i}$ denotes the penalty parameter for the distance between $x_{i}^{k+1}$ and $x_{i}^{k}$.

Again, the global convergence of Algorithm 3 is a straightforward consequence of the standard convergence results for the two-block ADMM.

The sequence \{\mbox{\boldmatht}^{k}\} converges to the dual optimal solution for problem (2.3).

The sequence $\left\{\sum\limits_{i=1}^{m}f_{i}^{*}(A_{i}^{T}\lambda_{i}^{k})-\langle\lambda^{k},b\rangle\right\}$ converges to the primal optimal value for problem (2.3).

The residual sequence $\{\lambda^{k}-\lambda_{i}^{k}\}$ converges to 0 for each $i=1,\cdots,m$.

For each $i=1,\cdots,m$, if the subproblem about $x_{i}$ is exactly solved, then the sequence $\{A_{i}x_{i}^{k}\}$ converges. If $A_{i}$ has full column rank, then $\{x_{i}^{k}\}$ converges to $x_{i}^{*}$, for all $i=1,\cdots,m$; the same is true if the subproblem about $x_{i}$ is solved using (2.9) with $\tau_{i}>\frac{\rho(A_{i}^{T}A_{i})}{\beta}$.

Proof. When the subproblems are solved exactly, Algorithm 3 corresponds to the classical two-block ADMM applied to solve the equivalent formulation (2.3). As a result, the first four conclusions follow directly from the classical analysis of the two-block ADMM (see, e.g., [2, Section 3.2]). In the last conclusion, the convergence of $\{A_{i}x_{i}^{k}\}$ follows from the convergence of $\{\lambda^{k}_{i}\}$, $\lambda^{t}$ and $\{t^{t}_{i}\}$; see (2.7). When the subproblems are solved inexactly with $\tau_{i}>\frac{\rho(A_{i}^{T}A_{i})}{\beta}$, there is no existing result which covers the convergence of the algorithm. We will provide a proof for this case in the Appendix. $\square$

Let us discuss some additional convergence properties of Algorithm 3. First of all, it is possible to derive the iteration complexity for both the exact and the inexact versions of the dual splitting ADMM algorithm. For the exact version, its iteration complexity based on variational inequalities follows from the existing results . The iteration complexity of the inexact dual splitting ADMM algorithm is not covered by any existing result. As a result, in Appendix, we provide a unified iteration complexity analysis for the dual splitting ADMM algorithm. Additionally, similar to Theorem 2.2, we have the following result that bounds the gap of objective value. The proof is similar to that of Theorem 2.2, thus we omit it for brevity.

From the equivalence relationship between the problems (1) and (2.3), we can readily claim that the primal optimal solution $\lambda^{*}$ of (2.3) is the dual optimal solution of (1). Recall that from the discussion following (2.6), $\{x_{i}\}$ is the set of primal variables for the original problem (1).

3. Discussions

Several existing methods for solving (1) are similar to the two algorithms (Algorithm 2 and Algorithm 3) proposed in this paper. In particular, Spingarn applied a method called partial inverse method to separable convex problems. This method can be directly applied to (1) as follows. Let us define two subspaces $A$ and $B$ as:

Then problem (1) is equivalent to the one that minimizes F({\mbox{\boldmathx}},{\mbox{\boldmathu}}) over $A$. Define two operators

To solve (1), partial inverse method generates a sequence of iterates \{({\mbox{\boldmathx}}^{k},{\mbox{\boldmathy}}^{k})\} and \{(0,{\mbox{\boldmath\lambda}}^{k})\}:

with positive sequence $\{c_{k}\}$ bounded away from zero.

From [38, Algorithm 2], it is known that the partial inverse method is the same as ADMM for minimizing F({\mbox{\boldmathx}},{\mbox{\boldmathy}}) over $A$. That is, the variables ({\mbox{\boldmathx}},{\mbox{\boldmathy}}) in the Partial Inverse are the same as those in Algorithm 2 when the subproblems about $x_{i}$ are solved exactly. However, the subproblem about $\lambda$ in the partial inverse method additionally requires every component of $\lambda$ to be equal, which is different from that in Algorithm 2. Notice that the results in do not apply to the case when the subproblems for $x_{i}$ are solved inexactly.

Algorithm 3 is related to the proximal decomposition method proposed in . The latter solves

by applying the two block ADMM to the following reformulation

where ${\mathcal{X}}$ is the common closed convex constrained set for all $\{x_{i}\}_{i=1}^{m}$. In this decomposition, a single variable $x$ is split into $m$ copies $\{x_{i}\}_{i=1}^{m}$, and the consistency among these copies are enforced using the linking variable $y$. This decomposition technique is also used in Algorithm 3, but for solving the dual reformulation of (1).

Algorithm 2 is closely related to the distributed sharing algorithm presented in [2, Chapter 7]. Consider the following sharing problem, in which $m$ agents jointly solve the following problem

where $f_{i}(\cdot)$ is the cost related to agent $i$; $g(\cdot)$ is the cost shared among all the agents. The distributed sharing algorithm introduces a set of extra variables $y_{i}=A_{i}x_{i}-\frac{b}{m},\ \forall\ i$, and applies the two-block ADMM to the following reformulation

To see the relationship between Algorithm 2 and the distributed sharing algorithm, we note that problem (1) is a special case of problem (2.14), with $g(\cdot)$ being the indicator function. Hence Algorithm 2 with the subproblems being solved exactly can be viewed as a special case of the distributed sharing algorithm.

Numerical Experiments

In this section, we test Algorithms 1, 2 and 3 on three problems: Basis Pursuit, Latent Variable Gaussian Graphical Model Selection and Robust Principal Component Analysis, and compare their performance with the Alternating Direction Method (ADM), Proximal Gradient based ADM (PGADM) , ADMM with Gaussian Back Substitution (ADMGBS) , and Variant Alternating Splitting Augmented Lagrangian Mehtod (VASALM) . Our codes were written in Matlab 7.14(R2012a) and all experiments were conducted on a laptop with Intel Core 2 Duo@2.40GHz CPU and 4GB of memory.

Consider the following basis pursuit (BP) problem

By letting {\mbox{\boldmathx}}=({\mbox{\boldmathx}}_{1},\cdots,{\mbox{\boldmathx}}_{m}), the BP problem can be viewed as a special case of problem (1) with $m$ block variables. If we set $m=p$ (i.e., each component $x_{i}$ is viewed as a single block variable), then Algorithm 1 can be used, with each of its primal iteration given by

Alternatively, when the number of blocks $m$ is chosen as $m<p$, the primal subproblem in Algorithm 1 cannot be solved exactly. In this case, the inexact version of Algorithm 2 and Algorithm 3 can be used, where the primal subproblems (2.4) and (2.9) are respectively given by

In the following, we compare Algorithm 1 (with $m=p$), and the inexact versions of Algorithm 2, Algorithm 3 (with $m=2,5,10,20,50,100,200$) with the ADM algorithm , which has been shown to be effective for solving BP. Algorithms 1–3 are denoted as MULTADMM, PSADMM and DSADMM, respectively.

In our experiment, the matrix $A$ is randomly generated using standard Gaussian distribution per element; the true solution $x^{*}$ is also generated using standard Gaussian distribution, with $6\%$ sparsity level, i.e., $94\%$ of the components are zero; the “observation” vector $b$ is computed by $b=Ax^{*}$. The stepsize $\beta$ is set to be $\frac{400}{\|b\|_{1}}$, $\frac{400}{\|b\|_{1}}$, 10 and $\frac{400}{\|b\|_{1}}$ for MULTADMM, PSADMM, DSADMM and ADM respectively. Note that the stepsize for the DSADMM is chosen differently because it is the ADMM applied to the dual of (1), while the rest of the algorithms applies directly to the primal version of (1). To ensure convergence, the proximal parameters are set to be $\tau_{i}=1.01\beta\times\rho(A_{i}^{T}A_{i})$ for PSADMM, $\tau_{i}=1.01\times\frac{\rho(A_{i}^{T}A_{i})}{\beta}$ for DSADMM, and $\tau=1.01\beta\times\rho(A^{T}A)$ for ADM, respectively.

Figure 1 shows the convergence progress of the different algorithms. The curves in the figure represent the relative error for different algorithms along the iterations averaged over 100 runs. For a given iterate $x^{k}$, the relative error is defined as

The left part of Figure 1 shows the performance of MULTADMM, PSADMM and ADM. We observe that PSADMM converges faster than ADM when the number of blocks is relatively small. When the number of blocks increases, PSADMM converges fast at the beginning, but becomes slower after about $200$ iterations. This is because larger number of blocks results in smaller proximal parameter $\tau_{i}$, hence larger stepsizes can be taken Note that the 2-norm of a matrix (i.e., its largest singular value) cannot decrease if some of its columns are removed.. On the other hand, it becomes increasingly difficult to simultaneously satisfy all the constraints. (Note that the number of constraint is the same as the number of block variables). We also observe that the MULTADMM performs much better than all other methods, in terms of both the convergence speed and the solution accuracy. The main reason for its superior performance is the exact solvability of the primal subproblems. Similar observations can be obtained from the right part of Figure 1, where the performance of MULTADMM, ADM and DSADMM are compared.

In Table 1 and Table 2, we report the performance of different algorithms for cases with $n=300,\ p=1000$ and $n=600,\ p=2000$, respectively. In these tables, ‘Iter’ denotes the iteration number with $2000$ as the default maximum iteration number; ‘Obj’ denotes the objective value; ‘Time’ denotes the CPU time used; ‘$\sim$’ indicates that the algorithm did not converge in 2000 iterations. Again we observe that to obtain the same accuracy, MULTADMM requires significantly fewer iterations and less CPU time than all other methods. In the meantime, PSADMM and DSADMM perform better than ADM when the number of blocks is not too large.

2. Latent Variable Gaussian Graphical Model Selection

The problem of latent variable Gaussian graphical model selection has been briefly introduced in Section 1. Recall that one of its equivalent reformulation is given by

This model can be viewed as a combination of dimensionality reduction (to identify latent variables) and graphical modeling (to capture remaining statistical structure that is not attributable to the latent variables). It consistently estimates both the number of hidden components and the conditional graphical model structure among the observed variables. In the following, we show that to solve (3.4), the primal subproblems for Algorithm 1, Algorithm 2 and Algorithm 3 can be solved exactly and efficiently. To this end, the following two lemmas which can be found in are needed.

Now we are ready to present the steps of different algorithms for solving (1.5), by using the previous two lemmas.

Algorithm 1: At the $k$th iteration, the update rule is given by:

where $U\hbox{diag}(\sigma)U^{T}$ is the eigenvalue decomposition of matrix $\frac{1}{\beta}\hat{\Sigma}_{X}\!-\!\frac{1}{\beta}\lambda^{k}\!-\!S^{k}\!+\!L^{k}$ and $\gamma_{i}=\left(-\sigma_{i}+\sqrt{\sigma_{i}^{2}+\frac{4}{\beta}}\right)/2,\ \forall\;i=1,\cdots,p.$

Algorithm 2: At the $k$th iteration, the update rule is given by:

where $U\hbox{diag}(\sigma)U^{T}$ is the eigenvalue decomposition of another matrix $\frac{1}{\beta}\hat{\Sigma}_{X}\!-\!(y_{1}^{k+1}\!+\!\frac{1}{\beta}\lambda_{1}^{k})$ and $\gamma_{i}=\left(-\sigma_{i}+\sqrt{\sigma_{i}^{2}+\frac{4}{\beta}}\right)/2,\ \forall\;i=1,\cdots,p.$

Algorithm 3: At the $k$th iteration, the update rule is given by:

where $U\hbox{diag}(\sigma)U^{T}$ is the eigenvalue decomposition of another new matrix $\frac{1}{\beta}\hat{\Sigma}_{X}\!-\!(\beta\lambda^{k+1}\!-\!t_{1}^{k})$ and $\gamma_{i}=\left(-\sigma_{i}+\sqrt{\sigma_{i}^{2}+\frac{4}{\beta}}\right)/2,\ \forall\;i=1,\cdots,p.$

In the following, all three methods are compared with PGADM , which is used to solve the same problem. The stopping criterion is set to be

where $\epsilon$ is some given error tolerance. All the variables are initialized as zero matrices, and the error tolerance $\epsilon$ is set to be $10^{-5}$ in all the experiments. The comparison results are presented in Table 3, in which ‘Iter’, ‘Obj’ and ‘Time’ denote respectively the iteration number, the objective function value and the CPU time.

while its low rank part is computed as $\Theta_{XY}\Theta_{Y}^{-1}\Theta_{YX}=\Theta_{X,Y}(1:p,p+1:p+r)\cdot\Theta_{X,Y}(p+1:p+r,p+1:p+r)^{-1}\cdot\Theta_{X,Y}(p+1:p+r,1:p)$.

We then draw $N=5p$ independent and identically distributed vectors, $Y_{1},\cdots,Y_{N}$, from the Gaussian distribution ${\mathcal{N}}(0,(\Theta_{X}-\Theta_{XY}\Theta_{Y}^{-1}\Theta_{YX})^{-1})$, and compute a sample covariance matrix of the observed variables according to $\hat{\Sigma}_{X}:=\frac{1}{N}\sum_{i=1}^{N}Y_{i}Y_{i}^{T}$. The parameter $\beta$ is set to be $0.1$ for MULTADMM, $0.01$ for PSADMM and DSADMM. For PGADM, following , $\beta$ is set to be $0.1$, and the parameter $\tau$ in PGADM is set to be $0.6$.

Table 3 compares the algorithms with different choices of penalty parameters $\alpha_{1}$ and $\alpha_{2}$. We can find that the proposed methods PSADMM and DSADMM appear to perform better than the state-of-the-art method PGADM: similar objective function values are achieved using significantly less computational time and fewer iterations. Furthermore, MULTADMM is much faster than all the remaining three algorithms, although no theoretical results have been proved yet.

2.2. Stock Dataset

We then test the algorithms using real stock dataset, which includes the monthly stock return data of companies in the S$\&$P100 index from January 1990 to December 2012. We disregard $26$ companies that were listed after 1990 and only use the data of the remaining $p=74$ companies. The number of months $n$ is equal to 276. Each correlation coefficient between two stocks is computed based on the $n$ dimensional S$\&$P100 index vector of each stock. As a result the estimated covariance matrix is a $p\times p$ matrix with diagonal elements all being $1$. We choose the penalty parameters as $\alpha_{1}=0.005$ and $\alpha_{2}=0.01$; choose $\beta=10$ for all four methods; choose $\tau=0.6$ for PGADM. The performance of the algorithm is presented in Table 4. The identified relationship among the companies, characterized by the sparse part of the inverse of the estimated concentration matrix $\hat{\Sigma}_{X}$, is shown in Figure 2.

From Table 4, we see that MULTADMM still performs the best with substantially less computational time and fewer number of iterations, while the PSADMM and DSADMM slightly outperform PGADM. Figure 2 shows the identified relationship between some of those $74$ chosen companies in S$\&$P100. Recall that two companies $i,j$ are related if the $(i,j)$th entry of the computed sparse part $S$ is nonzero. All these methods are applied to the same data-set, as a result all the results are almost the same. So that we can choose any one of the obtained $S$ to plot Figure 2, and the MULTADMM result is chosen in this paper. In Figure 2, three groups of the companies are shown, and it is clear that the identified groups are intuitively meaningful: they represent Information Technology Companies, Retail Companies and Oil and Gas Companies respectively.

3. Robust Principal Component Analysis

For detailed discussion about this model, we refer the interested readers to .

Next we specialize Algorithm 1, Algorithm 2 and Algorithm 3 to solve problem (3.6).

For Algorithm 1, at the $k$th iteration, the variables are updated according to

where $N^{k}:=M+\frac{\lambda^{k}}{\beta}-L^{k+1}-S^{k+1}$.

For Algorithm 2, at $k$th iteration, the variables are updated according to

where $y_{i}$ is the slack variables defined Algorithm 2; $N^{k}:=\frac{M}{3}+y_{3}^{k+1}+\frac{\lambda_{3}^{k}}{\beta}$.

As for Algorithm 3, the variables are updated according to

where $t_{i}$ is the slack variable defined in Algorithm 3; $N^{k}:=\beta\lambda^{k+1}-t_{3}^{k}$.

The test dataset is a video taken at the hall of an aiport Available at http://perception.i2r.a-star.edu.sg/$\hbox{bk}\_\hbox{model}$/$\hbox{bk}\_\hbox{index}$.html. that consists of 200 grayscale frames, each of the size 144$\times$176. As a result the matrix $M$ is of dimension 25,344$\times$200. The index set $\Omega$ is determined randomly with a fixed sampling ratio $sr=80\%$, meaning that $20\%$ of the pixels are missing. The proposed three methods (MULTADMM, PSADMM and DSADMM) are compared with ADMMGBS in and VASALM in . The parameters in model (1.2) is set as $\tau=\frac{1}{\sqrt{l}}$ and $\delta=10^{-2}$. We choose $\beta=\frac{0.002}{\|M\|_{1}}$ for MULTADMM and VASALM, $\beta=\frac{0.004}{\|M\|_{1}}$ for PSADMM, and $\beta=2\cdot\|M\|_{1}$ for DSADMM. All the algorithms are initialized using zero matrices. The stopping criterions are set to be

where $\left\{f^{k}:=\|L^{k}\|_{*}+\tau\|S^{k}\|_{1}|\right\}$ is the sequence of objective value of (1). In Table 5, we have defined $\hbox{Obj}_{error}$ as the relative successive difference of the objective, i.e.,

Moreover, “Iter”, “Time”, “nnz(S)”, and “rank(L)” respectively denote the iteration number, the computation time, the number of nonzero elements in $S$, and the rank of $L$.

From the results in Table 5, once again we observe that the MULTADMM is the most efficient algorithm for solving RPCA, although the margin of advantage over the remaining algorithms is small.

Conclusion

In this paper, we systematically study several ADMM based algorithms for solving multiple block separable convex minimization problems. Two algorithms are proposed, which apply the classical two-block ADMM to either a primal or a dual reformulation of the original multi-block problem. Various theoretical properties of these two algorithms, such as their global convergence, iteration complexity and linear convergence, are shown by leveraging existing results. Through extensive numerical experiments, we observe that these two methods are competitive with state-of-the-art algorithms. Somewhat surprisingly, we show that in most cases, the classical multiple block ADMM method is computationally much more efficient than the proposed algorithms. This observation suggests that it is important to investigate the theoretical properties of the multiple block ADMM methods.

Acknowledgment: The authors are grateful to Professor Bingsheng He of Nanjing University and Dr. Tsung-Hui Chang of National Taiwan University of Science and Technology for valuable suggestions and discussions.

Appendix

Proof of Theorem 2.2. First recall that at ($k+1$)-th iteration, the optimality conditions of each subproblem about $x_{i}$ for both exact and inexact versions are given by

where $P_{i}=0$ for exact version and $P_{i}=\tau_{i}I-\beta A_{i}^{T}A_{i}$ for inexact version. Further substitute $\lambda_{i}^{k+1}=\lambda_{i}^{k}-\beta\left(A_{i}x_{i}^{k+1}-\frac{b}{m}-y_{i}^{k+1}\right)$ into the above inequalities,

then combine with the convexity of $f_{i}(\cdot)$ and let $x_{i}=x_{i}^{*}\in{\mathcal{X}}_{i}$

At $(k+1)$-th iteration, the optimality condition of the subproblem about $y$ is given by

after letting {\mbox{\boldmathy}}={\mbox{\boldmathy}}^{*}\in{\mathcal{Y}} in (5.2), we obtain

Summing (5.1) from $i=1$ to $i=m$, and using the above inequality, we can bound the difference of the current and optimal objective value by

while the last inequality is obtained from the multiplier update step

Proof of Theorem 2.3 for inexact version. First recall the optimality conditions of each subproblems about $x_{i}$ and $\lambda$ in the ($k+1$)-th iteration, which are expressed as $\forall\ x_{i}\in{\mathcal{X}}_{i}$

where $\lambda_{i}$ and $t_{i}$ also satisfy the following two equalities,

After substituting these two equalities into the above optimality conditions, we can obtain

Adding these inequalities together, we have

Set $x_{i}=x_{i}^{*}$ and $\lambda=\lambda^{*}$ where ($x_{1}^{*},\cdots,x_{m}^{*},\lambda^{*}$) is any saddle point of the Lagrangian function of (1), then clearly the left hand side of the above inequality becomes no greater than zero. Note that the following equality is true

where ${\mathcal{Q}}=\hbox{diag}({\mathcal{Q}}_{1},\cdots,{\mathcal{Q}}_{m})$. We can get the following result by summing from $1$ to $\infty$,

Using the equality relationship (5.3) we can also obtain,

Proof of Corollary 2.4. Let us rewrite the optimality conditions of each subproblems in Algorithm 3,

where $G_{i}=\frac{1}{\beta}A_{i}^{T}A_{i}$ for exact version and $G_{i}=\tau_{i}I$ for inexact version. Using the Lagrangian function of (1) and for all $(x_{1},\cdots,x_{m},\lambda)$, we obtain

Next we bound part of the right hand side in the above inequality. Notice that

where the last inequality is due to the definitions of $G_{i}$ and the assumptions that $\tau_{i}>\frac{\rho(A_{i}^{T}A_{i})}{\beta}$ for inexact version. As a result, we have

Further summing from $k=0$ to $k=K-1$, we obtain

References