Structured Nonconvex and Nonsmooth Optimization: Algorithms and Iteration Complexity Analysis

Introduction

In this paper, we consider the following nonconvex and nonsmooth optimization problem with multiple block variables:

where $S$ is a convex and compact set. In this paper, we propose several first-order algorithms for computing an $\epsilon$-stationary solution (to be defined later) for (1.1) and (1.2), and analyze their iteration complexities. Throughout, we assume the following condition.

The sets of the stationary solutions for (1.1) and (1.2) are non-empty.

where $\mathcal{C}$ is the core tensor that has a smaller size than $\mathcal{Z}$, and $X_{i}$ are matrices with appropriate sizes, $i=1,\ldots,d$. In fact, the “low-rank” tensor in the above model corresponds to the tensor with a small core; however a recent work demonstrates that the CP-rank of the core regardless of its size could be as large as the original tensor. Therefore, if one wants to find the low CP-rank decomposition, then the following model is preferred:

where $\mathcal{T}(x_{1},x_{2},\cdots,x_{d})=\sum_{i_{1},\dots,i_{d}}\mathcal{T}_{i_{1},\dots,i_{d}}(x_{1})_{i_{1}}\cdots(x_{d})_{i_{d}}$.

The convergence and iteration complexity for various nonconvex and nonsmooth optimization problems have recently attracted considerable research attention; see e.g. . In this paper, we study several solution methods that use only the first-order information of the objective function, including a generalized conditional gradient method, variants of alternating direction method of multipliers, and a proximal block coordinate descent method, for solving (1.1) and (1.2). Specifically, we apply a generalized conditional gradient (GCG) method to solve (1.2). We prove that the GCG can find an $\epsilon$-stationary solution for (1.2) in $O(\epsilon^{-q})$ iterations under certain mild conditions, where $q$ is a parameter in the Hölder condition that characterizes the degree of smoothness for $f$. In other words, the convergence rate of the algorithm depends on the degree of “smoothness” of the objective function. It should be noted that a similar iteration bound that depends on the parameter $q$ was reported for convex problems , and for general nonconvex problem, analyzed the convergence results, but there was no iteration complexity result. Furthermore, we show that if $f$ is concave, then GCG finds an $\epsilon$-stationary solution for (1.2) in $O(1/\epsilon)$ iterations. For the affinely constrained problem (1.1), we propose two algorithms (called proximal ADMM-g and proximal ADMM-m in this paper), both can be viewed as variants of the alternating direction method of multipliers (ADMM). Recently, there has been an emerging research interest on the ADMM for nonconvex problems (see, e.g., ). However, the results in only show that the iterates produced by the ADMM converge to a stationary solution without providing an iteration complexity analysis. Moreover, the objective function is required to satisfy the so-called Kurdyka-Łojasiewicz (KL) property to enable those convergence results. In , Hong, Luo and Razaviyayn analyzed the convergence of the ADMM for solving nonconvex consensus and sharing problems. Note that they also analyzed the iteration complexity of the ADMM for the consensus problem. However, they require the nonconvex part of the objective function to be smooth, and nonsmooth part to be convex. In contrast, $r_{i}$ in our model (1.1) can be nonconvex and nonsmooth at the same time. Moreover, we allow general constraints $x_{i}\in\mathcal{X}_{i},i=1,\ldots,N-1$, while the consensus problem in only allows such constraint for one block variable. A very recent work of Hong discussed the iteration complexity of an augmented Lagrangian method for finding an $\epsilon$-stationary solution for the following problem:

under the assumption that $f$ is differentiable. We will compare our results with in more details in Section 3.

Before proceeding, let us first summarize:

We provide definitions of $\epsilon$-stationary solution for (1.1) and (1.2) using the variational inequalities. For (1.1), our definition of the $\epsilon$-stationary solution allows each $r_{i}$ to be nonsmooth and nonconvex.

We study a generalized conditional gradient method with a suitable line search rule for solving (1.2). We assume that the gradient of $f$ satisfies a Hölder condition, and analyze its iteration complexity for obtaining an $\epsilon$-stationary solution for (1.2). After we released the first version of this paper, we noticed there are several recent works that study the iteration complexity of conditional gradient method for nonconvex problems. However, our results are different from these. For example, the convergence rate given in is worse than ours, and only consider smooth nonconvex problem with Lipschitz continuous gradient, but our results cover nonsmooth models.

We study two ADMM variants (proximal ADMM-g and proximal ADMM-m) for solving (1.1), and analyze their iteration complexities for obtaining an $\epsilon$-stationary solution for nonconvex problem (1.1). In addition, the setup and the assumptions of our model are different from other recent works. For instance, considers a two-block nonconvex problem with an identity coefficient matrix for one block variable in the linear constraint, and requires the coerciveness of the objective or the boundedness of the domain. assumes that the objective function is coercive over the feasible set and the nonsmooth objective is restricted prox-regular or piece-wise linear. While our algorithm assumes the gradient of the smooth part of the objective function is Lipschitz continuous and the nonsmooth part does not involve the last block variable, which is weaker than the assumptions on the objective functions in .

As an extension, we also show how to use proximal ADMM-g and proximal ADMM-m to find an $\epsilon$-stationary solution for (1.1) without assuming any condition on $A_{N}$.

