Fast Alternating Linearization Methods for Minimizing the Sum of Two Convex Functions

Introduction

In this paper, we are interested in the following convex optimization problem:

In particular, we are specially interested in cases where solving (1.2) (or (1.3)) takes roughly the same effort as computing the gradient (or a subgradient) of $f(x)$ (or $g(x)$, respectively). Problems of this type arise in many applications of practical interest. The following are some interesting examples.

where $\rho>0$, which is of the form of (1.1) with $f(x)=\frac{1}{2}\|Ax-b\|_{2}^{2}$ and $g(x):=\rho\|x\|_{1}$. In this case, the two problems (1.2) and (1.3) are easy to solve. Specifically, (1.2) reduces to solving a linear system and (1.3) reduces to a vector shrinkage operation which requires $O(n)$ operations (see e.g., ). Depending on the size and structure of $A$ solving the system of linear equations required by (1.2) may be more expensive, less expensive or comparable to computing the gradient $A^{\top}(Ax-b)$ of $f(x)$. In the application we consider in Section 4 these computations are comparable due to the special structure of $A$.

Example 2. Nuclear norm minimization (NNM). The nuclear norm minimization problem, which seeks a low-rank solution of a linear system, can be cast as

Example 3. Robust principal component analysis (RPCA). The RPCA problem seeks to recover a low-rank matrix $X$ from a corrupted matrix $M$. This problem has many applications in computer vision, image processing and web data ranking (see e.g., ), and can be formulated as

which is of the form of (1.1). Moreover, the two problems (1.2) and (1.3) corresponding to (1.6) have closed-form solutions given respectively by a matrix shrinkage operation and a vector shrinkage operation. The matrix shrinkage operation requires a singular value decomposition (SVD) and is comparable in cost to computing a subgradient of $\|X\|_{*}$ or the gradient of the smoothed version of this function (see Section 5).

where $\rho>0$ and $\Sigma\in S_{+}^{n}$ (the set of symmetric positive semidefinite matrices) is the sample covariance matrix. Note that by defining $f(X):=-\log\det(X)+\langle\Sigma,X\rangle$ and $g(X):=\rho\|X\|_{1}$, (1.7) is of the form of (1.1). Moreover, it can be proved that the problem (1.2) has a closed-form solution, which is given by a spectral decomposition - a comparable effort to computing the gradient of $f(X)$, while the solution of problem (1.3) corresponds to a vector shrinkage operation.

Algorithms for solving problem (1.1) have been studied extensively in the literature. For large-scale problems, for which problems (1.2) and (1.3) are relatively easy to solve, the class of alternating direction methods that are based on variable splitting combined with the augmented Lagrangian method are particularly important. In these methods, one splits the variable $x$ into two variables, i.e., one introduces a new variable $y$ and rewrites Problem (1.1) as

Since Problem (1.8) is an equality constrained problem, the augmented Lagrangian method can be used to solve it. Given a penalty parameter $1/\mu$, at the $k$-th iteration, the augmented Lagrangian method minimizes the augmented Lagrangian function

with respect to $x$ and $y$, i.e., it solves the subproblem

and then updates the Lagrange multiplier $\lambda^{k}$ via:

Minimizing $\mathcal{L}_{\mu}(x,y;\lambda)$ with respect to $x$ and $y$ jointly is often not easy. In fact, it certainly is not any easier than solving the original problem (1.1). However, if one minimizes $\mathcal{L}_{\mu}(x,y;\lambda)$ with respect to $x$ and $y$ alternatingly, one needs to solve problems of the form (1.2) and (1.3), which as we have already discussed, is often easy to do. Such an alternating direction augmented Lagrangian method (ADAL) for solving (1.8) is given below as Algorithm 1.

Another important and related class of algorithms for solving (1.1) is based on operator-splitting. The aim of these algorithms is to find an $x$ such that

where $T_{1}$ and $T_{2}$ are maximal monotone operators. This is a more general problem than (1.1) and ADMs for it have been the focus of a substantial amount of research; e.g., see . Since the first-order optimality conditions for (1.1) are:

where $\partial f(x)$ denotes the subdifferential of $f(x)$ at the point $x$, a solution to Problem (1.1) can be obtained by solving Problem (1.12). For example, see and references therein for more information on this class of algorithms.

While global convergence results for various splitting and alternating direction algorithms have been established under appropriate conditions, our interest here is on iteration complexity bounds for such algorithms. By an iteration complexity bound we mean a bound on the number of iterations needed to obtain an $\epsilon$-optimal solution which is defined as follows.

Complexity bounds for first-order methods for solving convex optimization problems have been given by Nesterov and many others. In , Nesterov gave first-order algorithms for solving smooth unconstrained convex minimization problems with an iteration complexity of $O(\sqrt{L/\epsilon})$, where $L$ is the Lipschitz constant of the gradient of the objective function, and showed that this is the best complexity that is obtainable when only first-order information is used. These methods can be viewed as accelerated gradient methods where a combination of past iterates are used to compute the next iterate. Similar techniques were then applied to nonsmooth problems and corresponding optimal complexity results were obtained. The ISTA (Iterative Shrinkage/Thresholding Algorithm) and FISTA (Fast Iterative Shrinkage/Thresholding Algorithm) algorithms proposed by Beck and Teboulle in are designed for solving (1.1) when one of the functions (say $f(x)$) is smooth and the other is not. It is proved in that the number of iterations required by ISTA and FISTA to get an $\epsilon$-optimal solution to problem (1.1) are respectively $O(L(f)/\epsilon)$ and $O(\sqrt{L(f)/\epsilon})$, under the assumption that $\nabla f(x)$ is Lipschitz continuous with Lipschitz constant $L(f)$, i.e.,

ISTA computes a sequence $\{x^{k}\}$ via the iteration

while FISTA computes $\{x^{k}\}$ via the iteration

Note that ISTA and FISTA treat the functions $f(x)$ and $g(x)$ very differently. At each iteration they both linearize the function $f(x)$ but never directly minimize it, while they do minimize the function $g(x)$ in conjunction with the linearization of $f(x)$ and a proximal (penalty) term. These two methods have proved to be efficient for solving the CS problem (1.4) (see e.g., ) and the NNM problem (1.5) (see e.g., ). ISTA and FISTA work well in these areas because $f(x)$ is quadratic and is well approximated by linearization. However, for the RPCA problem (1.6) where two complicated functions are involved, ISTA and FISTA do not work well. As we shall show in Sections 2 and 3, our ADMs are very effective in solving RPCA problems. For the SICS problem (1.7), intermediate iterates $X^{k}$ may not be positive definite, and hence the gradient of $f(X)=-\log\det(X)+\langle\Sigma,X\rangle$ may not be well defined at $X^{k}$. Therefore, ISTA and FISTA cannot be used to solve the SICS problem (1.7). In , it is shown that SICS problems can be very efficiently solved by our ADM approach.

