A globally convergent algorithm for nonconvex optimization based on block coordinate update

Introduction

In this paper, we consider (nonconvex) optimization problems in the form of

Due to the lack of convexity, standard analysis tools such as convex inequalities and Fejér-monotonicity cannot be applied to establish the convergence of the iterate sequence. The case becomes more difficult when the problem is nonsmooth. In these cases, convergence analysis of existing algorithms is typically limited to objective convergence (to a possibly non-minimal value) or the convergence of a certain subsequence of iterates to a critical point. (Some exceptions will be reviewed below.) Although whole-sequence convergence is almost always observed, it is rarely proved. This deficiency abates some widely used algorithms. For example, KSVD only has nonincreasing monotonicity of its objective sequence, and iterative reweighted algorithms for sparse and low-rank recovery in only has subsequence convergence. Some other methods establish whole sequence convergence by assuming stronger conditions such as local convexity (on at least a part of the objective) and either unique or isolated limit points, which may be difficult to satisfy or to verify. In this paper, we aim to establish whole sequence convergence with conditions that are provably satisfied by a wide class of functions.

Block coordinate descent (BCD) (more precisely, block coordinate update) is very general and widely used for solving both convex and nonconvex problems in the form of (1) with multiple blocks of variables. Since only one block is updated at a time, it has a low per-iteration cost and small memory footprint. Recent literature has found BCD as a viable approach for “big data” problems.

In order to solve (1), we propose a block prox-linear (BPL) method, which updates a block of variables at each iteration by minimizing a prox-linear surrogate function. Specifically, at iteration $k$, a block $b_{k}\in\{1,\ldots,s\}$ is selected and ${\bf x}^{k}=({\bf x}_{1}^{k},\cdots,{\bf x}_{s}^{k})$ is updated as follows:

where $\alpha_{k}>0$ is a stepsize and $\hat{{\mathbf{x}}}_{i}^{k}$ is the extrapolation

While we can simply set $\omega_{k}=0$, appropriate $\omega_{k}>0$ can speed up the convergence; we will demonstrate this in the numerical results below. We can set the stepsize $\alpha_{k}=\frac{1}{\gamma L_{k}}$ with any $\gamma>1$, where $L_{k}>0$ is the Lipschitz constant of $\nabla_{{\mathbf{x}}_{i}}f({\mathbf{x}}_{\neq i}^{k-1},{\mathbf{x}}_{i})$ about ${\mathbf{x}}_{i}$. When $L_{k}$ is unknown or difficult to bound, we can apply backtracking on $\alpha_{k}$ under the criterion:

Special cases

When there is only one block, i.e., $s=1$, Algorithm 1 reduces to the well-known (accelerated) proximal gradient method (e.g., ). When the update block cycles from 1 through $s$, Algorithm 1 reduces to the cyclic block proximal gradient (Cyc-BPG) method in . We can also randomly shuffle the $s$ blocks at the beginning of each cycle. We demonstrate in section 3 that random shuffling leads to better numerical performance. When the update block is randomly selected following the probability $p_{i}>0$, where $\sum_{i=1}^{s}p_{i}=1$, Algorithm 1 reduces to the randomized block coordinate descent method (RBCD) (e.g., ). Unlike these existing results, we do not assume convexity.

In our analysis, we impose an essentially cyclic assumption — each block is selected for update at least once within every $T\geq s$ consecutive iterations — otherwise the order is arbitrary. Our convergence results apply to all the above special cases except RBCD, whose convergence analysis requires different strategies; see for the convex case and for the nonconvex case.

2 Kurdyka-Łojasiewicz property

To establish whole sequence convergence of Algorithm 1, a key assumption is the Kurdyka-Łojasiewicz (KL) property of the objective function $F$.

The KL property was introduced by Łojasiewicz for real analytic functions. Kurdyka extended it to functions of the $o$-minimal structure. Recently, the KL inequality (4) was further extended to nonsmooth sub-analytic functions . The work characterizes the geometric meaning of the KL inequality.

3 Related literature

There are many methods that solve general nonconvex problems. Methods in the papers , the books , and in the references therein, do not break variables into blocks. They usually have the properties of local convergence or subsequence convergence to a critical point, or global convergence in the terms of the violation of optimality conditions. Next, we review BCD methods.

BCD has been extensively used in many applications. Its original form, block coordinate minimization (BCM), which updates a block by minimizing the original objective with respect to that block, dates back to the 1950’s and is closely related to the Gauss-Seidel and SOR methods for linear equation systems. Its convergence was studied under a variety of settings (cf. and the references therein). The convergence rate of BCM was established under the strong convexity assumption for the multi-block case and under the general convexity assumption for the two-block case. To have even cheaper updates, one can update a block approximately, for example, by minimizing an approximate objective like was done in (2), instead of sticking to the original objective. The work is a block coordinate gradient descent (BCGD) method where taking a block gradient step is equivalent to minimizing a certain prox-linear approximation of the objective. Its whole sequence convergence and local convergence rate were established under the assumptions of a so-called local Lipschitzian error bound and the convexity of the objective’s nondifferentiable part. The randomized block coordinate descent (RBCD) method in randomly chooses the block to update at each iteration and is not essentially cyclic. Objective convergence was established , and the violation of the first-order optimization condition was shown to converge to zero . There is no iterate convergence result for RBCD.

Some special cases of Algorithm 1 have been analyzed in the literature. The work uses cyclic updates of a fixed order and assumes block-wise convexity; studies two blocks without extrapolation, namely, $s=2$ and $\hat{{\mathbf{x}}}_{i}^{k}={\mathbf{x}}_{i}^{k-1},\,\forall k$ in (2). A more general result is [5, Lemma 2.6], where three conditions for whole sequence convergence are given and are met by methods including averaged projection, proximal point, and forward-backward splitting. Algorithm 1, however, does not satisfy the three conditions in .

The extrapolation technique in (3) has been applied to accelerate the (block) prox-linear method for solving convex optimization problems (e.g., ). Recently, show that the (block) prox-linear iteration with extrapolation can still converge if the nonsmooth part of the problem is convex, while the smooth part can be nonconvex. Because of the convexity assumption, their convergence results do not apply to Algorithm 1 for solving the general nonconvex problem (1).

4 Contributions

We summarize the main contributions of this paper as follows.

We propose a block prox-linear (BPL) method for nonconvex smooth and nonsmooth optimization. Extrapolation is used to accelerate it. To our best knowledge, this is the first work of prox-linear acceleration for fully nonconvex problems (where both smooth and nonsmooth terms are nonconvex) with a convergence guarantee. However, we have not proved any improved convergence rate.