When the affine constraints are absent in model (1.1), as a by-product, we demonstrate that the iteration complexity of proximal block coordinate descent (BCD) method with cyclic order can be obtained directly from that of proximal ADMM-g and proximal ADMM-m. Although gives an iteration complexity result of nonconvex BCD, it requires the KL property, and the complexity depends on a parameter in the KL condition, which is typically unknown.

Organization. The rest of this paper is organized as follows. In Section 2 we introduce the notion of $\epsilon$-stationary solution for (1.2) and apply a generalized conditional gradient method to solve (1.2) and analyze its iteration complexity for obtaining an $\epsilon$-stationary solution for (1.2). In Section 3 we give two definitions of $\epsilon$-stationarity for (1.1) under different settings and propose two ADMM variants that solve (1.1) and analyze their iteration complexities to reach an $\epsilon$-stationary solution for (1.1). In Section 4 we provide some extensions of the results in Section 3. In particular, we first show how to remove some of the conditions that we assume in Section 3, and then we apply a proximal BCD method to solve (1.1) without affine constraints and provide an iteration complexity analysis. In Section 5, we present numerical results to illustrate the practical efficiency of the proposed algorithms.

A generalized conditional gradient method

In this section, we study a GCG method for solving (1.2) and analyze its iteration complexity. The conditional gradient (CG) method, also known as the Frank-Wolfe method, was originally proposed in , and regained a lot of popularity recently due to its capability in solving large-scale problems (see, ). However, these works focus on solving convex problems. Bredies et al. proved the convergence of a generalized conditional gradient method for solving nonconvex problems in Hilbert space. In this section, by introducing a suitable line search rule, we provide an iteration complexity analysis for this algorithm.

We make the following assumption in this section regarding (1.2).

In (1.2), $r(x)$ is convex and nonsmooth, and the constraint set $S$ is convex and compact. Moreover, $f$ is differentiable and there exist some $p>1$ and $\rho>0$ such that

The above inequality (2.1) is also known as the Hölder condition and was used in other works on first-order algorithms (e.g., ). It can be shown that (2.1) holds for a variety of functions. For instance, (2.1) holds for any $p$ when $f$ is concave, and is valid for $p=2$ when $\nabla f$ is Lipschitz continuous.

For smooth unconstrained problem $\min_{x}f(x)$, it is natural to define the $\epsilon$-stationary solution using the criterion $\|\nabla f(x)\|_{2}\leq\epsilon.$ Nesterov and Cartis et al. showed that the gradient descent type methods with properly chosen step size need $O({1}/{\epsilon^{2}})$ iterations to find such a solution. Moreover, Cartis et al. constructed an example showing that the $O({1}/{\epsilon^{2}})$ iteration complexity is tight for the steepest descent type algorithm. However, the case for the constrained nonsmooth nonconvex optimization is subtler. There exist some works on how to define $\epsilon$-optimality condition for the local minimizers of various constrained nonconvex problems . Cartis et al. proposed an approximate measure for smooth problem with convex set constraint. discussed general nonsmooth nonconvex problem in Banach space by using the tool of limiting Fréchet $\epsilon$-subdifferential. showed that under certain conditions $\epsilon$-KKT solutions can converge to a stationary solution as $\epsilon\to 0$. Here the $\epsilon$-KKT solution is defined by relaxing the complimentary slackness and equilibrium equations of KKT conditions. Ghadimi et al. considered the following notion of $\epsilon$-stationary solution for (1.2):

where $\gamma>0$ and $V$ is a prox-function. They proposed a projected gradient algorithm to solve (1.2) and proved that it takes no more than $O({1}/{\epsilon^{2}})$ iterations to find an $x$ satisfying

Our definition of an $\epsilon$-stationary solution for (1.2) is as follows.

We call $x$ an $\epsilon$-stationary solution ($\epsilon\geq 0$) for (1.2) if the following holds:

If $\epsilon=0$, then $x$ is called a stationary solution for (1.2).

Observe that if $r(\cdot)$ is continuous then any cluster point of $\epsilon$-stationary solutions defined above is a stationary solution for (1.2) as $\epsilon\to 0$. Moreover, the stationarity condition is weaker than the usual KKT optimality condition. To see this, we first rewrite (1.2) as the following equivalent unconstrained problem

where $\iota_{S}(x)$ is the indicator function of $S$. Suppose that $x$ is any local minimizer of this problem and thus also a local minimizer of (1.2). Since $f$ is differentiable, $r$ and $\iota_{S}$ are convex, Fermat’s rule yields

which further implies that there exists some $z\in\partial r(x)$ such that

Using the convexity of $r(\cdot)$, it is equivalent to

Therefore, (2.6) is a necessary condition for local minimum of (1.2) as well. Furthermore, we claim that $\psi_{S}(x)\geq-\epsilon$ implies $\|P_{S}(x,\gamma)\|_{2}^{2}\leq\epsilon/\gamma$ with the prox-function $V(y,x)=\|y-x\|_{2}^{2}/2$. In fact, (2.2) guarantees that

for some $z\in\partial r(x^{+})$. By choosing $y=x$ in (2.7) one obtains

Therefore, if $\psi_{S}(x)\geq-\epsilon$, then $\|P_{S}(x,\gamma)\|_{2}^{2}\leq\frac{\epsilon}{\gamma}$ holds.

2 The algorithm

For given point $z$, we define an approximation of the objective function of (1.2) to be:

which is obtained by linearizing the smooth part (function $f$) of $\Phi$ in (1.2). Our GCG method for solving (1.2) is described in Algorithm 1, where $\rho$ and $p$ are from Assumption 2.1.

Note that here we assumed that solving the subproblem in Step 1 of Algorithm 1 is relatively easy. That is, we assumed the following assumption.

All subproblems in Step 1 of Algorithm 1 can be solved relatively easily.

Assumption 2.3 is quite common in conditional gradient method. For a list of functions $r$ and sets $S$ such that Assumption 2.3 is satisfied, see .

It is easy to see that the sequence $\{\Phi(x^{k})\}$ generated by GCG is monotonically nonincreasing , which implies that any cluster point of $\{x^{k}\}$ cannot be a strict local maximizer.

3 An iteration complexity analysis

Before we proceed to the main result on iteration complexity of GCG, we need the following lemma that gives a sufficient condition for an $\epsilon$-stationary solution for (1.2). This lemma is inspired by , and it indicates that if the progress gained by minimizing (2.9) is small, then $z$ must already be close to a stationary solution for (1.2).

Denoting $\Phi^{*}$ to be the optimal value of (1.2), we are now ready to give the main result of the iteration complexity of GCG (Algorithm 1) for obtaining an $\epsilon$-stationary solution for (1.2).

where the third inequality is due to the convexity of function $r$ and the fact that $x^{k+1}-x^{k}=\alpha_{k}(y^{k}-x^{k})$, and the last inequality is due to (2.1). Furthermore, (2.3) immediately yields

For any integer $K>0$, summing (2.11) over $k=0,1,\ldots,K-1$, yields

Finally, if $f$ is concave, then the iteration complexity can be improved as $O(1/\epsilon)$.

Suppose that $f$ is a concave function. If we set $\alpha_{k}=1$ for all $k$ in GCG (Algorithm 1), then it returns an $\epsilon$-stationary solution for (1.2) within $\left\lceil\frac{\Phi(x^{0})-\Phi^{*}}{\epsilon}\right\rceil$ iterations.

Proof. By setting $\alpha_{k}=1$ in Algorithm 1 we have $x^{k+1}=y^{k}$ for all $k$. Since $f$ is concave, it holds that

Variants of ADMM for solving nonconvex problems with affine constraints

In this section, we study two variants of the ADMM (Alternating Direction Method of Multipliers) for solving the general problem (1.1), and analyze their iteration complexities for obtaining an $\epsilon$-stationary solution (to be defined later) under certain conditions. Throughout this section, the following two assumptions regarding problem (1.1) are assumed.

and $r_{i}^{*}=\inf\limits_{x_{i}\in\mathcal{X}_{i}}\ \{r_{i}(x_{i})\}$ for $i=1,2,\ldots,N-1$.

To characterize the optimality conditions for (1.1) when $r_{i}$ is nonsmooth and nonconvex, we need to recall the notion of the generalized gradient (see, e.g., ).

(ii). $v$ is a general subgradient of $h$ at $\bar{x}$, written $v\in\partial h(\bar{x})$, if there exist sequences $\{x^{k}\}$ and $\{v^{k}\}$ such that $x^{k}\to\bar{x}$ with $h(x^{k})\to h(\bar{x})$, and $v^{k}\in\hat{\partial}h(x^{k})$ with $v^{k}\to v$ when $k\to\infty$.

The following proposition lists some well-known facts about the lower semi-continuous functions.

In our analysis, we frequently use the following identity that holds for any vectors $a,b,c,d$,

2 An ϵitalic-ϵ\epsilon-stationary solution for problem (1.1)

where ${g}^{*}_{i}$ is a general subgradient of $r_{i}$ at point $x_{i}^{*}$. If $\epsilon=0$, we call ${\left(x_{1}^{*},\cdots,x_{N}^{*}\right)}$ a stationary solution for (1.1).

where ${\partial}r_{i}(x_{i}^{*})$ is the general subgradient of $r_{i}$ at $x_{i}^{*}$, $i=1,2,\ldots,N-1$. If $\epsilon=0$, we call ${\left(x_{1}^{*},\cdots,x_{N}^{*}\right)}$ to be a stationary solution for (1.1).

The two settings of problem (1.1) considered in this section and their corresponding definitions of $\epsilon$-stationary solution, are summarized in Table 1.

A very recent work of Hong proposes a definition of an $\epsilon$-stationary solution for problem (1.4), and analyzes the iteration complexity of a proximal augmented Lagrangian method for obtaining such a solution. Specifically, $(x^{*},\lambda^{*})$ is called an $\epsilon$-stationary solution for (1.4) in if $Q(x^{*},\lambda^{*})\leq\epsilon$, where

and $\mathcal{L}_{\beta}(x,\lambda):=f(x)-\lambda^{\top}\left(Ax-b\right)+\frac{\beta}{2}\left\|Ax-b\right\|^{2}$ is the augmented Lagrangian function of (1.4). Note that assumes that $f$ is differentiable and has bounded gradient in (1.4). It is easy to show that an $\epsilon$-stationary solution in is equivalent to an $O(\sqrt{\epsilon})$-stationary solution for (1.1) according to Definition 3.6 with $r_{i}=0$ and $f$ being differentiable. Note that there is no set constraint in (1.4), and so the notion of the $\epsilon$-stationarity in is not applicable in the case of Definition 3.5.