Our contribution. In this paper, we propose both basic and accelerated (i.e., fast) versions of first-order alternating linearization methods (ALMs) based on an alternating direction augmented Lagrangian approach for solving (1.1) and analyze their iteration complexities. Our basic methods require at most $O(L/\epsilon)$ iterations to obtain an $\epsilon$-optimal solution, while our fast methods require at most $O(\sqrt{L/\epsilon})$ iterations with only a very small increase in the computational effort required at each iteration. Thus, our fast methods are optimal first-order methods in terms of iteration complexity. For both types of methods, we present an algorithm that requires both functions to be continuously differentiable with Lipschitz constants for the gradients denoted by $L(f)$ and $L(g)$. In this case $L=\max\{L(f),L(g)\}$. We also present for each type of method, an algorithm that only needs one of the functions, say $f(x)$, to be smooth, in which case $L=L(f)$. These algorithms are related to the multiple splitting algorithms in a recent paper by Goldfarb and Ma . The algorithms in are Jacobi type methods since they do not use information from the current iteration to solve succeeding subproblems in that iteration, while the algorithms proposed in this paper are Gauss-Seidel type methods since information from the current iteration is used later in the same iteration. These algorithms can also be viewed as extensions of the ISTA and FISTA algorithms in . The complexity bounds we obtain for our algorithms are similar to (and as much as a factor of two better that) those in .

At each iteration, our algorithms alternatively minimize two different approximations to the original objective function, obtained by keeping one function unchanged and linearizing the other one. Our basic algorithm is similar in many ways to the alternating linearization method proposed by Kiwiel et al.. In particular, the approximate functions minimized at each step of Algorithm 3.1 in have the same form as those minimized in our algorithm. However, our basic algorithm differs from the one in in the way that the proximal terms are chosen, and our accelerated algorithms are very different. Moreover, no complexity bounds have been given for the algorithm in . To the best of our knowledge, the complexity results in this paper are the first ones that have been given for a Gauss-Seidel type alternating direction method After completion of an earlier version of the present paper, which is available on http://arxiv.org/abs/0912.4571, Monteiro and Svaiter gave an iteration complexity bound to achieve a desired closeness of the current iterate to the solution for ADMs for solving the more general problem (1.12).. Complexity results for related Jacobi type alternating direction methods are given in .

Organization. The rest of this paper is organized as follows. In Sections 2 and 3 we propose our alternating linearization methods based on alternating direction augmented Lagrangian methods and give convergence/complexity bounds for them. We compare the performance of our ALMs to other competing first-order algorithms using an image deblurring problem in Section 4. In Section 5, we apply our ALMs to solve very large RPCA problems arising from background extraction in surveillance video and matrix completion and report the numerical results. Finally, we make some conclusion in Section 6.

Alternating Linearization Methods

In iteration of the ADAL method, Algorithm 1, the Lagrange multiplier $\lambda$ is updated just once, immediately after the augmented Lagrangian is minimized with respect to $y$. Since the alternating direction approach is meant to be symmetric with respect to $x$ and $y$, it is natural to also update $\lambda$ after solving the subproblem with respect to $x$. By doing this, we get a symmetric version of the ADAL method. This algorithm is given below as Algorithm 2.

This ADAL variant is described and analyzed in . Moreover, it is shown in that Algorithms 1 and 2 are equivalent to the Douglas-Rachford and Peaceman-Rachford methods, respectively applied to the optimality condition (1.13) for problem (1.1). If we assume that both $f(x)$ and $g(x)$ are differentiable, it follows from the first order optimality conditions for the two subproblems in lines 3 and 5 of Algorithm 2 that

Substituting (2.1) into Algorithm 2, we get the following alternating linearization method (ALM) which is equivalent to the SADAL method (Algorithm 2) when both $f$ and $g$ are differentiable.

In Algorithm 3, $Q_{f}(u,v)$ is defined by (1.15) and

In Algorithm 3, we alternatively replace the functions $g$ and $f$ by their linearizations plus a proximal regularization term to get an approximation to the original function $F$. Thus, our ALM algorithm can also be viewed as a proximal point algorithm.

A drawback of Algorithm 3 is that it requires both $f$ and $g$ to be continuously differentiable. In many applications, however, one of these functions is nonsmooth, as in the examples given in Section 1. Although Algorithm 2 can be applied when $f(x)$ and $g(x)$ are nonsmooth, we are unable to provide a comparable complexity bound in this case. However, when only one of the functions of $f$ and $g$ is nonsmooth (say $g$ is nonsmooth), the following variant of Algorithm 3 applies, and for this algorithm, we have a complexity result.

We call Algorithm 4, ALM with skipping steps (ALM-S) because in line 4 of Algorithm 4, if

holds, we let $x^{k+1}:=y^{k}$, i.e., we skip the computation of $x^{k+1}$ in line 3. An alternative version of Algorithm 4 that has smaller average work per iteration is the following Algorithm 5.

Note that in Algorithm 5, when (2.3) does not hold, we switch to the SADAL algorithm, which updates $\lambda$ instead of computing $\nabla f(x^{k+1})$. Algorithm 5 is usually faster than Algorithm 4 when $\nabla f(x^{k+1})$ is costly to compute in addition to performing Step 3. Note also that when (2.3) holds, then the steps of the algorithm reduce to those of ISTA.

The following theorem gives conditions under which Algorithms 2, 3, 4 and 5 are equivalent.

(i) If both $f$ and $g$ are differentiable, and $\lambda_{0}$ is set to $-\nabla g(y^{0})$, then Algorithms 2 and 3 are equivalent. (ii) If in addition $g$ is Lipschitz continuous with Lipschitz constant $L(g)$, and $\mu\leq 1/L(g)$, then Algorithms 3 and 4 are equivalent. (iii) If $f$ is differentiable, then Algorithms 4 and 5 are equivalent.

When both $f$ and $g$ are differentiable and $\lambda_{0}=-\nabla g(y^{0})$, (2.1) holds for all $k\geq 0$, and it follows that $\mathcal{L}_{\mu}(x,y^{k};\lambda^{k})\equiv Q_{g}(x,y^{k})$ and $\mathcal{L}_{\mu}(x^{k+1},y;\lambda^{k+\frac{1}{2}})\equiv Q_{f}(y,x^{k+1}).$ This proves part (i). If $\nabla g(x)$ is Lipschitz continuous and $\mu\leq 1/L(g)$,