Assuming essentially cyclic updates of the blocks, we obtain the whole sequence convergence of BPL to a critical point with rate estimates, by first establishing subsequence convergence and then applying the Kurdyka-Łojasiewicz (KL) property. Furthermore, we tailor our convergence analysis to several existing algorithms, including non-convex regularized linear regression and nonnegative matrix factorization, to improve their existing convergence results.

We numerically tested BPL on nonnegative matrix and tensor factorization problems. At each cycle of updates, the blocks were randomly shuffled. We observed that BPL was very efficient and that random shuffling avoided local solutions more effectively than the deterministic cyclic order.

5 Notation and preliminaries

Since the update may be aperiodic, extra notation is used for when and how many times a block is updated. Let ${\mathcal{K}}[i,k]$ denote the set of iterations in which the $i$-th block has been selected to update till the $k$th iteration:

which is the number of times the $i$-th block has been updated till iteration $k$. For $k=1,\ldots,$ we have $\cup_{i=1}^{s}{\mathcal{K}}[i,k]=[k]\triangleq\{1,2,\ldots,k\}$ and $\sum_{i=1}^{s}d_{i}^{k}=k$.

The extrapolated point in (2) (for $i=b_{k}$) is computed from the last two updates of the same block:

for some weight $0\leq\omega_{k}\leq 1$. We partition the set of Lipschitz constants and the extrapolation weights into $s$ disjoint subsets as

Hence, for each block $i$, we have three sequences:

For simplicity, we take stepsizes and extrapolation weights as follows

We make the following definitions, which can be found in .

The limiting Fréchet subdifferential is denoted by ${\partial}F({\mathbf{x}})$ and defined as

where ${\mathcal{X}}_{1}\times{\mathcal{X}}_{2}$ denotes the Cartesian product of ${\mathcal{X}}_{1}$ and ${\mathcal{X}}_{2}$ .

A point ${\mathbf{x}}^{*}$ is called a critical point of $F$ if $\mathbf{0}\in\partial F({\mathbf{x}}^{*})$.

For a proper, lower semicontinuous function $r$, its proximal mapping ${\mathbf{prox}}_{r}(\cdot)$ is defined as

As $r$ is nonconvex, ${\mathbf{prox}}_{r}(\cdot)$ is generally set-valued. Using this notation, the update in (2) can be written as (assume $i=b_{k}$)

6 Organization

The rest of the paper is organized as follows. Section 2 establishes convergence results. Examples and applications are given in section 3, and finally section 4 concludes this paper.

Convergence analysis

In this section, we analyze the convergence of Algorithm 1. Throughout our analysis, we make the following assumptions.

Let $i=b_{k}$. $\nabla_{{\mathbf{x}}_{i}}f({\bf x}_{\neq i}^{k-1},{\bf x}_{i})$ has Lipschitz continuity constant $L_{k}$ with respect to ${\bf x}_{i}$, i.e.,

In Algorithm 1, within any $T$ consecutive iterations, every block is updated at least one time.

Our analysis proceeds with several steps. We first estimate the objective decrease after every iteration (see Lemma 5) and then establish a square summable result of the iterate differences (see Proposition 6). Through the square summable result, we show a subsequence convergence result that every limit point of the iterates is a critical point (see Theorem 7). Assuming the KL property (see Definition 1) on the objective function and the following monotonicity condition, we establish whole sequence convergence of our algorithm and also give estimate of convergence rate (see Theorems 11 and 13).

The weight $\omega_{k}$ is chosen so that $F({\mathbf{x}}^{k})\leq F({\mathbf{x}}^{k-1}),\,\forall k$.

We will show that a range of nontrivial $\omega_{k}>0$ always exists to satisfy Condition 2.1 under a mild assumption, and thus one can backtrack $\omega_{k}$ to ensure $F({\mathbf{x}}^{k})\leq F({\mathbf{x}}^{k-1}),\,\forall k$. Maintaining the monotonicity of $F({\mathbf{x}}^{k})$ can significantly improve the numerical performance of the algorithm, as shown in our numerical results below and also in . Note that subsequence convergence does not require this condition.

We begin our analysis with the following lemma. The proofs of all the lemmas and propositions are given in Appendix A.

Take $\alpha_{k}$ and $\omega_{k}$ as in (9). After each iteration $k$, it holds

where $c_{1}=\frac{1}{4},c_{2}=9$, $i=b_{k}$ and $j=d_{i}^{k}$.

Note that $d_{i}^{k}=d_{i}^{k-1}+1$ for $i=b_{k}$ and $d_{i}^{k}=d_{i}^{k-1},\forall i\neq b_{k}$. Adopting the convention that $\sum_{j=p}^{q}a_{j}=0$ when $q<p$, we can write (13) into

which will be used in our subsequent convergence analysis.

Using Lemma 5, we can have the following result, through which we show subsequence convergence of Algorithm 1.

Let $\{{\bf x}^{k}\}_{k\geq 1}$ be generated from Algorithm 1 with $\alpha_{k}$ and $\omega_{k}$ taken from (9). We have

Under Assumptions 1 through 3, let $\{{\bf x}^{k}\}_{k\geq 1}$ be generated from Algorithm 1 with $\alpha_{k}$ and $\omega_{k}$ taken from (9). Then any limit point $\bar{{\mathbf{x}}}$ of $\{{\bf x}^{k}\}_{k\geq 1}$ is a critical point of (1). If the subsequence $\{{\mathbf{x}}^{k}\}_{k\in\bar{{\mathcal{K}}}}$ converges to $\bar{{\mathbf{x}}}$, then

The existence of finite limit point is guaranteed if $\{{\mathbf{x}}^{k}\}_{k\geq 1}$ is bounded, and for some applications, the boundedness of $\{{\mathbf{x}}^{k}\}_{k\geq 1}$ can be satisfied by setting appropriate parameters in Algorithm 1; see examples in section 3. If $r_{i}$’s are continuous, (16) immediately holds. Since we only assume lower semi-continuity of $r_{i}$’s, $F({\mathbf{x}})$ may not converge to $F(\bar{{\mathbf{x}}})$ as ${\mathbf{x}}\to\bar{{\mathbf{x}}}$, so (16) is not obvious.