Consider the $\epsilon$-stationary solution in Definition 3.6 applied to problem (1.4), i.e., one block variable and $r_{i}(x)=0$. Then $x^{*}$ is a $\gamma_{1}\sqrt{{\epsilon}}$-stationary solution in Definition 3.6, with Lagrange multiplier $\lambda^{*}$ and $\gamma_{1}=1/(\sqrt{2\beta^{2}\|A\|_{2}^{2}+3})$, implies $Q(x^{*},\lambda^{*})\leq\epsilon$. On the contrary, if $Q(x^{*},\lambda^{*})\leq\epsilon$, then $x^{*}$ is a $\gamma_{2}\sqrt{{\epsilon}}$-stationary solution from Definition 3.6 with Lagrange multiplier $\lambda^{*}$, where $\gamma_{2}=\sqrt{2(1+\beta^{2}\|A\|_{2}^{2})}$.

Proof. Suppose $x^{*}$ is a $\gamma_{1}\sqrt{{\epsilon}}$-stationary solution as defined in Definition 3.6. We have $\|\nabla f(x^{*})-A^{\top}\lambda^{*}\|\leq\gamma_{1}\sqrt{\epsilon}$ and $\|Ax^{*}-b\|\leq\gamma_{1}\sqrt{\epsilon}$, which implies that

On the other hand, if $Q(x^{*},\lambda^{*})\leq\epsilon$, then we have $\|\nabla f(x^{*})-A^{\top}\lambda^{*}+\beta A^{\top}(Ax^{*}-b)\|^{2}\leq\epsilon$ and $\|Ax^{*}-b\|^{2}\leq\epsilon$. Therefore,

The desired result then follows immediately. $\Box$

In the following, we introduce two variants of ADMM, to be called proximal ADMM-g and proximal ADMM-m, that solve (1.1) under some additional assumptions on $A_{N}$. In particular, proximal ADMM-g assumes $A_{N}=I$, and proximal ADMM-m assumes $A_{N}$ to have full row rank.

3 Proximal gradient-based ADMM (proximal ADMM-g)

Our proximal ADMM-g solves (1.1) under the condition that $A_{N}=I$. In this case, the problem reduces to a so-called sharing problem in the literature which has the following form

For applications of the sharing problem, see . Our proximal ADMM-g for solving (1.1) with $A_{N}=I$ is described in Algorithm 2. It can be seen from Algorithm 2 that proximal ADMM-g is based on the framework of augmented Lagrangian method, and can be viewed as a variant of the ADMM. The augmented Lagrangian function of (1.1) is defined as

where $\lambda$ is the Lagrange multiplier associated with the affine constraint, and $\beta>0$ is a penalty parameter. In each iteration, proximal ADMM-g minimizes the augmented Lagrangian function plus a proximal term for block variables $x_{1},\ldots,x_{N-1}$, with other variables being fixed; and then a gradient descent step is conducted for $x_{N}$, and finally the Lagrange multiplier $\lambda$ is updated. The interested readers are referred to for gradient-based ADMM and its various stochastic variants for convex optimization.

which guarantee the convergence rate of the algorithm as shown in Lemma 3.9 and Theorem 3.12.

Before presenting the main result on the iteration complexity of proximal ADMM-g, we need some lemmas.

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 2. The following inequality holds

Proof. Note that Steps 2 and 3 of Algorithm 2 yield that

We now define the following function, which will play a crucial role in our analysis:

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 2, where the parameters $\beta$ and $\gamma$ are taken according to (3.8) and (3.9) respectively. Then $\Psi_{G}(x_{1}^{k+1},\cdots,x_{N}^{k+1},\lambda^{k+1},x_{N}^{k})$ monotonically decreases over $k\geq 0$.

Proof. From Step 1 of Algorithm 2 it is easy to see that

where the inequality follows from (3.2) and (3.3). Moreover, the following equality holds trivially

Combining (3.13), (3.14), (3.15) and (3.10) yields that

It is easy to verify that when $\beta>\frac{18\sqrt{3}+6}{13}L$, then $\gamma$ defined as in (3.9) ensures that $\gamma>0$ and

Therefore, choosing $\beta>\max\left(\frac{18\sqrt{3}+6}{13}L,\;\max\limits_{i=1,2,\ldots,N-1}\frac{6L^{2}}{\sigma_{\min}(H_{i})}\right)$ and $\gamma$ as in (3.9) guarantees that $\Psi_{G}(x_{1}^{k+1},\cdots,x_{N}^{k+1},\lambda^{k+1},x_{N}^{k})$ monotonically decreases over $k\geq 0$. In fact, (3.17) can be verified as follows. By denoting $z=\beta-\frac{1}{\gamma}$, (3.17) is equivalent to

which holds when $\beta>\frac{18\sqrt{3}+6}{13}L$ and $\frac{-\beta-\sqrt{13\beta^{2}-12\beta L-72L^{2}}}{12}<z<\frac{-\beta+\sqrt{13\beta^{2}-12\beta L-72L^{2}}}{12}$, i.e.,

which holds when $\gamma$ is chosen as in (3.9). $\Box$

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 2. Under the same conditions as in Lemma 3.10, for any $k\geq 0$, we have