holds (see e.g., ). This implies that (2.3) does not hold and hence, $x^{k+1}:=\arg\min_{x}\mathcal{L}_{\mu}(x,y^{k};\lambda^{k})$, and the equivalence of Algorithms 3 and 4 follows. This proves part (ii). The optimality of $x^{k+1}$ in line 3 of Algorithm 5 implies that $\lambda^{k+\frac{1}{2}}=\nabla f(x^{k+1})$ when (2.3) does not hold and hence that $\mathcal{L}_{\mu}(x^{k+1},y;\lambda^{k+\frac{1}{2}})\equiv Q_{f}(y,x^{k+1})$. This proves part (iii). ∎

We show in the following that the iteration complexity of Algorithm 4 is $O(1/\epsilon)$ for obtaining an $\epsilon$-optimal solution for (1.1). First, we need the following generalization of Lemma 2.3 in .

where $\gamma_{\psi}(v)$ is any subgradient in the subdifferential $\partial\psi(v)$ of $\psi(v)$ at the point $v$. Let $\Phi(\cdot)=\phi(\cdot)+\psi(\cdot)$. For any $v$, if

Since $\phi$ and $\psi$ are convex we have

where $\gamma_{\phi}(\cdot)$ is a subgradient of $\phi(\cdot)$ and $\gamma_{\phi}(p_{\psi}(v))$ satisfies the first-order optimality conditions for (2.4), i.e.,

Therefore, from (2.9), (2.12) and (2.13) it follows that

Assume $\nabla f(\cdot)$ is Lipschitz continuous with Lipschitz constant $L(f)$. For $\mu\leq 1/L(f)$, the iterates $y^{k}$ in Algorithm 4 satisfy

where $x^{*}$ is an optimal solution of (1.1) and $k_{n}$ is the number of iterations until the $k$-th for which $F(x^{k+1})\leq\mathcal{L}_{\mu}(x^{k+1},y^{k};\lambda^{k})$, i.e., the number of iterations when no skipping step occurs. Thus, the sequence $\{F(y^{k})\}$ produced by Algorithm 4 converges to $F(x^{*})$. Moreover, if $1/(\beta L(f))\leq\mu\leq 1/L(f)$ where $\beta\geq 1$, the number of iterations needed to obtain an $\epsilon$-optimal solution is at most $\lceil C/\epsilon\rceil$, where $C=\beta L(f)\|x^{0}-x^{*}\|^{2}/2$.

Let $I$ be the set of all iteration indices until $k-1$-st for which no skipping occurs and let $I_{c}$ be its complement. Let $I=\{n_{i}\},\ i=0,\ldots,k_{n}-1$. It follows that for all $n\in I_{c}$, $x^{n+1}=y^{n}$.

For $n\in I$ we can apply Lemma 3 to obtain the following inequalities. In (2.6), by letting $\psi=f$, $\phi=g$, $u=x^{*}$ and $v=x^{n+1}$, we get $p_{\psi}(v)=y^{n+1}$, $\Phi=F$ and

Similarly, by letting $\psi=g$, $\phi=f$, $u=x^{*}$ and $v=y^{n}$ in (2.6) we get $p_{g}(v)=x^{n+1}$, $\Phi=F$ and

Taking the summation of (2.16) and (2.17) we get

For $n\in I_{c}$, (2.16) holds. Then since $x^{n+1}=y^{n}$ we get

Summing (2.18) and (2.19) over $n=0,1,\ldots,k-1$ we get

For any $n$, since Lemma 3 holds for any $u$, letting $u=x^{n+1}$ instead of $x^{*}$ we get from (2.16) that

Thus we get $F(y^{n})\leq F(x^{n}),\forall n.$

Similarly, for $n\in I$ by letting $u=y^{n}$ instead of $x^{*}$ we get from (2.17) that

On the other hand, for $n\in I_{c}$, (2.23) holds trivially because $x^{n+1}=y^{n}$; thus (2.23) holds for all $n$.

Adding (2.21) and (2.23) and adding (2.22) and (2.23), respectively, yield

The inequalities (2.24) show that the sequences of function values $F(y^{n})$ and $F(x^{n})$ are non-increasing. Thus we have,

which gives us the desired result (2.15). ∎

Assume $\nabla f$ and $\nabla g$ are both Lipschitz continuous with Lipschitz constants $L(f)$ and $L(g)$, respectively. For $\mu\leq\min\{1/L(f),1/L(g)\}$, Algorithm 3 satisfies

where $x^{*}$ is an optimal solution of (1.1). Thus sequence $\{F(y^{k})\}$ produced by Algorithm 3 converges to $F(x^{*})$. Moreover, if $1/(\beta\max\{L(f),L(g)\})\leq\mu\leq 1/\max\{L(f),L(g)\}$ where $\beta\geq 1$, the number of iterations needed to get an $\epsilon$-optimal solution is at most $\lceil C/\epsilon\rceil$, where $C=\beta\max\{L(f),L(g)\}\|x^{0}-x^{*}\|^{2}/4$.

The conclusion follows from Theorems 2 and 4 and $k_{n}=k$. ∎

The complexity bound in Corollary 5 is smaller than the analogous bound for ISTA in by a factor of two. It is easy to see that the bound in Theorem 4 is also an improvement over the bound in as long as $F(x^{k+1})\leq\mathcal{L}_{\mu}(x^{k+1},y^{k};\lambda^{k})$ holds for at least one value of $k$. It is reasonable then to ask if the per-iteration cost of Algorithms 3 and 4 are comparable to that of ISTA. It is indeed the case when the assumption holds that minimizing $\mathcal{L}_{\mu}(x,y^{k};\lambda^{k})$ has comparable cost (and often involves the same computations) as computing the gradient $\nabla f(y^{k})$.

are easily handled by our approach, since one can express (2.30) as

If a convex constraint $x\in\mathcal{C}$, where $\mathcal{C}$ is a convex set is added to problem (1.1), and we impose this constraint in the two subproblems in Algorithms 3 and 4, i.e., we impose $x\in\mathcal{C}$ in the subproblems with respect to $x$ and $y\in\mathcal{C}$ in the subproblems with respect to $y$, the complexity results in Theorem 4 and Corollary 5 continue to hold. The only changes in the proof are in Lemma 3. If there is a constraint $x\in\mathcal{C}$, then (2.6) holds for any $u\in\mathcal{C}$ and $v\in\mathcal{C}$. Also in the proof of Lemma 3, the first equality in (2.14) becomes a “$\geq$” inequality due to the fact that the optimality conditions (2.12) become

Although Algorithms 3 and 4 assume that the Lipschitz constants are known, and hence that an upper bound for $\mu$ is known, this can be relaxed by using the backtracking technique in to estimate $\mu$ at each iteration.

Fast Alternating Linearization Methods

In this section, we propose a fast alternating linearization method (FALM) which computes an $\epsilon$-optimal solution to problem (1.1) in $O(\sqrt{L/\epsilon})$ iterations, while keeping the work at each iteration almost the same as that required by ALM.