Assume $\bar{{\bf x}}$ is a limit point of $\{{\bf x}^{k}\}_{k\geq 1}$. Then there exists an index set ${\mathcal{K}}$ so that the subsequence $\{{\bf x}^{k}\}_{k\in{\mathcal{K}}}$ converging to $\bar{{\bf x}}$. From (15), we have $\|{\mathbf{x}}^{k-1}-{\mathbf{x}}^{k}\|\to 0$ and thus $\{{\bf x}^{k+\kappa}\}_{k\in{\mathcal{K}}}\to\bar{{\bf x}}$ for any $\kappa\geq 0$. Define

Take an arbitrary $i\in\{1,\ldots,s\}$. Note ${\mathcal{K}}_{i}$ is an infinite set according to Assumption 3. Taking another subsequence if necessary, $L_{k}$ converges to some $\bar{L}_{i}$ as ${\mathcal{K}}_{i}\ni k\to\infty$. Note that since $\alpha_{k}=\frac{1}{2L_{k}},\forall k$, for any $k\in{\mathcal{K}}_{i}$,

Note from (15) and (6) that $\hat{{\mathbf{x}}}_{i}^{k}\to\bar{{\mathbf{x}}}_{i}$ as ${\mathcal{K}}_{i}\ni k\to\infty$. Since $f$ is continuously differentiable and $r_{i}$ is lower semicontinuous, letting ${\mathcal{K}}_{i}\ni k\to\infty$ in (17) yields

and $\bar{{\bf x}}_{i}$ satisfies the first-order optimality condition:

Since (18) holds for arbitrary $i\in\{1,\ldots,s\}$, $\bar{{\bf x}}$ is a critical point of (1).

Taking limit superior on both sides of the above inequality over $k\in{\mathcal{K}}_{i}$ gives $\underset{{\mathcal{K}}_{i}\ni k\to\infty}{\limsup}r_{i}({\bf x}_{i}^{k})\leq r_{i}(\bar{{\bf x}}_{i}).$ Since $r_{i}$ is lower semi-continuous, it holds $\underset{{\mathcal{K}}_{i}\ni k\to\infty}{\liminf}r_{i}({\bf x}_{i}^{k})\geq r_{i}(\bar{{\bf x}}_{i})$, and thus

Noting that $f$ is continuous, we complete the proof. ∎

2 Whole sequence convergence and rate

In this subsection, we establish the whole sequence convergence and rate of Algorithm 1 by assuming Condition 2.1. We first show that under mild assumptions, Condition 2.1 holds for certain $\omega_{k}>0$.

Let $i=b_{k}$. Assume ${\mathbf{prox}}_{\alpha_{k}r_{i}}$ is single-valued near ${\mathbf{x}}_{i}^{k-1}-\alpha_{k}\nabla_{{\mathbf{x}}_{i}}f({\mathbf{x}}^{k-1})$ and

namely, progress can still be made by updating the $i$-th block. Then, there is $\bar{\omega}_{k}>0$ such that for any $\omega_{k}\in[0,\bar{\omega}_{k}]$, we have $F({\mathbf{x}}^{k})\leq F({\mathbf{x}}^{k-1})$.

By this proposition, we can find $\omega_{k}>0$ through backtracking to maintain the monotonicity of $F({\mathbf{x}}^{k})$. All the examples in section 3 satisfy the assumptions of Proposition 8. The proof of Proposition 8 involves the continuity of ${\mathbf{prox}}_{\alpha_{k}r_{i}}$ and is deferred to Appendix A.4.

Under Condition 2.1 and the KL property of $F$ (Definition 1), we show that the sequence $\{{\bf x}^{k}\}$ converges as long as it has a finite limit point. We first establish a lemma, which has its own importance and together with the KL property implies Lemma 2.6 of .

The result in Lemma 9 below is very general because we need to apply it to Algorithm 1 in its general form. To ease understanding, let us go over its especial cases. If $s=1$, $n_{1,m}=m$ and $\beta=0$, then (21) below with $\alpha_{1,m}=\alpha_{m}$ and $A_{1,m}=A_{m}$ reduces to $\alpha_{m+1}A_{m+1}^{2}\leq B_{m}A_{m}$, which together with Young’s inequality gives $\sqrt{\underline{\alpha}}A_{m+1}\leq\frac{\sqrt{\underline{\alpha}}}{2}A_{m}+\frac{1}{2\sqrt{\underline{\alpha}}}B_{m}$. Hence, if $\{B_{m}\}_{m\geq 1}$ is summable, so will be $\{A_{m}\}_{m\geq 1}$. This result can be used to analyze the prox-linear method. The more general case of $s>1$, $n_{i,m}=m,\,\forall i$ and $\beta=0$ applies to the cyclic block prox-linear method. In this case, (21) reduces to $\sum_{i=1}^{s}\alpha_{i,m+1}A_{i,m+1}^{2}\leq B_{m}\sum_{i=1}^{s}A_{i,m},$ which together with the Young’s inequality implies

where $\tau$ is sufficiently large so that $\frac{1}{\tau}<\sqrt{\underline{\alpha}}$. Less obviously but still, if $\{B_{m}\}_{m\geq 1}$ is summable, so will be $\{A_{i,m}\}_{m\geq 1},\,\forall i$. Finally, we will need $\beta>0$ in (21) to analyze the accelerated block prox-linear method.

For nonnegative sequences $\{A_{i,j}\}_{j\geq 0},\{\alpha_{i,j}\}_{j\geq 0},\,i=1,\ldots,s$, and $\{B_{m}\}_{m\geq 0}$, if

where $0\leq\beta<1$, and $\{n_{i,m}\}_{m\geq 0},\forall i$ are nonnegative integer sequences satisfying: $n_{i,m}\leq n_{i,m+1}\leq n_{i,m}+N,\forall i,m$, for some integer $N>0$. Then we have

In addition, if $\sum_{m=1}^{\infty}B_{m}<\infty$, $\lim_{m\to\infty}n_{i,m}=\infty,\forall i$, and (21) holds for all $m$, then we have

The proof of this lemma is given in Appendix A.5

Let $\{{\bf x}^{k}\}$ be generated from Algorithm 1. For a specific iteration $k\geq 3T$, assume ${\mathbf{x}}^{\kappa}\in{\mathcal{B}}_{\rho}(\bar{{\mathbf{x}}}),\,\kappa=k-3T,k-3T+1,\ldots,k$ for some $\bar{{\mathbf{x}}}$ and $\rho>0$. If for each $i$, $\nabla_{{\mathbf{x}}_{i}}f({\mathbf{x}})$ is Lipschitz continuous with constant $L_{G}$ within $B_{4\rho}(\bar{{\mathbf{x}}})$ with respect to ${\mathbf{x}}$, i.e.,