where $r_{i}^{*}$ and $f^{*}$ are defined in Assumption 3.2.

where the first inequality follows from (3.2), and the second inequality is due to $\beta\geq 3L/2$. The desired result follows from the definition of $\Psi_{G}$ in (3.12). $\Box$

Now we are ready to give the iteration complexity of Algorithm 2 for finding an $\epsilon$-stationary solution of (1.1).

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 2. Furthermore, suppose that $\beta$ satisfies (3.8) and $\gamma$ satisfies (3.9). Denote

Then to get an $\epsilon$-stationary solution, the number of iterations that the algorithm runs can be upper bounded by:

and we can further identify one iteration $\hat{k}\in\mathop{\rm argmin}\limits_{2\leq k\leq K+1}\sum_{i=1}^{N}\left(\|x_{i}^{k}-x_{i}^{k+1}\|^{2}+\|x_{i}^{k-1}-x_{i}^{k}\|^{2}\right)$ such that $(x_{1}^{\hat{k}},\cdots,x_{N}^{\hat{k}})$ is an $\epsilon$-stationary solution for optimization problem (1.1) with Lagrange multiplier $\lambda^{\hat{k}}$ and $A_{N}=I$, for Settings 1 and 2 respectively.

By summing (3.3) over $k=1,\ldots,K$, we obtain that

where $\tau$ is defined in (3.18). By invoking Lemmas 3.10 and 3.11, we get

We now derive upper bounds on the terms in (3.5) and (3.6) through $\theta_{k}$. Note that (3.11) implies that

From Step 3 of Algorithm 2 and (3.10) it is easy to see that

We now derive upper bounds on the terms in (3.4) and (3.7) under the two settings in Table 1, respectively.

By combining (3.3), (3.22) and (3.27) we conclude that Algorithm 2 returns an $\epsilon$-stationary solution for (1.1) according to Definition 3.6 under the conditions of Setting 2 in Table 1.

where $g_{i}\in\partial r_{i}(x_{i}^{k+1})$ is a general subgradient of $r_{i}$ at $x_{i}^{k+1}$. By combining (3.3), (3.22) and (3.27) we conclude that Algorithm 2 returns an $\epsilon$-stationary solution for (1.1) according to Definition 3.5 under the conditions of Setting 1 in Table 1.

Note that the potential function $\Psi_{G}$ defined in (3.12) is related to the augmented Lagrangian function. The augmented Lagrangian function has been used as a potential function in analyzing the convergence of nonconvex splitting and ADMM methods in . See for a more detailed discussion on this.

In Step 1 of Algorithm 2, we can also replace the function

so that the subproblem can be solved by computing the proximal mappings of $r_{i}$, with some properly chosen matrix $H_{i}$ for $i=1,\ldots,N-1$, and the same iteration bound still holds.

4 Proximal majorization ADMM (proximal ADMM-m)

Our proximal ADMM-m solves (1.1) under the condition that $A_{N}$ has full row rank. In this section, we use $\sigma_{N}$ to denote the smallest eigenvalue of $A_{N}A_{N}^{\top}$. Note that $\sigma_{N}>0$ because $A_{N}$ has full row rank. Our proximal ADMM-m can be described as follows

In Algorithm 3, $U(x_{1},\cdots,x_{N-1},x_{N},\lambda,\bar{x})$ is defined as

to guarantee the convergence rate of the algorithm shown in Lemma 3.16 and Theorem 3.18.

It is worth noting that the proximal ADMM-m and proximal ADMM-g differ only in Step 2: Step 2 of proximal ADMM-g takes a gradient step of the augmented Lagrangian function with respect to $x_{N}$, while Step 2 of proximal ADMM-m requires to minimize a quadratic function of $x_{N}$.

We provide some lemmas that are useful in analyzing the iteration complexity of proximal ADMM-m for solving (1.1).

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 3. The following inequality holds

Proof. From the optimality conditions of Step 2 of Algorithm 3, we have

where the second equality is due to Step 3 of Algorithm 3. Therefore, we have

We define the following function that will be used in the analysis of proximal ADMM-m:

Similar to the function used in proximal ADMM-g, we can prove the monotonicity and boundedness of function $\Psi_{L}$.

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 3, where$\beta$ is chosen according to (3.30). Then $\Psi_{L}(x^{k+1},\cdots,x_{N}^{k+1},\lambda^{k+1},x_{N}^{k})$ monotonically decreases over $k>0$.

Proof. By Step 1 of Algorithm 3 one observes that

where the first inequality is due to (3.2) and (3.3). Moreover, from (3.31) we have

Combining (3.32), (3.40) and (3.4) yields that

where the second inequality is due to (3.30). This completes the proof. $\Box$

The following lemma shows that the function $\Psi_{L}$ is lower bounded.

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by Algorithm 3. Under the same conditions as in Lemma 3.16, the sequence $\{\Psi_{L}(x^{k+1},\cdots,x_{N}^{k+1},\lambda^{k+1},x_{N}^{k})\}$ is bounded from below.

Proof. From Step 3 of Algorithm 3 we have

where the third equality follows from (3.3). Summing this inequality over $k=0,1,\ldots,K-1$ for any integer $K\geq 1$ yields that

Lemma 3.16 stipulates that $\{\Psi_{L}(x_{1}^{k+1},\cdots,x_{N}^{k+1},\lambda^{k+1},x_{N}^{k})\}$ is a monotonically decreasing sequence; the above inequality thus further implies that the entire sequence is bounded from below. $\Box$