FALM is an accelerated version of ALM for solving (1.1), or equivalently (1.8), when $f(x)$ and $g(x)$ are both differentiable, and is given below as Algorithm 6. Clearly, FALM is also a Gauss-Seidel type algorithm. In fact, it is a successive over-relaxation type algorithm since $(t_{k}-1)/t_{k+1}>0,\forall k\geq 2.$

Algorithm 6 requires both $f$ and $g$ to be continuously differentiable. To develop an algorithm that can be applied to problems where one of the functions is non-differentiable, we use a skipping technique as in Algorithm 4. FALM with skipping steps (FALM-S), which does not require $g(x)$ to be smooth, is given below as Algorithm 7.

The following theorem gives conditions under which Algorithms 6 and 7 are equivalent.

If both $f(x)$ and $g(x)$ are differentiable and $\nabla g(x)$ is Lipschitz continuous with Lipschitz constant $L(g)$, and $\mu\leq 1/L(g)$, then Algorithms 6 and 7 are equivalent.

As in Theorem 2, if $f$ and $g$ are differentiable, $\mathcal{L}_{\mu}(x,z^{k};\lambda^{k})\equiv Q_{g}(x,z^{k})$. The conclusion then follows from the fact that $F(x^{k})\leq\mathcal{L}_{\mu}(x^{k},z^{k};\lambda^{k})$ always holds when $\mu\leq 1/L(g)$ and thus there are no skipping steps. ∎

To prove that Algorithm 7 requires $O(\sqrt{L(f)/\epsilon})$ iterations to obtain an $\epsilon$-optimal solution, we need the following lemmas. We call $k$-th iteration a skipping step if $x^{k}=z^{k}$, and a regular step if $x^{k}\neq z^{k}$.

The sequence $\{x^{k},y^{k}\}$ generated by Algorithm 7 satisfies

where $u^{k}:=t_{k}y^{k}-(t_{k}-1)y^{k-1}-x^{*}$ and $v_{k}:=2F(y^{k})-2F(x^{*})$ if iteration $k$ is a regular step and $v_{k}:=F(y^{k})-F(x^{*})$ if iteration $k$ is a skipping step.

There are four cases to consider: (i) both the $k$-th and the $(k+1)$-st iterations are regular steps; (ii) the $k$-th iteration is a regular step and the $(k+1)$-st iteration is a skipping step; (iii) both the $k$-th and the $(k+1)$-st iterations are skipping steps; (iv) the $k$-th iteration is a skipping step and the $(k+1)$-st iteration is a regular step. We will prove that the following inequality holds for all the four cases:

The proof of (3.1) and hence, the lemma, then follows from the fact that the right hand side of inequality (3.2) equals

where we have used the fact that $t_{k+1}z^{k+1}:=t_{k+1}y^{k}+t_{k}(y^{k}-y^{k-1})-(y^{k}-y^{k-1}).$

Case (i): Let us first consider the case when both the $k$-th and $(k+1)$-st iterations are regular steps. In (2.6), by letting $\psi=f$, $\phi=g$, $u=y^{k}$ and $v=x^{k+1}$, we get $p_{\psi}(v)=y^{k+1}$, $\Phi=F$ and

In (2.6), by letting $\psi=g$, $\phi=f$, $u=y^{k}$, $v=z^{k+1}$, we get $p_{\psi}(v)=x^{k+1}$, $\Phi=F$ and

Summing (3.3) and (3.4), and using the fact that $F(y^{k+1})\leq F(x^{k+1})$, we obtain,

Again, in (2.6), by letting $\psi=g$, $\phi=f$, $u=x^{*}$, $v=z^{k+1}$, we get $p_{\psi}(v)=x^{k+1}$, $\Phi=F$ and

In (2.6), by letting $\psi=f$, $\phi=g$, $u=x^{*}$, $v=x^{k+1}$, we get $p_{\psi}(v)=y^{k+1}$, $\Phi=F$ and

Summing (3.6) and (3.7), and again using the fact that $F(y^{k+1})\leq F(x^{k+1})$, we obtain,

If we multiply (3.5) by $t_{k}^{2}$, and (3.8) by $t_{k+1}$, and take the sum of the resulting two inequalities, we get (3.2) by using the fact that $t_{k}^{2}=t_{k+1}(t_{k+1}-1)$.

Case (ii): By letting $\psi=f$, $\phi=g$, $u=y^{k}$ and $v=z^{k+1}$ in (2.6), we get $p_{\psi}(v)=y^{k+1}$, $\Phi=F$ and

Since the steps taken in the $k$-th and $(k+1)$-st iterations are regular and skipping, respectively, we have

Also by letting $\psi=f$, $\phi=g$, $u=x^{*}$ and $v=z^{k+1}$ in (2.6), we get $p_{\psi}(v)=y^{k+1}$, $\Phi=F$ and

Then multiplying (3.10) by $2t_{k}^{2}$, (3.11) by $t_{k+1}$, summing the resulting two inequalities and using the fact that in this case $2t_{k}^{2}=t_{k+1}(t_{k+1}-1)$, we obtain (3.2).

Case (iii): This case reduces to two consecutive FISTA steps and the proof above applies with $t_{k}^{2}=t_{k+1}(t_{k+1}-1)$ and inequality (3.10) replaced by

Case (iv): In this case, (3.5) in the proof of case (i) is replaced by

which when multiplied by $t_{k}^{2}/2$ and combined with (3.8) multiplied by $t_{k+1}$, and the fact that in this case $t_{k}^{2}/2=t_{k+1}(t_{k+1}-1)$, yields (3.2). ∎

The following lemma gives lower bounds for the sequence of scalars $\{t_{k}\}$ generated by Algorithm 7.

For all $k\geq 1$ the sequence $\{t_{k}\}$ generated by Algorithm 7 satisfies:

where $r(k)$ and $s(k)$ are the number of steps among the first $k$ steps that are regular and skipping steps, respectively, and $\alpha\equiv\sqrt{2}-1$ and $\hat{\alpha}\equiv\frac{1}{\sqrt{2}}-1$.

Consider the case where the first iteration of Algorithm 7 is a skipping step. Clearly, the sequence of iterations follows a pattern of alternating blocks of one or more skipping steps and one or more regular steps. Let the index of the first iteration in the $i$-th block be denoted by $n_{i}$. Since it is assumed that the first iteration is a skipping step, iterations $n_{1},n_{3},n_{5}\ldots$ are skipping steps ($n_{1}=1$) and $n_{2},n_{4},n_{6}\ldots$ are regular steps. Note that the statement of the lemma in this case corresponds to

We first note that it follows from the updating rules and formulas for $t_{k}$ that

Consider $j=1$. Clearly (3.18) holds for all iterations $n_{1}=1\leq k\leq n_{2}-1$, since $t_{k}\geq\frac{k+1}{2}$ holds trivially for $t_{1}=1$, and for $1<k\leq n_{2}-1$, $t_{k}\geq\frac{k-1}{2}+t_{1}=\frac{k+1}{2}$.