We are now ready to present and show the whole sequence convergence of Algorithm 1.

Suppose that Assumptions 1 through 3 and Condition 2.1 hold. Let $\{{\bf x}^{k}\}_{k\geq 1}$ be generated from Algorithm 1. Assume

$\{{\bf x}^{k}\}_{k\geq 1}$ has a finite limit point $\bar{{\bf x}}$;

$F$ satisfies the KL property (4) around $\bar{{\bf x}}$ with parameters $\rho$, $\eta$ and $\theta$.

For each $i$, $\nabla_{{\mathbf{x}}_{i}}f({\bf x})$ is Lipschitz continuous within $B_{4\rho}(\bar{{\mathbf{x}}})$ with respect to ${\bf x}$.

Before proving the theorem, let us remark on the conditions 1–3. The condition 1 can be guaranteed if $\{{\mathbf{x}}^{k}\}_{k\geq 1}$ has a bounded subsequence. The condition 2 is satisfied for a broad class of applications as we mentioned in section 1.2. The condition 3 is a weak assumption since it requires the Lipschitz continuity only in a bounded set.

From (16) and Condition (2.1), we have $F({\mathbf{x}}^{k})\to F(\bar{{\mathbf{x}}})$ as $k\to\infty$. We consider two cases depending on whether there is an integer $K_{0}$ such that $F({\mathbf{x}}^{K_{0}})=F(\bar{{\mathbf{x}}})$.

Case 1: Assume $F({\mathbf{x}}^{k})>F(\bar{{\mathbf{x}}}),\,\forall k$.

Since $\bar{{\bf x}}$ is a limit point of $\{{\mathbf{x}}^{k}\}$ and according to (15), one can choose a sufficiently large $k_{0}$ such that the points ${\bf x}^{k_{0}+\kappa},\kappa=0,1,\ldots,3T$ are all sufficiently close to $\bar{{\bf x}}$ and in ${\mathcal{B}}_{\rho}(\bar{{\bf x}})$, and also the differences $\|{\bf x}^{k_{0}+\kappa}-{\bf x}^{k_{0}+\kappa+1}\|,\kappa=0,1,\ldots,3T$ are sufficiently close to zero. In addition, note that $F({\mathbf{x}}^{k})\to F(\bar{{\mathbf{x}}})$ as $k\to\infty$, and thus both $F({\mathbf{x}}^{3(k_{0}+1)T})-F(\bar{{\mathbf{x}}})$ and $\phi(F({\mathbf{x}}^{3(k_{0}+1)T})-F(\bar{{\mathbf{x}}}))$ can be sufficiently small. Since $\{{\mathbf{x}}^{k}\}_{k\geq 0}$ converges if and only if $\{{\mathbf{x}}^{k}\}_{k\geq k_{0}}$ converges, without loss of generality, we assume $k_{0}=0$, which is equivalent to setting ${\mathbf{x}}^{k_{0}}$ as a new starting point, and thus we assume

Assume that ${\mathbf{x}}^{3mT}\in{\mathcal{B}}_{\rho}(\bar{{\bf x}})$ and $F({\mathbf{x}}^{3mT})<F(\bar{{\mathbf{x}}})+\eta,\,m=0,\ldots,M$ for some $M\geq 1$. Note that from (25), we can take $M=1$. Letting $k=3mT$ in (24) and using KL inequality (4), we have

where $L_{G}$ is a uniform Lipschitz constant of $\nabla_{{\mathbf{x}}_{i}}f({\mathbf{x}}),\forall i$ within ${\mathcal{B}}_{4\rho}(\bar{{\mathbf{x}}})$. In addition, it follows from (14) that

Let $\phi_{m}=\phi(F({\bf x}^{3mT})-F(\bar{{\bf x}}))$. Note that

where $C$ is given in (26). Letting $N=1$ in the above inequality, we have

Case 2: Assume $F({\mathbf{x}}^{K_{0}})=F(\bar{{\mathbf{x}}})$ for a certain integer $K_{0}$.

Since $F({\mathbf{x}}^{k})$ is nonincreasingly convergent to $F(\bar{{\mathbf{x}}})$, we have $F({\mathbf{x}}^{k})=F(\bar{{\mathbf{x}}}),\,\forall k\geq K_{0}$. Take $M_{0}$ such that $3M_{0}T\geq K_{0}$. Then $F({\mathbf{x}}^{3mT})=F({\mathbf{x}}^{3(m+1)T})=F(\bar{{\mathbf{x}}}),\,\forall m\geq M_{0}$. Summing up (28) from $m=M\geq M_{0}$ gives

and thus ${\mathbf{x}}^{k}$ converges to the limit point $\bar{{\mathbf{x}}}$. This completes the proof. ∎

In addition, we can show convergence rate of Algorithm 1 through the following lemma.

For nonnegative sequence $\{A_{k}\}_{k=1}^{\infty}$, if $A_{k}\leq A_{k-1}\leq 1,\,\forall k\geq K$ for some integer $K$, and there are positive constants $\alpha,\beta$ and $\gamma$ such that

If $\gamma\geq 1$, then A_{k}\leq\big{(}\frac{\alpha+\beta}{1+\alpha+\beta}\big{)}^{k-K}A_{K},\,\forall k\geq K;

If $0<\gamma<1$, then $A_{k}\leq\nu(k-K)^{-\frac{\gamma}{1-\gamma}},\,\forall k\geq K,$ for some positive constant $\nu$.

Under the assumptions of Theorem 11, we have:

If $\theta\in[0,\frac{1}{2}]$, $\|{\bf x}^{k}-\bar{{\bf x}}\|\leq C\alpha^{k},\forall k$, for a certain $C>0,~{}\alpha\in[0,1)$;

If $\theta\in(\frac{1}{2},1)$, $\|{\bf x}^{k}-\bar{{\bf x}}\|\leq Ck^{-(1-\theta)/(2\theta-1)},\forall k$, for a certain $C>0$.

For $k>k_{0}$, since $F({\mathbf{x}}^{k-1})=F({\mathbf{x}}^{k})$, and noting that in (14) all terms but one are zero under the summation over $i$, we have

Note $\|{\mathbf{x}}^{m-1}-\bar{{\mathbf{x}}}\|\leq B_{m}$. Hence, choosing a sufficiently large $C>0$ gives the result in item 1 for $\theta=0$.