We are now ready to give the iteration complexity of proximal ADMM-m, whose proof is similar to that of Theorem 3.12.

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},\lambda^{k})\}$ is generated by proximal ADMM-m (Algorithm 3), and $\beta$ satisfies (3.30). Denote

Then to get an $\epsilon$-stationary solution, the number of iterations that the algorithm runs can be upper bounded by:

and we can further identify one iteration $\hat{k}\in\mathop{\rm argmin}\limits_{2\leq k\leq K+1}\sum_{i=1}^{N}\left(\|x_{i}^{k}-x_{i}^{k+1}\|^{2}+\|x_{i}^{k-1}-x_{i}^{k}\|^{2}\right)$, such that $(x_{1}^{\hat{k}},\cdots,x_{N}^{\hat{k}})$ is an $\epsilon$-stationary solution for (1.1) with Lagrange multiplier $\lambda^{\hat{k}}$ and $A_{N}$ being full row rank, for Settings 1 and 2 respectively.

Proof. By summing (3.4) over $k=1,\ldots,K$, we obtain that

where $\tau$ is defined in (3.44). From Lemma 3.17 we know that there exists a constant $\Psi_{L}^{*}$ such that $\Psi(x_{1}^{k+1},\cdots,x_{N}^{k+1},\lambda^{k+1},x_{N}^{k})\geq\Psi_{L}^{*}$ holds for any $k\geq 1$. Therefore,

where $\theta_{k}$ is defined in (3.20), i.e., for $K$ defined as in (3.45), $\theta_{\hat{k}}=O(\epsilon^{2})$.

We now give upper bounds to the terms in (3.5) and (3.6) through $\theta_{k}$. Note that Step 2 of Algorithm 3 implies that

By Step 3 of Algorithm 3 and (3.31) we have

The remaining proof is to give upper bounds to the terms in (3.4) and (3.7). Since the proof steps are almost the same as Theorem 3.12, we shall only provide the key inequalities below.

Setting 2. Under conditions in Setting 2 in Table 1, the inequality (3.3) becomes

By combining (3.50), (3.48) and (3.4) we conclude that Algorithm 3 returns an $\epsilon$-stationary solution for (1.1) according to Definition 3.6 under the conditions of Setting 2 in Table 1.

Setting 1. Under conditions in Setting 1 in Table 1, the inequality (3.3) becomes

By combining (3.4), (3.48) and (3.4) we conclude that Algorithm 3 returns an $\epsilon$-stationary solution for (1.1) according to Definition 3.5 under the conditions of Setting 1 in Table 1.

In Step 1 of Algorithm 3, we can replace the function $f(x_{1}^{k+1},\cdots,x_{i-1}^{k+1},x_{i},x_{i+1}^{k},\cdots,x_{N}^{k})$ by its linearization

Under the same conditions as in Remark 3.14, the same iteration bound follows by slightly modifying the analysis above.

Extensions

It is noted that in (1.1), we have some restrictions on the last block variable $x_{N}$, i.e., $r_{N}\equiv 0$ and $A_{N}=I$ or is full row rank. In this subsection, we show how to remove these restrictions and consider the more general problem

Before proceeding, we make the following assumption on (4.1).

holds uniformly for all $(x_{1},\cdots,x_{N})\in S$, where $A=[A_{1},\ldots,A_{N}]$.

We introduce the following problem that is closely related to (4.1):

where $\epsilon>0$ is the target tolerance, and $\mu(\epsilon)$ is a function of $\epsilon$ which will be specified later. Now, proximal ADMM-m is ready to be used for solving (4.2) because $A_{N+1}=I$ and $y$ is unconstrained. We have the following iteration complexity result for proximal ADMM-m to obtain an $\epsilon$-stationary solution of (4.1); proximal ADMM-g can be analyzed similarly.

Consider problem (4.1) under Setting 2 in Table 1. Suppose that Assumption 4.1 holds, and the objective in (4.1), i.e., $f+\sum_{i=1}^{N}r_{i}$, has a bounded level set. Furthermore, suppose that $f$ has a Lipschitz continuous gradient with Lipschitz constant $L$, and $A$ is of full row rank. Now let the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k},y^{k},\lambda^{k})\}$ be generated by proximal ADMM-m for solving (4.2) with initial iterates $y^{0}=\lambda^{0}=0$, and $(x_{1}^{0},\cdots,x_{N}^{0})$ such that $\sum_{i=1}^{N}A_{i}x_{i}^{0}=b$. Assume that the target tolerance $\epsilon$ satisfies

Then in no more than $O(1/\epsilon^{4})$ iterations we will reach an iterate $(x_{1}^{{\hat{K}}+1},\cdots,x_{N}^{{\hat{K}}+1},y^{{\hat{K}}+1})$ that is an $\epsilon$-stationary solution for (4.2) with Lagrange multiplier $\lambda^{{\hat{K}}+1}$. Moreover, $(x_{1}^{{\hat{K}}+1},\cdots,x_{N}^{{\hat{K}}+1})$ is an $\epsilon$-stationary solution for (4.1) with Lagrange multiplier $\lambda^{\hat{K}+1}$.

Proof. Denote the penalty parameter as $\beta(\epsilon)$. The augmented Lagrangian function of (4.2) is given by