Now assume that (3.18) holds for all $j<\bar{j}$. If $\bar{j}$ is even, iterations $k=n_{\bar{j}}$ and $k-1$ are, respectively, regular and skipping iterations. Hence, from (3.22), we have that

Since the remaining $p\equiv n_{{\bar{j}}+1}-1-n_{\bar{j}}$ iterations before iteration $n_{{\bar{j}}+1}$ are all regular iterations ($p$ may be zero), we have from (3.22) that, for $n_{\bar{j}}<k\leq n_{\bar{j}+1}-1$,

If $\bar{j}$ is odd, iteration $k=n_{\bar{j}}$ is a skipping iteration. Hence, from (3.22), we have that $t_{k}\geq\frac{1}{2}+\sqrt{2}t_{k-1}\geq\frac{1}{2}+\sqrt{2}\frac{1}{2\sqrt{2}}(k+\alpha r(k-1))=\frac{1}{2}(k+1+\alpha r(k))$. Since the remaining $p\equiv n_{{\bar{j}}+1}-1-n_{\bar{j}}$ iterations before iteration $n_{{\bar{j}}+1}$ are all skipping iterations (again $p$ may be zero), we have from (3.22) that, for $n_{\bar{j}}<k\leq n_{\bar{j}+1}-1$, $t_{k}\geq\frac{k-n_{\bar{j}}}{2}+t_{n_{\bar{j}}}\geq\frac{k-n_{\bar{j}}}{2}+\frac{1}{2}(n_{\bar{j}}+1+\alpha r(n_{\bar{j}}))=\frac{1}{2}(k+1+\alpha r(k))$. This concludes the induction.

Since the proof for the case that the first step is a regular step is totally analogous to the above proof, we leave this to the reader. ∎

Now we are ready to give the complexity of Algorithm 7.

Let $\alpha=\sqrt{2}-1$ and $r(k)$ be the number of steps among the first $k$ steps that are regular steps. Assuming $\nabla f(\cdot)$ is Lipschitz continuous with Lipschitz constant $L(f)$, if $\mu\leq 1/L(f)$, the sequence $\{y^{k}\}$ generated by Algorithm 7 satisfies:

where $\hat{r}(k)=r(k)$ if the first step is a skipping step, and $\hat{r}(k)=r(k)+1$ if the first step is a regular step.

Hence, the sequence $\{F(y^{k})\}$ produced by Algorithm 4 converges to $F(x^{*})$. Moreover, if $1/(\beta L(f))\leq\mu\leq 1/L(f)$ where $\beta\geq 1$, the number of iterations required by Algorithm 7 to get an $\epsilon$-optimal solution to (1.1) is at most $\lfloor\sqrt{C/\epsilon}\rfloor$, where $C=2\beta L(f)\|x^{0}-x^{*}\|^{2}$.

Using the same notation as in Lemmas 11 and 12, (3.8) and (3.11) imply that

holds whether the first iteration is a skipping step or not. Thus we have

From Lemma 11 we know that the sequence $\{2\mu t_{k}^{2}v_{k}+\|u^{k}\|^{2}\}$ is non-increasing. Therefore, we have

where the equality follows from the facts that $t_{1}=1$ and $u^{1}=y^{1}-x^{*}$, and the last inequality is from (3.25).

Recall the definition of $v_{k}$ in Lemma 11 and the bounds of $t_{k}$ in Lemma 12. We get (i) and (ii) from (3.26). Keeping in mind that $v_{k}$ has a different expression depending on whether the $k$-th step is a skipping or a regular step, it follows that the sequence $\{y^{k}\}$ generated by Algorithm 7 satisfies:

(i) if the first step is a skipping step, then

(ii) if the first step is a regular step, then

It is easy to check that these bounds are equivalent to (3.24), and that the worst case bound on the number of iterations follows from (3.24). ∎

Assume $\nabla f$ and $\nabla g$ are both Lipschitz continuous with Lipschitz constants $L(f)$ and $L(g)$, respectively. For $\mu\leq\min\{1/L(f),1/L(g)\}$, Algorithm 6 satisfies

where $x^{*}$ is an optimal solution of (1.1). Hence, the sequence $\{F(y^{k})\}$ produced by Algorithm 6 converges to $F(x^{*})$, and if $1/(\beta\max\{L(f),L(g)\})\leq\mu\leq 1/\max\{L(f),L(g)\}$, where $\beta\geq 1$, the number of iterations needed to get an $\epsilon$-optimal solution is at most $\lceil\sqrt{C/\epsilon}-1\rceil$, where $C=\beta\max\{L(f),L(g)\}\|x^{0}-x^{*}\|^{2}$.

Note that since $\mu\leq\min\{1/L(f),1/L(g)\}$, from Theorem 10 we know that Algorithms 6 and 7 are equivalent. That is, every step in Algorithm 7 is a regular step. Therefore, case (ii) in Theorem 13 holds and $s(k)=0$, which leads to (3.27). ∎

The complexity bound in Corollary 14 is smaller than the analogous bound for FISTA in by a factor of $\sqrt{2}$. It is easy to see that the bound in Theorem 13 is also an improvement over the bound in as long as $F(x^{k+1})\leq\mathcal{L}_{\mu}(x^{k+1},y^{k};\lambda^{k})$ holds for at least one value of $k$. As in the case of Algorithms 3 and 4, the per-iteration cost of Algorithms 6 and 7 are comparable to that of FISTA.

Line 6 in Algorithm 6 and Line 15 in Algorithm 7 can be changed to:

where $w^{k}:=\alpha x^{k}+(1-\alpha)y^{k},\alpha\in(0,1)$, and Theorem 13 and Corollary 14 still hold.

Although Algorithms 6 and 7, assume that the Lipschitz constants are known, and hence that an upper bound for $\mu$ is known, this can be relaxed by using the backtracking technique in to estimate $\mu$ at each iteration.

Comparison of ALM, FALM, ISTA, FISTA, SADAL and SALSA

In this section we compare the performance of our basic and fast ALMs, with and without skipping steps, against ISTA, FISTA, SADAL (Algorithm 2) and an alternating direction augmented Lagrangian method SALSA described in on a benchmark wavelet-based image deblurring problem from . In this problem, the original image is the well-known Cameraman image of size $256\times 256$ and the observed image is obtained after imposing a uniform blur of size $9\times 9$ (denoted by the operator $R$) and Gaussian noise (generated by the function randn in MATLAB with a seed of 0 and a standard deviation of $0.56$). Since the coefficient of the wavelet transform of the image is sparse in this problem, one can try to reconstruct the image $u$ from the observed image $b$ by solving the problem:

It is easy to show that the optimal solution $z_{\sigma}(x)$ of (4.2) is