When $0<\theta<1$, if for some $k_{0}$, $F({\mathbf{x}}^{k_{0}})=F(\bar{{\mathbf{x}}})$, we have (36) by the same arguments as above and thus obtain linear convergence. Below we assume $F({\mathbf{x}}^{k})>F(\bar{{\mathbf{x}}}),\,\forall k$. Let

In addition, letting $N=m$ in (30), we have

where $C$ is the same as that in (30). Letting $M\to\infty$ in the above inequality, we have

where the second inequality is from (37). Since $A_{m-1}-A_{m}\leq 1$ as $m$ is sufficiently large and $\|{\mathbf{x}}^{m}-\bar{{\mathbf{x}}}\|\leq A_{\lfloor\frac{m}{3T}\rfloor}$, the results in item 2 for $\theta\in(0,\frac{1}{2}]$ and item 3 now immediately follow from Lemma 12. ∎

Before closing this section, let us make some comparison to the recent work . The whole sequence convergence and rate results in this paper are the same as those in . However, the results here cover more applications. We do not impose any convexity assumption on (1) while requires $f$ to be block-wise convex and every $r_{i}$ to be convex. In addition, the results in only apply to cyclic block prox-linear method. Empirically, a different block-update order can give better performance. As demonstrated in , random shuffling can often improve the efficiency of the coordinate descent method for linear support vector machine, and shows that for the Tucker tensor decomposition (see (47)), updating the core tensor more frequently can be better than cyclicly updating the core tensor and factor matrices.

Applications and numerical results

In this section, we give some specific examples of (1) and show the whole sequence convergence of some existing algorithms. In addition, we demonstrate that maintaining the nonincreasing monotonicity of the objective value can improve the convergence of accelerated gradient method and that updating variables in a random order can improve the performance of Algorithm 1 over that in the cyclic order.

FISTA is an accelerated proximal gradient method for solving composite convex problems. It is a special case of Algorithm 1 with $s=1$ and specific $\omega_{k}$’s. For the readers’ convenience, we present the method in Algorithm 2, where both $f$ and $g$ are convex functions, and $L_{f}$ is the Lipschitz constant of $\nabla f({\mathbf{x}})$. The algorithm reaches the optimal order of convergence rate among first-order methods, but in general, it does not guarantee monotonicity of the objective values. A restarting scheme is studied in that restarts FISTA from ${\mathbf{x}}^{k}$ whenever $F({\mathbf{x}}^{k+1})>F({\mathbf{x}}^{k})$ occursAnother restarting option is tested based on gradient information. It is demonstrated that the restarting FISTA can significantly outperform the original one. In this subsection, we show that FISTA with backtracking extrapolation weight can do even better than the restarting one.

We test the algorithms on solving the following problem

2 Coordinate descent method for nonconvex regression

As the number of predictors is larger than sample size, variable selection becomes important to keep more important predictors and obtain a more interpretable model, and penalized regression methods are popularly used to achieve variable selection. The work considers the linear regression with nonconvex penalties: the minimax concave penalty (MCP) and the smoothly clipped absolute deviation (SCAD) penalty . Specifically, the following model is considered

The cyclic coordinate descent method used in performs the update from $j=1$ through $p$

which can be equivalently written into the form of (2) by

Note that the data has been standardized such that $\|{\mathbf{x}}_{j}\|^{2}=n$. Hence, if $\gamma>1$ in (40) and $\gamma>2$ in (41), it is easy to verify that the objective in (42) is strongly convex, and there is a unique minimizer. From the convergence results of , it is concluded in that any limit pointIt is stated in that the sequence generated by (42) converges to a coordinate-wise minimizer of (38). However, the result is obtained directly from , which only guarantees subsequence convergence. of the sequence $\{\boldsymbol{\beta}^{k}\}$ generated by (42) is a coordinate-wise minimizer of (38). Since $r_{\lambda,\gamma}$ in both (40) and (41) is piecewise polynomial and thus semialgebraic, it satisfies the KL property (see Definition 1). In addition, let $f(\boldsymbol{\beta})$ be the objective of (38). Then

where $\mu$ is the strong convexity constant of the objective in (42). Hence, according to Theorem 11 and Remark 2.1, we have the following convergence result.

Assume ${\mathbf{X}}$ is standardized as in (39). Let $\{\boldsymbol{\beta}^{k}\}$ be the sequence generated from (42) or by the following update with random shuffling of coordinates

where $(\pi^{k}_{1},\ldots,\pi^{k}_{p})$ is any permutation of $(1,\ldots,p)$, and $r_{\lambda,\gamma}$ is given by either (40) with $\gamma>1$ or (41) with $\gamma>2$. If $\{\boldsymbol{\beta}^{k}\}$ has a finite limit point, then $\boldsymbol{\beta}^{k}$ converges to a coordinate-wise minimizer of (38).

3 Rank-one residue iteration for nonnegative matrix factorization

The nonnegative matrix factorization can be modeled as

In the literature, most existing algorithms for solving (43) update ${\mathbf{X}}$ and ${\mathbf{Y}}$ alternatingly; see the review paper and the references therein. The work partitions the variables in a different way: $({\mathbf{x}}_{1},{\mathbf{y}}_{1},\ldots,{\mathbf{x}}_{p},{\mathbf{y}}_{p})$, where ${\mathbf{x}}_{j}$ denotes the $j$-th column of ${\mathbf{X}}$, and proposes the rank-one residue iteration (RRI) method. It updates the variables cyclically, one column at a time. Specifically, RRI performs the updates cyclically from $i=1$ through $p$,

where ${\mathbf{X}}_{>i}^{k}=({\mathbf{x}}_{i+1}^{k},\ldots,{\mathbf{x}}_{p}^{k})$. It is a cyclic block minimization method, a special case of . The advantage of RRI is that each update in (44) has a closed form solution. Both updates in (44) can be written in the form of (2) by noting that they are equivalent to

Since $f({\mathbf{X}},{\mathbf{Y}})+r_{1}({\mathbf{X}})+r_{2}({\mathbf{Y}})$ is semialgebraic and has the KL property, directly from Theorem 11, we have the following whole sequence convergence, which is stronger compared to the subsequence convergence in .

Let $\{({\mathbf{X}}^{k},{\mathbf{Y}}^{k})\}_{k=1}^{\infty}$ be the sequence generated by (44) or (45) from any starting point $({\mathbf{X}}^{0},{\mathbf{Y}}^{0})$. If $\{{\mathbf{x}}_{i}^{k}\}_{i,k}$ and $\{{\mathbf{y}}_{i}^{k}\}_{i,k}$ are uniformly bounded and away from zero, then $({\mathbf{X}}^{k},{\mathbf{Y}}^{k})$ converges to a critical point of (43).