From (4.3) we have $\mu(\epsilon)>L$. This implies that the Lipschitz constant of the smooth part of the objective of (4.2) is equal to $\mu(\epsilon)$. Then from the optimality conditions of Step 2 of Algorithm 3, we have $\mu(\epsilon)y^{k}=\lambda^{k},\forall k\geq 1$.

Similar to Lemma 3.16, we can prove that $\mathcal{L}_{\beta(\epsilon)}(x_{1}^{k},\ldots,x_{N}^{k},y^{k},\lambda^{k})$ monotonically decreases. Specifically, since $\mu(\epsilon)y^{k}=\lambda^{k}$, combining (3.32), (3.40) and the equality in (3.4) yields,

where the last inequality is due to (4.4).

Similar to Lemma 3.17, we can prove that $\mathcal{L}_{\beta(\epsilon)}(x_{1}^{k},\cdots,x_{N}^{k},y^{k},\lambda^{k})$ is bounded from below, i.e., the exists a constant $\mathcal{L}^{*}=f^{*}+\sum_{i=1}^{N}r_{i}^{*}$ such that

Actually the following inequalities lead to the above fact:

where the second equality is from $\mu(\epsilon)y^{k}=\lambda^{k}$, and the last inequality is due to (4.4). Moreover, denote $\mathcal{L}^{0}\equiv\mathcal{L}_{\beta(\epsilon)}(x_{1}^{0},\cdots,x_{N}^{0},y^{0},\lambda^{0})$, which is a constant independent of $\epsilon$.

Furthermore, for any integer $K\geq 1$, summing (4.5) over $k=0,\ldots,K$ yields

where $\theta_{k}:=\sum_{i=1}^{N}\|x_{i}^{k}-x_{i}^{k+1}\|^{2}+\|y^{k}-y^{k+1}\|^{2}$. Note that (4.7) and (4.6) imply that

Similar to (3.3), it can be shown that for $i=1,\ldots,N$,

Set $K=1/\epsilon^{4}$ and denote $\hat{K}=\mathop{\rm argmin}_{0\leq k\leq K}\theta_{k}$. Then we know $\theta_{\hat{K}}=O(\epsilon^{4})$. As a result,

Note that (4.6) also implies that $f(x_{1}^{k},\cdots,x_{N}^{k})+\sum_{i=1}^{N}r_{i}(x_{i}^{k})$ is upper-bounded by a constant. Thus, from the assumption that the level set of the objective is bounded, we know $(x_{1}^{k},\cdots,x_{N}^{k})$ is bounded. Then Assumption 4.1 implies that $\lambda^{k}$ bounded, which results in $\|y^{k}\|=O(\epsilon)$. Therefore, from (4.12) we have

which combining with (4.11) yields that $(x_{1}^{\hat{K}+1},\cdots,x_{N}^{\hat{K}+1})$ is an $\epsilon$-stationary solution for (4.1) with Lagrange multiplier $\lambda^{\hat{K}+1}$, according to Definition 3.6. $\Box$

Without Assumption 4.1, we can still provide an iteration complexity of proximal ADMM-m, but the complexity bound is worse than $O(1/\epsilon^{4})$. To see this, note that because $\mathcal{L}_{\beta(\epsilon)}(x_{1}^{k},\cdots,x_{N}^{k},y^{k},\lambda^{k})$ monotonically decreases, the first inequality in (4.6) implies that

Therefore, by setting $K=1/\epsilon^{6}$, $\mu(\epsilon)=1/\epsilon^{2}$ and $\beta(\epsilon)=3/\epsilon^{2}$ instead of (4.4), and combining (4.11) and (4.13), we conclude that $(x_{1}^{\hat{K}+1},\cdots,x_{N}^{\hat{K}+1})$ is an $\epsilon$-stationary solution for (4.1) with Lagrange multiplier $\lambda^{\hat{K}+1}$, according to Definition 3.6.

2 Proximal BCD (Block Coordinate Descent)

In this section, we apply a proximal block coordinate descent method to solve the following variant of (1.1) and present its iteration complexity:

Similar to the settings in Table 1, depending on the properties of $r_{i}$ and $\mathcal{X}_{i}$, the $\epsilon$-stationary solution for (4.14) is as follows.

$(x_{1}^{*},\ldots,x_{N}^{*},\lambda^{*})$ is called an $\epsilon$-stationary solution for (4.14), if

$r_{i}$ is Lipschitz continuous, $\mathcal{X}_{i}$ is convex and compact, and for any $x_{i}\in\mathcal{X}_{i}$, $i=1,\ldots,N$, it holds that ($g_{i}={\partial}r_{i}(x_{i}^{*})$ denotes a generalized subgradient of $r_{i}$)

It is easy to see that applying proximal ADMM-g to solve (4.16) (with $x_{N+1}$ being the last block variable) reduces exactly to Algorithm 4. Hence, we have the following iteration complexity result of Algorithm 4 for obtaining an $\epsilon$-stationary solution of (4.14).

Suppose the sequence $\{(x_{1}^{k},\cdots,x_{N}^{k})\}$ is generated by proximal BCD (Algorithm 4). Denote