According to Theorem 1 in , the gradient of $g_{\sigma}$ is given by $\nabla g_{\sigma}(x)=z_{\sigma}(x)$ and is Lipschitz continuous with Lipschitz constant $L(g_{\sigma})=1/\sigma$. After smoothing $g$, we can apply Algorithms 3 and 6 to solve the smoothed problem:

We have the following theorem about the $\epsilon$-optimal solutions of problems (4.1) and (4.4).

Let $\sigma=\frac{\epsilon}{n\rho^{2}}$ and $\epsilon>0$. If $x(\sigma)$ is an $\epsilon/2$-optimal solution to (4.4), then $x(\sigma)$ is an $\epsilon$-optimal solution to (4.1).

Proof. Let $D_{g}:=\max\{\frac{1}{2}\|z\|_{2}^{2}:\|z\|_{\infty}\leq\rho\}=\frac{1}{2}n\rho^{2}$ and $x^{*}$ and $x^{*}(\sigma)$ be optimal solution to problems (4.1) and (4.4), respectively. Note that

Using the inequalities in (4.5) and the facts that $x(\sigma)$ is an $\epsilon/2$-optimal solution to (4.4) and $\sigma D_{g}=\frac{\epsilon}{2}$, we have

Thus, to find an $\epsilon$-optimal solution to (4.1), we can apply Algorithms 3 and 6 to find an $\epsilon/2$-optimal solution to (4.4) with $\sigma=\frac{\epsilon}{n\rho^{2}}$. The iteration complexity results in Corollaries 5 and 14 hold since the gradient of $g_{\sigma}$ is Lipschitz continuous. However, the numbers of iterations needed by ALM and FALM to obtain an $\epsilon$-optimal solution to (4.1) become $O(1/\epsilon^{2})$ and $O(1/\epsilon)$, respectively, due to the fact that the Lipschitz constant $L(g_{\sigma})=1/\sigma=\frac{n\rho^{2}}{\epsilon}=O(1/\epsilon)$.

When Algorithms 4 and 7 are applied to solve (4.1), the subproblems (1.2) and (1.3) are easy to solve. Specifically, (1.2) corresponds to solving a linear system which is particularly easy to do because of the special structures of $R$ and $W$ (see ); (1.3) corresponds to a vector shrinkage operation. When Algorithms 3 and 6 are applied to solve (4.4), (1.3) with $g$ replaced by $g_{\sigma}$ is also easy to solve; its optimal solution is

Since ALM is equivalent to SADAL when both functions are smooth, we implemented ALM as SADAL when we solved (4.4). We also applied SADAL to the nonsmooth problem (4.1). We also implemented ALM-S as Algorithm 5 since the latter was usually faster. In all algorithms, we set the initial points $x^{0}=y^{0}=\mathbf{0}$, and in FALM and FALM-S we set $z^{1}=\mathbf{0}$. MATLAB codes for SALSA, FISTA and ISTA (modified from FISTA) were downloaded from http://cascais.lx.it.pt/$\sim$mafonso/salsa.html and their default inputs were used. Moreover, $\lambda^{0}$ was set to $\mathbf{0}$ in algorithms 2, 4, 5 and 7 since $\nabla g_{\sigma}(x^{0})=\mathbf{0}$ and $\mathbf{0}\in-\partial g(x^{0})$ when $x^{0}=\mathbf{0}$. Also, whenever $g(x)$ was smoothed, we set $\sigma=10^{-6}$. $\mu$ was set to $1$ in all the algorithms since the Lipschitz constant of the gradient of function $\frac{1}{2}\|RW(\cdot)-b\|_{2}^{2}$ was known to be $1$. We set $\mu$ to 1 even for the smoothed problems. Although this violates the requirement $\mu\leq\frac{1}{L(g_{\sigma})}$ in Corollaries 5 and 14, we see from our numerical results reported below that ALM and FALM still work very well. All of the algorithms tested were terminated after 1000 iterations. The (nonsmoothed) objective function values in (4.1) produced by these algorithms at iterations: 10, 50, 100, 200, 500, 800 and 1000 for different choices of $\rho$ are presented in Tables 1 and 2. The CPU times (in seconds) and the number of iterations required to reduce the objective function value to below $1.04e+5$ and $8.60e+5$ are reported respectively in the last columns of Tables 1 and 2.

All of our codes were written in MATLAB and run in MATLAB 7.3.0 on a Dell Precision 670 workstation with an Intel Xeon(TM) 3.4GHZ CPU and 6GB of RAM.

From Tables 1 and 2 we see that in terms of the value of the objective function achieved after a specified number of iterations, the performance of FALM-S and FALM is always slightly better than that of FISTA and is much better than the performance of the other algorithms. On the two test problems, since FALM-S and FALM are always better than ALM-S and ALM, and FISTA is always better than ISTA, we can conclude that the Nesterov-type acceleration technique greatly speeds up the basic algorithms on these problems. Moreover, although sometimes in the early iterations FISTA (ISTA) is better than FALM-S and FALM (ALM-S and ALM), it is always worse than the latter two algorithms when the iteration number is large. We also illustrate our comparisons graphically by plotting in Figure 1 the objective function value versus the number of iterations taken by these algorithms for solving (4.1) with $\rho=0.1$. From Figure 1 we see clearly that for this problem, ALM outperforms ISTA, FALM outperforms FISTA, FALM outperforms ALM and FALM-S outperforms ALM-S.

From the CPU times and the iteration numbers in the last columns of Tables 1 and 2 we see that, the fast versions are always much better than the basic versions of the algorithms. Since iterations of FISTA cost less than those of FALM-S (and FALM as well), we see that although FISTA takes 35% (4%) more iterations than FALM-S in the last column of Table 1 (2) it takes only 7% more time (20% less time).

We note that for the problems with $\rho=0.01$ and $\rho=0.1$, (i.e., for the results given in Tables 1 and 2), 891 and 981, respectively, of the first 1000 iterations performed by FALM-S were skipping steps. In contrast, none of the steps performed by ALM-S on either of these problems were skipping steps. While the latter result is somewhat surprising, the fact that FALM-S performs many skipping steps is not, since the Nesterov-like acceleration approach is an over-relaxation approach that generates points that extrapolate beyond the previous point and the one produced by the ALM algorithm.

Applications

In this section, we describe how ALM and FALM can be applied to problems that can be formulated as RPCA problems to illustrate the use of Nesterov-type smoothing when the functions $f$ and $g$ do not satisfy the smoothness conditions required by the theorems in Sections 2 and 3. Our numerical results show that our methods are able to solve huge problems that arise in practice; e.g., one problem involving roughly 40 million variables and 20 million linear constraints is solved in about three-quarters of an hour. We alse describe application of our methods to the SICS problem.

It is easy to show that the optimal solution $W_{\sigma}(X)$ of (5.1) is