However, during the iterations of RRI, it may happen that some columns of ${\mathbf{X}}$ and ${\mathbf{Y}}$ become or approach to zero vector, or some of them blow up, and these cases fail the assumption of Theorem 15. To tackle with the difficulties, we modify the updates in (44) and improve the RRI method as follows.

Our first modification is to require each column of ${\mathbf{X}}$ to have unit Euclidean norm; the second modification is to take the Lipschitz constant of $\nabla_{{\mathbf{x}}_{i}}f({\mathbf{X}}_{<i}^{k+1},{\mathbf{x}}_{i},{\mathbf{X}}_{>i}^{k},{\mathbf{Y}}_{<i}^{k+1},{\mathbf{Y}}_{\geq i}^{k})$ to be $L_{i}^{k}=\max(L_{\min},\|{\mathbf{y}}_{i}^{k}\|^{2})$ for some $L_{\min}>0$; the third modification is that at the beginning of the $k$-th cycle, we shuffle the blocks to a permutation $(\pi_{1}^{k},\ldots,\pi_{p}^{k})$. Specifically, we perform the following updates from $i=1$ through $p$,

Note that if $\pi_{i}^{k}=i$ and $L_{i}^{k}=\|{\mathbf{y}}_{i}^{k}\|^{2}$, the objective in (46a) is the same as that in (45a). Both updates in (46) have closed form solutions; see Appendix B. Using Theorem 11, we have the following theorem, whose proof is given in Appendix C.1. Compared to the original RRI method, the modified one automatically has bounded sequence and always has the whole sequence convergence.

Let $\{({\mathbf{X}}^{k},{\mathbf{Y}}^{k})\}_{k=1}^{\infty}$ be the sequence generated by (46) from any starting point $({\mathbf{X}}^{0},{\mathbf{Y}}^{0})$. Then $\{{\mathbf{Y}}^{k}\}$ is bounded, and $({\mathbf{X}}^{k},{\mathbf{Y}}^{k})$ converges to a critical point of (43).

Numerical tests. We tested (45) and (46) on randomly generated data and also the Swimmer dataset . We set $L_{\min}=0.001$ in the tests and found that (46) with $\pi_{i}^{k}=i,\forall i,k$ produced the same final objective values as those by (45) on both random data and the Swimmer dataset. In addition, (46) with random shuffling performed almost the same as those with $\pi_{i}^{k}=i,\forall i$ on randomly generated data. However, random shuffling significantly improved the performance of (46) on the Swimmer dataset. There are 256 images of resolution $32\times 32$ in the Swimmer dataset, and each image (vectorized to one column of ${\mathbf{M}}$) is composed of four limbs and the body. Each limb has four different positions, and all images have the body at the same position; see Figure 2. Hence, each of these images is a nonnegative combination of 17 images: one with the body and each one of another 16 images with one limb. We set $p=17$ in our test and ran (45) and (46) with/without random shuffling to 100 cycles. If the relative error ${\|{\mathbf{X}}^{out}({\mathbf{Y}}^{out})^{\top}-{\mathbf{M}}\|_{F}}/{\|{\mathbf{M}}\|_{F}}$ is below $10^{-3}$, we regard the factorization to be successful, where $({\mathbf{X}}^{out},{\mathbf{Y}}^{out})$ is the output. We ran the three different updates for 50 times independently, and for each run, they were fed with the same randomly generated starting point. Both (45) and (46) without random shuffling succeed 20 times, and (46) with random shuffling succeeds 41 times. Figure 3 plots all cases that occur. Every plot is in terms of running time (sec), and during that time, both methods run to 100 cycles. Since (45) and (46) without random shuffling give exactly the same results, we only show the results by (46). From the figure, we see that (46) with fixed cyclic order and with random shuffling has similar computational complexity while the latter one can more frequently avoid bad local solutions.

4 Block prox-linear method for nonnegative Tucker decomposition

The nonnegative Tucker decomposition is to decompose a given nonnegative tensor (multi-dimensional array) into the product of a core nonnegative tensor and a few nonnegative factor matrices. It can be modeled as

where ${\mathbf{A}}=({\mathbf{A}}_{1},\ldots,{\mathbf{A}}_{N})$ and $\boldsymbol{{\mathcal{X}}}\times_{i}{\mathbf{Y}}$ denotes tensor-matrix multiplication along the $i$-th mode (see for example). The cyclic block proximal gradient method for solving (47) performs the following updates cyclically

Here, $f(\boldsymbol{{\mathcal{C}}},{\mathbf{A}})=\frac{1}{2}\|\boldsymbol{{\mathcal{C}}}\times_{1}{\mathbf{A}}_{1}\ldots\times_{N}{\mathbf{A}}_{N}-\boldsymbol{{\mathcal{M}}}\|_{F}^{2}$, $L_{c}^{k}$ and $L_{i}^{k}$ (chosen no less than a positive $L_{\min}$) are gradient Lipschitz constants with respect to $\boldsymbol{{\mathcal{C}}}$ and ${\mathbf{A}}_{i}$ respectively, and $\hat{\boldsymbol{{\mathcal{C}}}}^{k}$ and $\hat{{\mathbf{A}}}_{i}^{k}$ are extrapolated points:

where $\omega_{k}$ is the same as that in Algorithm 2. Our setting of extrapolated points exactly follows . Figure 4 shows that the extrapolation technique significantly accelerates the convergence speed of the method. Note that the block-prox method with no extrapolation reduces to the block coordinate gradient method in .

Since the core tensor $\boldsymbol{{\mathcal{C}}}$ interacts with all factor matrices, the work proposes to update $\boldsymbol{{\mathcal{C}}}$ more frequently to improve the performance of the block proximal gradient method. Specifically, at each cycle, it performs the following updates sequentially from $i=1$ through $N$

It was demonstrated that (51) numerically performs better than (48). Numerically, we observed that the performance of (51) could be further improved if the blocks of variables were randomly shuffled as in (46), namely, we performed the updates sequentially from $i=1$ through $N$

where $(\pi^{k}_{1},\pi^{k}_{2},\ldots,\pi^{k}_{N})$ is a random permutation of $(1,2,\ldots,N)$ at the $k$-th cycle. Note that both (48) and (52) are special cases of Algorithm 1 with $T=N+1$ and $T=2N+2$ respectively. If $\{(\boldsymbol{{\mathcal{C}}}^{k},{\mathbf{A}}^{k})\}$ is bounded, then so are $L_{c}^{k},L_{c}^{k,i}$ and $L_{i}^{k}$’s. Hence, by Theorem 11, we have the convergence result as follows.