with $\tau$ being defined in (3.18), and $\hat{K}:=\min\limits_{1\leq k\leq K}\sum_{i=1}^{N}\left(\|x_{i}^{k}-x_{i}^{k+1}\|^{2}\right)$, we have that $(x_{1}^{\hat{K}},\cdots,x_{N}^{\hat{K}})$ is an $\epsilon$-stationary solution for problem (4.14).

Proof. Note that $A_{1}=\cdots=A_{N}=0$ and $A_{N+1}=I$ in problem (4.16). By applying proximal ADMM-g with $\beta>\max\left\{18L,\;\max\limits_{1\leq i\leq N}\left\{\frac{6L^{2}}{\sigma_{\min}(H_{i})}\right\}\right\}$, Theorem 3.12 holds. In particular, (3.3) and (3.3) are valid in different settings with $\beta\sqrt{N}\max\limits_{i+1\leq j\leq N+1}\left[\|A_{j}\|_{2}\right]\|A_{i}\|_{2}=0$ for $i=1,\ldots,N$, which leads to the choices of $\kappa_{5}$ and $\kappa_{6}$ in the above. Moreover, we do not need to consider the optimality with respect to $x_{N+1}$ and the violation of the affine constraints, thus $\kappa_{1}$ and $\kappa_{2}$ in Theorem 3.12 are excluded in the expression of $K$, and the conclusion follows. $\Box$

Numerical Experiments

The following identities are useful for our presentation later:

where ${Z}_{(i)}$ stands for the mode-$i$ unfolding of tensor $\mathcal{Z}$ and $\odot$ stands for the Khatri-Rao product of matrices.

Note that there are six block variables in (5.1), and we choose $\mathcal{B}$ as the last block variable. A typical iteration of proximal ADMM-g for solving (5.1) can be described as follows (we chose $H_{i}=\delta_{i}I$, with $\delta_{i}>0,i=1,\ldots,5$):

where $\circ$ is the matrix Hadamard product and $\mathcal{S}$ stands for the soft shrinkage operator. The updates in proximal ADMM-m are almost the same as proximal ADMM-g except ${B}_{(1)}$ is updated as

On the other hand, note that (5.1) can be equivalently written as

which can be solved by the classical BCD method as well as our proximal BCD (Algorithm 4).

In the following we shall compare the numerical performance of BCD, proximal BCD, proximal ADMM-g and proximal ADMM-m for solving (5.1). We let $\alpha=2/\max\{\sqrt{I_{1}},\sqrt{I_{2}},\sqrt{I_{3}}\}$ and $\alpha_{\mathcal{N}}=1$ in model (5.1). We apply proximal ADMM-g and proximal ADMM-m to solve (5.1), and apply BCD and proximal BCD to solve (5.2). In all the four algorithms we set the maximum iteration number to be $2000$, and the algorithms are terminated either when the maximum iteration number is reached or when $\theta_{k}$ as defined in (3.20) is less than $10^{-6}$. The parameters used in the two ADMM variants are specified in Table 2.

In the experiment, we randomly generate $20$ instances for fixed tensor dimension and CP-rank. Suppose the low-rank part $\mathcal{Z}^{0}$ is of rank $R_{CP}$. It is generated by

where vectors $a^{i,r}$ are generated from standard Gaussian distribution for $i=1,2,3$, $r=1,\dots,R_{CP}$. Moreover, a sparse tensor $\mathcal{E}^{0}$ is generated with cardinality of $0.001\cdot I_{1}I_{2}I_{3}$ such that each nonzero component follows from standard Gaussian distribution. Finally, we generate noise $\mathcal{B}^{0}=0.001*\hat{\mathcal{B}}$, where $\hat{\mathcal{B}}$ is a Gaussian tensor. Then we set $\mathcal{T}=\mathcal{Z}^{0}+\mathcal{E}^{0}+\mathcal{B}^{0}$ as the observed data in (5.1). A proper initial guess $R$ of the true rank $R_{CP}$ is essential for the success of our algorithms. We can borrow the strategy in matrix completion , and start from a large $R$ ($R\geq R_{CP}$) and decrease it aggressively once a dramatic change in the recovered tensor $\mathcal{Z}$ is observed. We report the average performance of 20 instances of the four algorithms with initial guess $R=R_{CP}$, $R=R_{CP}+1$ and $R=R_{CP}+\lceil 0.2*R_{CP}\rceil$ in Tables 3, 4 and 5, respectively.

In Tables 3, 4 and 5, “Err.” denotes the averaged relative error $\frac{\|\mathcal{Z}^{*}-\mathcal{Z}^{0}\|_{F}}{\|\mathcal{Z}^{0}\|_{F}}$ of the low-rank tensor over 20 instances, where $\mathcal{Z}^{*}$ is the solution returned by the corresponding algorithm; “Iter.” denotes the averaged number of iterations over 20 instances; “Num” records the number of solutions (out of 20 instances) that have relative error less than $0.01$.

Tables 3, 4 and 5 suggest that BCD mostly converges to a local solution rather than the global optimal solution, while the other three methods are much better in finding the global optimum. It is interesting to note that the results presented in Table 5 are better than that of Table 4 and Table 3 when a larger basis is allowed in tensor factorization. Moreover, in this case, the proximal BCD usually consumes less number of iterations than the two ADMM variants.

Acknowledgements

We would like to thank Professor Renato D. C. Monteiro and two anonymous referees for their insightful comments, which helped improve this paper significantly.

References