where $U\mbox{Diag}(\gamma)V^{\top}$ is the singular value decomposition (SVD) of $X/\sigma.$ According to Theorem 1 in , the gradient of $f_{\sigma}$ is given by $\nabla f_{\sigma}(X)=W_{\sigma}(X)$ and is Lipschitz continuous with Lipschitz constant $L(f_{\sigma})=1/\sigma$. After smoothing $f$ and $g$, we can apply Algorithms 3 and 6 to solve the following smoothed problem:

We have the following theorem about $\epsilon$-optimal solutions of problems (1.6) and (5.3).

Let $\sigma=\frac{\epsilon}{2\max\{\min\{m,n\},mn\rho^{2}\}}$ and $\epsilon>0$. If $(X(\sigma),Y(\sigma))$ is an $\epsilon/2$-optimal solution to (5.3), then $(X(\sigma),Y(\sigma))$ is an $\epsilon$-optimal solution to (1.6).

Proof. Let $D_{f}:=\max\{\frac{1}{2}\|W\|_{F}^{2}:\|W\|\leq 1\}=\frac{1}{2}\min\{m,n\}$, $D_{g}:=\max\{\frac{1}{2}\|Z\|_{F}^{2}:\|Z\|_{\infty}\leq\rho\}=\frac{1}{2}mn\rho^{2}$ and $(X^{*},Y^{*})$ and $(X^{*}(\sigma),Y^{*}(\sigma))$ be optimal solution to problems (1.6) and (5.3), respectively. Note that

Using the inequalities in (5.4) and (5.5) and the facts that $(X(\sigma),Y(\sigma))$ is an $\epsilon/2$-optimal solution to (5.3) and $\sigma\max\{D_{f},D_{g}\}=\frac{\epsilon}{4}$, we have

Thus, according to Theorem 19, to find an $\epsilon$-optimal solution to (1.6), we need to find an $\epsilon/2$-optimal solution to (5.3) with $\sigma=\frac{\epsilon}{2\max\{\min\{m,n\},mn\rho^{2}\}}$. We can either apply Algorithms 3 and 6 to solve (5.3), or apply Algorithms 4 and 7 to solve (5.3) with only one functions (say $f(x)$) smoothed. The iteration complexity results in Theorems 4 and 13 hold since the gradients of $f_{\sigma}$ is Lipschitz continuous. However, the numbers of iterations needed by ALM and FALM to obtain an $\epsilon$-optimal solution to (1.6) become $O(1/\epsilon^{2})$ and $O(1/\epsilon)$, respectively, due to the fact that the Lipschitz constant $L(f_{\sigma})=1/\sigma=\frac{2\max\{\min\{m,n\},mn\rho^{2}\}}{\epsilon}=O(1/\epsilon)$.

The two subproblems at iteration $k$ of Algorithm 3 when applied to (5.3) reduce to

The first-order optimality conditions for (5.6) are:

where $W_{\sigma}(X)$ and $Z_{\sigma}(Y)$ are defined in (5.2) and (4.3). It is easy to check that

satisfies (5.8), where $U\mbox{Diag}(\gamma)V^{\top}$ is the SVD of the matrix $\mu Z_{\sigma}(Y^{k})-Y^{k}+M$. Thus, solving the subproblem (5.6) corresponds to an SVD. If we define $B:=\mu W_{\sigma}(X^{k+1})-X^{k+1}+M,$ it is easy to verify that

satisfies the first-order optimality conditions for (5.7): $-W_{\sigma}(X^{k+1})+\frac{1}{\mu}(X^{k+1}+Y-M)+Z_{\sigma}(Y)=0.$ Thus, solving the subproblem (5.7) can be done very cheaply. The two subproblems at the $k$-th iteration of Algorithm 6 can be done in the same way and the main computational effort in each iteration of both ALM and FALM corresponds to an SVD.

2 RPCA with Missing Data

In some applications of RPCA, some of the entries of $M$ in (1.6) may be missing (e.g., in low-rank matrix completion problems where the matrix is corrupted by noise). Let $\Omega$ be the index set of the entries of $M$ that are observable and define the projection operator $\mathcal{P}_{\Omega}$ as: $(\mathcal{P}_{\Omega}(X))_{ij}=X_{ij}$, if $(i,j)\in\Omega$ and $(\mathcal{P}_{\Omega}(X))_{ij}=0$ otherwise. It has been shown under some randomness hypotheses that the low rank $\bar{X}$ and sparse $\bar{Y}$ can be recovered with high probability by solving (see Theorem 1.2 in ),

To solve (5.11) by ALM or FALM, we need to transform it into the form of (1.6). For this we have

$(\bar{X},\mathcal{P}_{\Omega}(\bar{Y}))$ is an optimal solution to (5.11) if

Suppose $(X^{*},Y^{*})$ is an optimal solution to (5.11). We claim that $Y^{*}_{ij}=0,\forall(i,j)\notin\Omega$. Otherwise, $(X^{*},\mathcal{P}_{\Omega}(Y^{*}))$ is feasible to (5.11) and has a strictly smaller objective function value than $(X^{*},Y^{*})$, which contradicts the optimality of $(X^{*},Y^{*})$. Thus, $\|\mathcal{P}_{\Omega}(Y^{*})\|_{1}=\|Y^{*}\|_{1}$. Now suppose that $(\bar{X},\mathcal{P}_{\Omega}(\bar{Y}))$ is not optimal to (5.11); then we have

which contradicts the optimality of $(\bar{X},\bar{Y})$ to (5.12). Therefore, $(\bar{X},\mathcal{P}_{\Omega}(\bar{Y}))$ is optimal to (5.11). ∎

The only differences between (1.6) and (5.12) lie in that the matrix $M$ is replaced by $\mathcal{P}_{\Omega}(M)$ and $g(Y)=\rho\|Y\|_{1}$ is replaced by $\rho\|\mathcal{P}_{\Omega}(Y)\|_{1}$. A smoothed approximation $g_{\sigma}(Y)$ to $g(Y):=\rho\|\mathcal{P}_{\Omega}(Y)\|_{1}$ is given by

According to Theorem 1 in , $\nabla g_{\sigma}(Y)$ is Lipschitz continuous with $L_{\sigma}(g)=1/\sigma$. Thus the convergence and iteration complexity results in Theorems 4 and 13 apply. The only changes in Algorithms 3 and 4 and Algorithms 6 and 7 are: replacing $M$ by $\mathcal{P}_{\Omega}(M)$ and computing $Y^{k+1}$ using (5.10) with $B$ is replaced by $\mathcal{P}_{\Omega}(B)$.

3 Numerical Results on RPCA Problems