The sequence $\{(\boldsymbol{{\mathcal{C}}}^{k},{\mathbf{A}}^{k})\}$ generated from (48) or (52) is either unbounded or converges to a critical point of (47).

We tested (51) and (52) on the $32\times 32\times 256$ Swimmer dataset used above and set the core size to $24\times 17\times 16$. We ran them to 500 cycles from the same random starting point. If the relative error $\|\boldsymbol{{\mathcal{C}}}^{out}\times_{1}{\mathbf{A}}_{1}^{out}\ldots\times_{N}{\mathbf{A}}_{N}^{out}-\boldsymbol{{\mathcal{M}}}\|_{F}/\|\boldsymbol{{\mathcal{M}}}\|_{F}$ is below $10^{-3}$, we regard the decomposition to be successful, where $(\boldsymbol{{\mathcal{C}}}^{out},{\mathbf{A}}^{out})$ is the output. Among 50 independent runs, (52) with random shuffling succeeds 21 times while (51) succeeds only 11 times. Figure 5 plots all cases that occur. Similar to Figure 3, every plot is in terms of running time (sec), and during that time, both methods run to 500 iterations. From the figure, we see that (52) with fixed cyclic order and with random shuffling has similar computational complexity while the latter one can more frequently avoid bad local solutions.

Conclusions

Acknowledgements

The authors would like to thank three anonymous referees for their careful reviews and constructive comments.

Appendix A Proofs of key lemmas

In this section, we give proofs of the lemmas and also propositions we used.

Since ${\mathbf{x}}_{i}^{k}$ is the minimizer of (2), then

Summing (53) and (54) and noting that ${\mathbf{x}}_{j}^{k+1}={\mathbf{x}}_{j}^{k},\forall j\neq i$, we have

A.2 Proof of the claim in Remark 2.2

Assume $b_{k}=i$ and $\alpha_{k}=\frac{1}{L_{k}}$. When $f$ is block multi-convex and $r_{i}$ is convex, from Lemma 2.1 of , it follows that

A.3 Proof of Proposition 6

Summing (14) over $k$ from $1$ to $K$ gives

A.4 Proof of Proposition 8

From Corollary 5.20 and Example 5.23 of , we have that if ${\mathbf{prox}}_{\alpha_{k}r_{i}}$ is single valued near ${\mathbf{x}}_{i}^{k-1}-\alpha_{k}\nabla_{{\mathbf{x}}_{i}}f({\mathbf{x}}^{k-1})$, then ${\mathbf{prox}}_{\alpha_{k}r_{i}}$ is continuous at ${\mathbf{x}}_{i}^{k-1}-\alpha_{k}\nabla_{{\mathbf{x}}_{i}}f({\mathbf{x}}^{k-1})$. Let $\hat{{\mathbf{x}}}^{k}_{i}(\omega)$ explicitly denote the extrapolated point with weight $\omega$, namely, we take $\hat{{\mathbf{x}}}^{k}_{i}(\omega_{k})$ in (6). In addition, let {\mathbf{x}}^{k}_{i}(\omega)={\mathbf{prox}}_{\alpha_{k}r_{i}}\big{(}\hat{{\mathbf{x}}}_{i}^{k}(\omega)-\alpha_{k}\nabla_{{\mathbf{x}}_{i}}f({\bf x}_{\neq i}^{k-1},\hat{{\bf x}}_{i}^{k}(\omega))\big{)}. Note that (14) implies

From the optimality of ${\mathbf{x}}^{k}_{i}(\omega)$, it holds that

Taking limit superior on both sides of the above inequality, we have

which implies $\underset{\omega\to 0^{+}}{\limsup}\,r_{i}({\bf x}_{i}^{k}(\omega))\leq r_{i}({\bf x}_{i}^{k}(0))$. Since $r_{i}$ is lower semicontinuous, $\underset{\omega\to 0^{+}}{\liminf}\,r_{i}({\bf x}_{i}^{k}(\omega))\geq r_{i}({\bf x}_{i}^{k}(0))$. Hence, $\underset{\omega\to 0^{+}}{\lim}r_{i}({\bf x}_{i}^{k}(\omega))=r_{i}({\bf x}_{i}^{k}(0))$, and thus $\underset{\omega\to 0^{+}}{\lim}F({\bf x}^{k}(\omega))=F({\bf x}^{k}(0))$. Together with (63), we conclude that there exists $\bar{\omega}_{k}>0$ such that $F({\mathbf{x}}^{k-1})-F({\bf x}^{k}(\omega))\geq 0,\,\forall\omega\in[0,\bar{\omega}_{k}]$. This completes the proof.

A.5 Proof of Lemma 9

Let ${\mathbf{a}}_{m}$ and ${\mathbf{u}}_{m}$ be the vectors with their $i$-th entries

By the Cauchy-Schwarz inequality and noting $n_{i,m+1}-n_{i,m}\leq N,\forall i,m$, we have

Summing the above inequality over $m$ from $M_{1}$ through $M_{2}\leq M$ and arranging terms gives

which together with $\sum_{i=1}^{s}A_{i,n_{i,m+1}}\leq\sqrt{s}\|{\mathbf{u}}_{m+1}\|$ gives

where we have used $\|{\mathbf{u}}_{M_{1}}\|\leq\sum_{i=1}^{s}A_{i,n_{i,M_{1}}}$, and

where the inequality can be verified by noting $(1-\beta^{2})(4-(1+\beta)^{2})-(1-\beta)^{2}$ is decreasing with respect to $\beta$ in $$. Thus from (78) and (83), we have$C_{2}=\frac{1}{2C_{1}},\,C_{3}=\frac{C_{1}}{2},\,C_{4}=\beta\sqrt{\overline{\alpha}}+\frac{\sqrt{s}C_{1}}{2}$. Hence, from (81), we complete the proof of (22).

If $\lim_{m\to\infty}n_{i,m}=\infty,\forall i$, $\sum_{m=1}^{\infty}B_{m}<\infty$, and (21) holds for all $m$, letting $M_{1}=1$ and $M_{2}\to\infty$, we have (23) from (81).

A.6 Proof of Proposition 10