In this section, we report numerical results obtained using the ALM method to solve RPCA problems with both complete and incomplete data matrices $M$. We compare the performance of ALM with the exact ADM (EADM) and the inexact ADM (IADM) methods in . The MATLAB codes of EADM and IADM were downloaded from $http://watt.csl.illinois.edu/\sim perceive/matrix-rank/sample\_code.html$ and their default settings were used. To further accelerate ALM, we adopted the continuation strategy used in EADM and IADM. Specifically, we set $\mu_{k+1}:=\max\{\bar{\mu},\eta\mu_{k}\}$, where $\mu_{0}=\|M\|/1.25,\bar{\mu}=10^{-6}$ and $\eta=2/3$ in our numerical experiments. Although in some iterations this violates the requirement $\mu\leq\min\{\frac{1}{L(f_{\sigma})},\frac{1}{L(g_{\sigma})}\}$ in Corollaries 5 and 14, we see from our numerical results reported below that ALM and FALM still work very well. We also found that by adopting this updating rule for $\mu$, there was not much difference between the performance of ALM and that of FALM. So we only compare ALM with EADM and IADM. As in Section 4, since we applied ALM to a smoothed problem, we implemented ALM as SADAL. The initial point in ALM was set to $(X^{0},Y^{0})=(M,\mathbf{0})$ and the initial Lagrange multiplier was set to $\Lambda^{0}=-\nabla g_{\sigma}(Y^{0})$. We set the smoothness parameter $\sigma=10^{-6}$. Solving subproblem (5.6) requires computing an SVD (see (5.9)). However, we do not have to compute the whole SVD, as only the singular values that are larger than the threshold $\tau=\frac{\mu\gamma}{\max\{\gamma,\mu+\sigma\}}$ and the corresponding singular vectors are needed. We therefore use PROPACK , which is also used in EADM and IADM, to compute these singular values and corresponding singular vectors. To use PROPACK, one has to specify the number of leading singular values (denoted by $sv_{k}$) to be computed at iteration $k$. We here adopt the strategy suggested in for EADM and IADM. This strategy starts with $sv_{0}=100$ and updates $sv_{k}$ via:

where $d=\min\{m,n\}$ and $svp_{k}$ is the number of singular values that are larger than the threshold $\tau$.

In all our experiments $\rho$ was chosen equal to $1/\sqrt{m}$. We stopped ALM, EADM and IADM when the relative infeasibility was less than $10^{-7}$, i.e., $\|X+Y-M\|_{F}<10^{-7}\|M\|_{F}$.

Extracting the almost still background from a sequence frames of video is a basic task in video surveillance. This problem is difficult due to the presence of moving foregrounds in the video. Interestingly, as shown in , this problem can be formulated as a RPCA problem (1.6). By stacking the columns of each frame into a long vector, we get a matrix $M$ whose columns correspond to the sequence of frames of the video. This matrix $M$ can be decomposed into the sum of two matrices $M:=\bar{X}+\bar{Y}$. The matrix $\bar{X}$, which represents the background in the frames, should be of low rank due to the correlation between frames. The matrix $\bar{Y}$, which represents the moving objects in the foreground in the frames, should be sparse since these objects usually occupy a small portion of each frame. We apply ALM to solve (1.6) for two videos introduced in .

3.2 Random Matrix Completion Problems with Grossly Corrupted Data

4 Sparse Inverse Covariance Selection

In ALM method was successfully applied to the Sparse Inverse Covariance Selection problem:

where $f(X)=-\log\det(X)+\langle S,X\rangle$ and $g(X)=\rho\|X\|_{1}$.

Note that in our case $f(X)$ does not have Lipschitz continuous gradient in general. Moreover, $f(X)$ is only defined for positive definite matrices while $g(X)$ is defined everywhere. These properties of the objective function make the SICS problem especially challenging for optimization methods. Nevertheless, we can still apply Algorithm 4 and obtain the complexity bound in Theorem 4 as follows. As proved in , the optimal solution $X^{*}$ of (5.18) satisfies $X\succeq\alpha I$, where $\alpha=\frac{1}{\|S\|+n\rho},$ (see Proposition 3.1 in ). Therefore, the SICS problem (5.18) can be formulated as:

where $\mathcal{C}:=\{X\in S^{n}:X\succeq\frac{\alpha}{2}I\}$. We can apply Algorithm 4 and Theorem 4 as per Remark 8. The difficulty arises, however, when performing minimization in $Y$ (Step 5 of Algorithm 4) with the constraint $Y\in\mathcal{C}$. Without this constraint, the minimization is obtained by a matrix shrinkage operation. However, the problem becomes harder to solve with this additional constraint. Minimization in $X$ (Step 3 of Algorithm 4) with or without the constraint $X\in\mathcal{C}$ is accomplished by performing an SVD of the current iterate $Y^{k}$. Hence the constraint can be easily imposed. Also note that once the SVD is computed both $\nabla f(X^{k+1})$ and $\nabla f(Y^{k})$ are readily available (see for details). This implies that either skipping or nonskipping iterations of Algorithm 4 can be performed at the same cost as one ISTA iteration.

Instead of imposing constraint $Y\in\mathcal{C}$ in Step 5 of Algorithm 4 we can obtain feasible solutions by a line search on $\mu$. We know that the constraint $X\succeq\frac{\alpha}{2}I$ is not tight at the solution. Hence if we start the algorithm with $X\succeq\alpha I$ and restrict the step size $\mu$ to be sufficiently small then the iterates of the method will remain in $\mathcal{C}$. Similarly, one can apply ISTA with small steps to remain in $\mathcal{C}$. Note however, that the bound on the Lipschitz constant of the gradient of $f(X)$ is $1/\alpha^{2}$ and hence can be very large. It is not practical to restrict $\mu$ in the algorithm to be smaller than $\alpha^{2}$, since $\mu$ determines the step size at each iteration. The advantage of ALM methods over ISTA in this case is that as soon as the $Y\in\mathcal{C}$ is relaxed ISTA can no longer be applied, while ALM/SADAL can be applied and indeed works very well. The theory in this case only applies once certain proximity to the optimal solution has been reached. But as shown in , the SADAL method is computationally superior to other state-of-the-art methods for SICS.

We have also applied the FALM method to the SICS problem, but we have not observed any advantage over ALM for this particular application.

Conclusion

In this paper, we proposed both basic and accelerated versions of alternating linearization methods for minimizing the sum of two convex functions. Our basic methods require at most $O(1/\epsilon)$ iterations to obtain an $\epsilon$-optimal solution, while our accelerated methods require at most $O(1/\sqrt{\epsilon})$ iterations with only a small additional amount of computational effort at each iteration. Numerical results on image deblurring, background extraction from surveillance video and matrix completion with grossly corrupted data are reported. These results demonstrate the efficiency and the practical potential of our algorithms.

Acknowledgement

We would like to thank Dr. Zaiwen Wen for insightful discussions on the topic of this paper.

References