Note that ${\mathbf{x}}_{i}$ may be updated to ${\mathbf{x}}_{i}^{k}$ not at the $k$-th iteration but at some earlier one, which must be between $k-T$ and $k$ by Assumption 3. In addition, for each pair $(i,j)$, there must be some $\kappa_{i,j}$ between $k-2T$ and $k$ such that

By triangle inequality, $({\mathbf{y}}_{\neq i}^{(i)},{\mathbf{z}}_{i})\in B_{4\rho}(\bar{{\mathbf{x}}})$ for all $i$. Therefore, it follows from (10) and (84) that

A.7 Proof of Lemma 12

The proof follows that of Theorem 2 of . When $\gamma\geq 1$, since $0\leq A_{k-1}-A_{k}\leq 1,\forall k\geq K$, we have $(A_{k-1}-A_{k})^{\gamma}\leq A_{k-1}-A_{k}$, and thus (34) implies that for all $k\geq K$, it holds that $A_{k}\leq(\alpha+\beta)(A_{k-1}-A_{k})$, from which item 1 immediately follows.

When $\gamma<1$, we have $(A_{k-1}-A_{k})^{\gamma}\geq A_{k-1}-A_{k}$, and thus (34) implies that for all $k\geq K$, it holds that $A_{k}\leq(\alpha+\beta)(A_{k-1}-A_{k})^{\gamma}$. Letting $h(x)=x^{-1/\gamma}$, we have for $k\geq K$,

where we have used nonincreasing monotonicity of $h$ in the second inequality. Hence,

Let $\mu$ be the positive constant such that

Note that the above equation has a unique solution $0<\mu<1$. We claim that

It obviously holds from (92) and (93) if \big{(}\frac{A_{k}}{A_{k-1}}\big{)}^{1/\gamma}\geq\mu. It also holds if \big{(}\frac{A_{k}}{A_{k-1}}\big{)}^{1/\gamma}\leq\mu from the arguments

where the last inequality is from $A_{k-1}^{1-1/\gamma}\geq 1$. Hence, (94) holds, and summing it over $k$ gives

which immediately gives item 2 by letting $\nu=(\mu^{\gamma-1}-1)^{\frac{\gamma}{\gamma-1}}$.

Appendix B Solutions of (46)

In this section, we give closed form solutions to both updates in (46). First, it is not difficult to have the solution of (46b):

Secondly, since $L_{\pi_{i}}^{k}>0$, it is easy to write (46a) in the form of

which ${\mathbf{c}}={\mathbf{a}}-{\mathbf{b}}$. Next we give solution to (95) in three different cases.

Case 1: ${\mathbf{c}}<0$. Let $i_{0}=\operatorname*{arg\,max}_{i}c_{i}$ and $c_{\max}=c_{i_{0}}<0$. If there are more than one components equal $c_{\max}$, one can choose an arbitrary one of them. Then the solution to (95) is given by $x_{i_{0}}=1$ and $x_{i}=0,\forall i\neq i_{0}$ because for any ${\mathbf{x}}\geq 0$ and $\|{\mathbf{x}}\|=1$, it holds that

Case 2: ${\mathbf{c}}\leq 0$ and ${\mathbf{c}}\not<0$. Let ${\mathbf{c}}=({\mathbf{c}}_{I_{0}},{\mathbf{c}}_{I_{-}})$ where ${\mathbf{c}}_{I_{0}}=\mathbf{0}$ and ${\mathbf{c}}_{I_{-}}<0$. Then the solution to (95) is given by ${\mathbf{x}}_{I_{-}}=\mathbf{0}$ and ${\mathbf{x}}_{I_{0}}$ being any vector that satisfies ${\mathbf{x}}_{I_{0}}\geq 0$ and $\|{\mathbf{x}}_{I_{0}}\|=1$ because ${\mathbf{c}}^{\top}{\mathbf{x}}\leq 0$ for any ${\mathbf{x}}\geq 0$.

Case 3: ${\mathbf{c}}\not\leq 0$. Let ${\mathbf{c}}=({\mathbf{c}}_{I_{+}},{\mathbf{c}}_{I_{+}^{c}})$ where ${\mathbf{c}}_{I_{+}}>0$ and ${\mathbf{c}}_{I_{+}^{c}}\leq 0$. Then (95) has a unique solution given by ${\mathbf{x}}_{I_{+}}=\frac{{\mathbf{c}}_{I_{+}}}{\|{\mathbf{c}}_{I_{+}}\|}$ and ${\mathbf{x}}_{I_{+}^{c}}=\mathbf{0}$ because for any ${\mathbf{x}}\geq 0$ and $\|{\mathbf{x}}\|=1$, it holds that

where the second inequality holds with equality if and only if ${\mathbf{x}}_{I_{+}}$ is collinear with ${\mathbf{c}}_{I_{+}}$, and the third inequality holds with equality if and only if ${\mathbf{x}}_{I_{+}^{c}}=\mathbf{0}$.

Appendix C Proofs of convergence of some examples

In this section, we give the proofs of the theorems in section 3.

Through checking the assumptions of Theorem 11, we only need to verify the boundedness of $\{{\mathbf{Y}}^{k}\}$ to show Theorem 16. Let ${\mathbf{E}}^{k}={\mathbf{X}}^{k}({\mathbf{Y}}^{k})^{\top}-{\mathbf{M}}$. Since every iteration decreases the objective, it is easy to see that $\{{\mathbf{E}}^{k}\}$ is bounded. Hence, $\{{\mathbf{E}}^{k}+{\mathbf{M}}\}$ is bounded, and

Let $y_{ij}^{k}$ be the $(i,j)$-th entry of ${\mathbf{Y}}^{k}$. Thus the columns of ${\mathbf{E}}^{k}+{\mathbf{M}}$ satisfy

where ${\mathbf{x}}_{j}^{k}$ is the $j$-th column of ${\mathbf{X}}^{k}$. Since $\|{\mathbf{x}}_{j}^{k}\|=1$, we have $\|{\mathbf{x}}_{j}^{k}\|_{\infty}\geq 1/\sqrt{m},\,\forall j$. Note that (96) implies each component of $\sum_{j=1}^{p}y_{ij}^{k}{\mathbf{x}}_{j}^{k}$ is no greater than $a$. Hence from nonnegativity of ${\mathbf{X}}^{k}$ and ${\mathbf{Y}}^{k}$ and noting that at least one entry of ${\mathbf{x}}_{j}^{k}$ is no less than $1/\sqrt{m}$, we have $y_{ij}^{k}\leq a\sqrt{m}$ for all $i,j$ and $k$. This completes the proof.

References