Accelerated Stochastic Algorithms for Nonconvex Finite-sum and Multi-block Optimization

Introduction

Nonconvex optimization plays a fundamental role in modern statistics and machine learning, e.g., for empirical risk minimization with either nonconvex loss () or regularization (), as well as the training of deep neural networks (). In this paper, we consider two classes of nonconvex optimization problems that are widely used in statistical learning. The first class of problems intends to minimize the summation of many terms:

Moreover, we assume that there exists $0<\mu\leq L$ such that (s.t.)

Clearly, (1.2) implies (1.3) (with $\mu=L$). While in the classical nonlinear programming setting one only assumes (1.2), by using both conditions (1.2) and (1.3) we can explore more structural information for the design of solution methods of problem (1.1). In particular, we intend to develop more efficient algorithms to solve problems where the condition number $L/\mu$ associated with problem (1.1) is large. As an example, consider the nonconvex composite problem arising from variable selection in statistics : $f(x)=\textstyle{\tfrac{1}{m}\sum_{i=1}^{m}}h_{i}(x)+\rho p(x),$ where $h_{i}$’s are smooth convex functions, $p$ is a nonconvex function, and $\rho>0$ is a relatively small penalty parameter. Note that some examples of the nonconvex penalties are given by minimax concave penalty (MCP) or smoothly clipped absolute deviation (SCAD) (see ). It can be shown that the condition number for these problems is usually larger than $m$ (see Section 4 for more details).

In addition to (1.1), we consider an important class of nonconvex multi-block optimization problems with linearly coupled constraints, i.e.,

A different line of research aims to incorporate Nesterov’s acceleration (momentum) into nonconvex optimization. Ghaidmi and Lan first established the convergence of the accelerated gradient method for nonconvex optimization and show that it can improve the complexity of GD if the problem has a large condition number (i.e., $L/\mu$ is large). Their results were further improved in , , and . Currently the best complexity result, in terms of total gradient computations, for these methods is given by $\mathcal{O}\left(m\sqrt{L\mu}/\epsilon\right)$ for unconstrained problems . However, it remains unknown if the complexity of such accelerated algorithms can be further improved in terms of the dependence on $m$, especially when one needs to maintain the ${\cal O}(1/\epsilon)$ complexity bound on gradient computations. Note that some nonconvex stochastic accelerated gradient methods have been discussed in but they all exhibit a worse ${\cal O}(L\sigma^{2}/\epsilon^{2})$ complexity bounds.

While stochastic and randomized methods are being intensively explored for solving problem (1.1), most existing studies for the nonconvex multi-block problem in (1) have been mainly focused on deterministic batch methods only. Many of these studies aim at the generalization of the alternating direction method of multipliers (ADMM) method for nonconvex optimization. For example, In , Hong et al. established the complexity for a variant of ADMM for nonconvex multi-block problems, see also and for some previous work on the asymptotic analysis of ADMM for nonconvex optimization. In , Melo and Monteiro presented a linearized proximal multiblock ADMM with complexity $\mathcal{O}\left(1/\epsilon\right)$ to attain a nearly feasible $\epsilon$-stationary solution, but all the blocks have to be updated in each iteration. Later, they proposed a Jacobi-type ADMM in with similar complexity bound, which shows benefits if parallel computing is available. While the idea of randomly selecting blocks in nonconvex ADMM has been explored recently, these studies focus on the asymptotical convergence of these schemes (e.g., in ). To the best of our knowledge, there does not exist any complexity analysis regarding randomized methods for solving the nonconvex multi-block problem in (1) in the literature and as a consequence, it remains unclear whether stochastic or randomized methods are more advantageous over batch ones or not.

Our contribution in this paper mainly exists in the following several aspects. Firstly, we develop a new randomized algorithm, namely the randomized accelerated proximal gradient (RapGrad) method for solving problem (1.1) and show that it can significantly improve the complexity of existing algorithms especially for problems with a large condition number. More specifically, we show that RapGrad requires totally $\mathcal{O}(\mu(m+\sqrt{mL/\mu})/\epsilon)$ gradient computations in order to find a stochastic $\epsilon$-stationary point. For problems with $L/\mu\geq m$, this bound reduces to $\mathcal{O}(\sqrt{mL\mu}/\epsilon)$, which dominates the best-known batch accelerated gradient methods by a factor of $\sqrt{m}$ , and outperforms those variance-reduced stochastic algorithms by a factor of $m^{\frac{1}{6}}L^{\frac{1}{2}}/\mu^{\frac{1}{2}}$ (at least $m^{\frac{2}{3}}$). In fact, our complexity bound will be better than the latter algorithms as long as $L/\mu\log(L/\mu)>m^{\frac{1}{3}}$. Therefore, we provide some affirmative answers regarding whether the complexity bounds of variance reduced algorithms and accelerated gradient methods for nonconvex optimization can be further improved, especially in terms of their dependence on $m$. To the best of our knowledge, all these complexity results seem to be new in the literature for nonconvex finite-sum optimization. It is worth noting that some improvement over variance-reduced stochastic algorithms under the region $m\geq L/\mu$ (i.e., $L/\mu$ is small) has been presented recently in . RapGrad is a proximal-point type method which iteratively transforms the original nonconvex problem into a series of convex subproblems. In RapGrad, we incorporate a modified optimal randomized incremental gradient method, namely the randomized primal-dual gradient (see ) to solve these convex subproblems, and as a consequence, each iteration of RapGrad requires gradient computation for only one randomly selected component function. In comparison with existing nonconvex proximal-point type methods, the design and analysis of RapGrad appear to be more complicated. In particular, RapGrad does not require the computation of full gradients throughout its entire procedure by properly initializing a few intertwined search points and gradients using information obtained from the previous subproblems. This comes with the price of requiring additional storage (memory) for maintaining $\mathcal{O}(m)$ variables (e.g., $\underline{x}^{t}$). Moreover, the analysis of RapGrad requires us to show the convergence for some auxiliary sequences where the gradients are computed, which has not been established for the original randomized primal-dual gradient method.

primal block updates, where $\bar{A}=\max_{i\in[m-1]}\|A_{i}\|$, $\|\mathbf{A}\|^{2}=\sum_{i=1}^{m-1}\|A_{i}\|^{2}$, and $\mathcal{D}^{0}:=\sum_{i=1}^{m}[f_{i}(\bar{x}_{i}^{0})-f_{i}(x_{i}^{*})]$. Moreover, we demonstrate that the total number primal block updates that RapGrad requires can be much smaller, up to a factor of ${\cal O}(\sqrt{m})$, than its batch counterpart. To the best of our knowledge, this is the first time that the complexity of randomized methods for solving this special class of nonconvex multi-block optimization has been established and their possible advantages over batch methods are quantified in the literature.

Thirdly, we perform some numerical experiments on both RapGrad and RapDual for solving nonconvex finite-sum and multi-block problems in (1.1) and (1) and demonstrate their potential advantages over some existing algorithms.

This paper is organized as follows. In Section 2, we present our algorithm RapGrad, and its convergence properties for solvinwg the nonconvex finite-sum problem in (1.1). RapDual for nonconvex finite-sum optimization with linear constraints (1) and its convergence analysis are included in Section 3. Section 4 is devoted to some numerical experiments of our algorithms for the above two types of problems. Finally some concluding remarks are made in Section 5.

Nonconvex finite-sum optimization

In this section, we develop a randomized accelerated proximal gradient (RapGrad) method for solving the nonconvex finite-sum optimization problem in (1.1) and demonstrate that it can significantly improve the existing rates of convergence for solving these problems, especially when their objective functions have large condition numbers. We will describe this algorithm and establish its convergence in Subsections 2.1 and 2.2, respectively.

Before establishing the convergence of the RapGrad method, we first need to define an approximate stationary point for problem (1.1). A point $x\in X$ is called an approximate stationary point if it sits within a small neighborhood of a point $\hat{x}\in X$ which approximately satisfies the first-order optimality condition.

A point $x\in X$ is called an $(\epsilon,\delta)$-solution of (1.1) if there exists some $\hat{x}\in X$ such that

A stochastic $(\epsilon,\delta)$-solution of (1.1) is one such that

Moreover, if $X$ is a compact set and $x\in X$ is an $(\epsilon,\delta)$-solution, we can bound strong gap as follows:

where $D_{X}:=\max_{x_{1},x_{2}\in X}\|x_{1}-x_{2}\|$. In comparison with the two well-known criterions in (2.18) and (2.19), the criterion given in Definition 1 seems to be applicable to a wider class of problems and is particularly suitable for proximal-point type methods (see for a related notion).

We are now ready to state the main convergence properties for RapGrad.

In view of Theorem 2, we can bound the total number of gradient evaluations required by RapGrad to yield a stochastic $(\epsilon,\delta)$-solution of (1.1). Indeed, observe that the full gradient is computed only once in the first outer loop, and that for each subproblem (1.1), we only need to compute $s$ gradients with

Hence, the total number of gradient evaluations performed by RapGrad can be bounded by

where $D^{0}:=f(\bar{x}^{0})-f(x^{*})$. As a comparison, the batch version of this algorithm, obtained by viewing $\textstyle{\tfrac{1}{m}\sum_{i=1}^{m}}f_{i}(x)$ as a single component, would update all the $\underline{x}_{i}^{t}$ and $y_{i}^{t}$ for $i=1,\ldots,m$, in (2.12) and (2.15) at each iteration, and hence would require

gradient evaluations to compute an $(\epsilon,\delta)$-solution of (1.1).

It is worth noting that both complexity bounds $N(\epsilon,\delta)$ and $\hat{N}(\epsilon,\delta)$ will go to $+\infty$ as $\mu$ tends to . This indicates that the proximal-point method should not be applied to the situation when $\mu$ is too small, no matter whether the randomized or batch version is used. However, if $\mu$ is indeed small, say $\mu<\epsilon$, we can see that the assumption in (1.3) also holds for $\mu=\epsilon$. Therefore, we can assume that $\mu\geq\epsilon$ when applying RapGrad. A closer examination to this issue reveals that the algorithm (i.e., RaGrad) that we used to solve the strongly convex subproblems is not optimal if the strongly convex modulus $\mu$ is too small. Hence, a different remedy would be to apply an optimal randomized incremental gradient method which can yield the best possible complexity bound even if $\mu$ is small. Using this latter approach, we possibly do not need to fix $\epsilon$ a priori as in the former one. However, this type of optimal randomized incremental gradient method has not been developed until recently (see ), about one year later after our paper was initially released.

Theorem 2 only shows the convergence of RapGrad in expectation. Similarly to the nonconvex SGD methods in , we can establish and then further improve the convergence of RapGrad with overwhelming probability by using a two-phase procedure, where one computes a short list of candidate solutions in the optimization phase by either taking a few independent runs of RapGrad or randomly selecting a few solutions from the trajectory of RapGrad, and then chooses the best solution, e.g., in terms of either (2.18) and (2.19), in the post-optimization phase.

2 Convergence analysis for RapGrad

In this section, we will first develop the convergence results for Algorithm 2 applied to the convex finite-sum subproblem (2.8), and then using them to establish the convergence of RapGrad. Observe that the component functions $\psi_{i}$ and $\varphi$ in (2.8) satisfy:

$\tfrac{\mu}{2}\|x-y\|^{2}\leq\psi_{i}(x)-\psi_{i}(y)-\langle\nabla\psi_{i}(y),x-y\rangle\leq\tfrac{\hat{L}}{2}\|x-y\|^{2},\ \ \forall x,y\in X,\quad i=1,\ldots,m,$

$\varphi(x)-\varphi(y)-\langle\nabla\varphi(y),x-y\rangle\geq\tfrac{\mu}{2}\|x-y\|^{2},\ \ \forall x,y\in X,$

We first state some simple relations about the iterations generated by Algorithm 2.

Lemma 4 below describes an important result about Algorithm 2, which improves Lemma 7 of by showing the convergence of ${\underline{x}}^{s}_{i}$. The proof of this result is more involved and will be deferred in Appendix A.

Let the iterates $x^{t}$ and $y^{t}$, for $t=1,\ldots,s$, be generated by Algorithm 2 and $x^{*}$ be an optimal solution of (2.8). If the parameters in Algorithm 2 satisfy for all $t=1,\ldots,s-1$,

With the help of Lemma 4, we now establish the main convergence properties of Algorithm 2.

Let $x^{*}$ be an optimal solution of (2.8), and suppose that the parameters $\{\alpha_{t}\}$, $\{\tau_{t}\}$, $\{\eta_{t}\}$ and $\{\gamma_{t}\}$ are set as in (2.20) and (2.21). If $\varphi(x)=\tfrac{\mu}{2}\|x-z\|^{2}$, for some $z\in X$, then, for any $s\geq 1$, we have

It is easy to check that (2.20) and (2.21) satisfy conditions (2.24), (2.25), (2.26) (2.27), (2.28), and (2.29). Then by Lemma 4, we have

Since $\varphi(x)=\tfrac{\mu}{2}\|x-z\|^{2}$, we have $V_{\varphi}(x^{*},x^{s})=\tfrac{\mu}{2}\|x^{*}-x^{s}\|^{2}$, and $V_{\varphi}(x^{0},x^{s})=\tfrac{\mu}{2}\|x^{*}-x^{0}\|^{2}$. Plugging into (2.30), we obtain the following two relations:

In view of the above two relations, we have

In view of Theorem 5, Algorithm 2 applied to subproblem (2.8) exhibits a fast linear rate of convergence. Actually, as shown below we do not need to solve the subproblem too accurately, and a constant number of iteration of Algorithm 2 for each subproblem is enough to guarantee the convergence of Algorithm 1.

By induction on (2.32) and noting $\bar{\underline{x}}^{0}_{i}=\bar{x}^{0}$, $i=1,\ldots,m$, we have

Using (2.2), (2.2) and the condition on $s$, we have

Now we are ready to prove Theorem 2 using all the previous results we have developed.

From the optimality condition of (2.17), we obtain

Using the above relation and Lemma 6, we have

Nonconvex multi-block optimization with linear constraints

In this section, we present a randomized accelerated proximal dual (RapDual) algorithm for solving the nonconvex multi-block optimization problem in (1) and show the potential advantages in terms of the total number of block updates.

As mentioned in Section 1, we assume the inverse of the last block of the constraint matrix is easily computable. Hence, denoting $\mathbf{A}_{i}=A_{m}^{-1}A_{i}$, $i=1,\ldots,m-1$ and $\mathbf{b}=A_{m}^{-1}$, we can reformulate problem (1) as

where $f({\textbf{x}}):=\sum_{i=1}^{m-1}f_{i}(x_{i})$, $X=X_{1}\times\ldots\times X_{m-1}$, $\mathbf{A}=[\mathbf{A}_{1},\ldots,\mathbf{A}_{m-1}]$, and ${\textbf{x}}=(x_{1},\ldots,x_{m-1})$. It should be noted that except for some special cases, the computation of $A_{m}^{-1}$ requires up to ${\cal O}(n^{3})$ arithmetic operations, which will be a one-time computational cost added on top of the overall computational cost of our algorithm (see Remark 9 below for more discussions).

One may also reformulate problem (3) in the form of (1.1) and directly apply Algorithm 1 to solve it. More specifically, substituting $x_{m}$ with $\mathbf{b}-\mathbf{A}{\textbf{x}}$ in the objective function of (3), we obtain

where $B_{i}=(\textbf{0},\ldots,I,\ldots,\textbf{0})$ with the $i$-th block given a $d_{i}\times d_{i}$ identity matrix and hence $x_{i}=B_{i}{\textbf{x}}$. However, this method will be inefficient since we enlarge the dimension of each $f_{i}$ from $d_{i}$ to $\sum_{i=1}^{m-1}d_{i}$ and as a result, every block has to be updated in each iteration. One may also try to apply a nonconvex randomized block coordinate descent method to solve the above reformulation. However, such methods do not apply to the case when $f_{i}$ are both nonconex and nonsmooth. This motivates us to design the new RapDual method which requires to update only a single block at a time, applies to the case when $f_{i}$ is nonsmooth and achieves an accelerated rate of convergence when $f_{i}$ is smooth.

RaDual (c.f. Algorithm 4) can be viewed as a randomized primal-dual type method. Indeed, by the method of multipliers and Fenchel conjugate duality, we have

Each iteration of Algorithm 4 updates only a randomly selected block $i_{t}$ in (3.49), making it especially favorable when the number of blocks $m$ is large. However, similar difficulty as mentioned in Section 2.1 also appears when we integrate this algorithm with proximal-point type method to yield the final RapDual method in Algorithm 3. Firstly, Algorithm 4 also keeps a few intertwined primal and dual sequences, thus we need to carefully decide the input and output of Algorithm 4 so that information from previous iterations of RapDual is fully used. Secondly, the number of iterations performed by Algorithm 4 to solve each subproblem plays a vital role in the convergence rate of RapDual, which should be carefully predetermined.

We first define an approximate stationary point for problem (1) before establishing the convergence of RapDual.

A stochastic counterpart is one that satisfies

Let $\mathcal{D}^{0}:=f(\bar{\textbf{x}}^{0})+f_{m}(\bar{x}_{m}^{0})-[f({\textbf{x}}^{*})+f_{m}(x_{m}^{*})]$. It can be seen that the total number of primal block updates required to obtain a stochastic $(\epsilon,\delta,\sigma)$-solution can be bounded by

As a comparison, the batch version of this algorithm would update all the $x_{i}^{t}$ for $i=1,\ldots,m$, in (3.49), and thus would require

primal block updates to obtain an $(\epsilon,\delta,\sigma)$-solution of (1). Therefore, the benefit of randomization comes from the difference between $\|\mathbf{A}\|$ and $\bar{A}$. Obviously we always have $\|\mathbf{A}\|>\bar{A}$, and the relative gap between $\|\mathbf{A}\|$ and $\bar{A}$ can be large when all the matrix blocks have close norms. In the case when all the blocks are identical, i.e. $\mathbf{A}_{1}=\mathbf{A}_{2}=\ldots=\mathbf{A}_{m-1}$, we immediately have $\|\mathbf{A}\|=\sqrt{m-1}\bar{A}$, which means that RapDual can potentially save the number of primal block updates by a factor of ${\cal O}(\sqrt{m})$ than its batch counterpart.

In this paper, we assume that $A_{m}^{-1}$ is easily computable. One natural question is whether we can avoid the computation of $A_{m}^{-1}$ by directly solving (1) instead of its reformulation (3). To do so, we can iteratively solve the following saddle-point subproblems in place of the ones in (3.1):

2 Convergence analysis for RapDual

In this section, we first show the convergence of Algorithm 4 for solving the convex multi-block subproblem (5) with

$\psi_{i}(x)-\psi_{i}(y)-\langle\nabla\psi_{i}(y),x-y\rangle\geq\tfrac{\mu}{2}\|x-y\|^{2},\ \ \forall x,y\in X_{i},\quad i=1,\ldots,m-1,$

Some simple relations about the iterations generated by the Algorithm 4 are characterized in the following lemma, and the proof follows directly from the definition of $\hat{\textbf{x}}$ in (3.54), thus has been omitted.

Let $\hat{\emph{{x}}}^{0}=\emph{{x}}^{0}$ and $\hat{\emph{{x}}}^{t}$ for $t=1,\ldots,s$ be defined as follows:

where $\emph{{x}}^{t}$ and $y^{t}$ are obtained from (3.46)-(3.49), then we have

The following lemma 11 builds some connections between the input and output of Algorithm 4 in terms of both primal and dual variables, and the proof can be found in Appendix B.

Let the iterates $\emph{{x}}^{t}$ and $y^{t}$ for $t=1,\ldots,s$ be generated by Algorithm 4 and $(\emph{{x}}^{*},y^{*})$ be a saddle point of (3.1). Assume that the parameters in Algorithm 4 satisfy for all $t=1,\ldots,s-1$

where $\bar{A}=\max_{i\in[m-1]}\|\mathbf{A}_{i}\|$. Then we have

Now we present the main convergence result of Algorithm 4 in Theorem 12, which eliminates the dependence on dual variables and relates directly the successive searching points of RapDual.

It is easy to check that (3.50) and (3.51) satisfy conditions (3.57), (3.58), (3.59) (3.60), and (3.61) when $\mu,\bar{\mu}>0$. Then we have

Therefore, by plugging in those values in (3.50) and (3.51), we have

Since $h(y)$ has $1/\mu$-Lipschitz continuous gradients and is $1/L$-strongly convex, we obtain

Combining (3.63), (3.64) and (3.65), we have

The above theorem shows that subproblem (5) can be solved efficiently by Algorithm 4 with a linear rate of convergence. In fact, we need not solve it too accurately. With a fixed and relatively small number of iterations $s$ Algorithm 4 can still converge, as shown by the following lemma.

Combining (3.66) and (3.2) and noticing that $({\textbf{x}}^{0}_{*},x_{m^{*}}^{0})$ = $(\bar{\textbf{x}}^{0},\bar{x}_{m}^{0})$, we have

Now we are ready to prove the results in Theorem 8 with all the results proved above.

Similarly, due to (3.71) and Lemma 13, we have

Numerical experiments

In this section, we report some preliminary numerical results for both RapGrad and RapDual and demonstrate their potential advantages in Subsection 4.1 and 4.2, respectively.

We consider the least square problem with the smoothly clipped absolute deviation (SCAD) penalty as a testing problem. SCAD has been proved in to be efficient in variable selection. While the original SCAD $p_{\lambda,\gamma}$ defined below does not have smooth gradient at $x=0$, we can bypass this potential problem by using a small positive number $\epsilon$ to obtain a smooth approximation $p_{\lambda,\gamma,\epsilon}$:

where $\gamma>2$, $\lambda>0$, and $\epsilon>0$ are given. Using $p_{\lambda,\gamma,\epsilon}$, our problem of interest is given by

which can be viewed as a special case of problem (1.1) with $f_{i}(x)=\tfrac{1}{2}(a_{i}^{\top}x-b_{i})^{2}+\tfrac{\rho}{2}\sum_{i=1}^{n}p_{\lambda,\gamma,\epsilon}(x_{i})$. Here $a_{i}$ denotes the $i$-th row of $A$. It is easy to see that assumptions (1.2) and (1.3) are satisfied with $\mu=\rho/[2(\gamma-1)]$ and $L=\rho\lambda\epsilon^{-1/2}/2+\max_{1\leq i\leq m}\|a_{i}\|^{2}$. Thus the condition number $L/\mu$ usually dominates $m$.

We test Algorithm 1 (RapGrad), both randomized and deterministic batch versions, on some randomly generated data sets with dimension $m=1000$, $n=100$. Note all entries of matrix $A$ and $20$ uniformly chosen components of $\hat{x}$ are i.i.d from N(0,1). The remaining variables of $\hat{x}$ are set to and $b$ is computed by $b=A\hat{x}$. The parameters used in $p_{\lambda,\gamma,\epsilon}$ are $\epsilon=${10}^{-3}$$,$\lambda=2$,$\gamma=4$and the penalty$\rho$is set to$0.01$. Notice that the$x$-axis represents the number of gradient evaluations divided by$m$, which is counted as number of passes to the dataset, i.e., each iteration of randomized gradient computation is$1/m$pass and the full-gradient computation of counts as$1$ pass. As is shown by Figure 1, randomized version can reduce both objective value and norm of gradient faster than its batch counterpart in terms of number of gradient evaluations. We also observe that there are some zig-zag pattern for the batch algorithm in Figure 1.b), which might have been related some numerical stability issue.

We also compare our RapGrad with the full SVRG in non-convex setting (Algorithm 2 in ) and Accelerated Gradient method (AG) in . Notice that, both RapGrad and SVRG are randomized algorithms, while AG is a deterministic batch method. For the sake of fairness in comparison, all the parameters in the three algorithms mentioned above are set to their theoretical values without any tuning in our experiments. Figure 2 shows our Algorithm 1 not only reduces the function value as well as gradient norm faster than both SVRG and AG.

In fact, our estimate on $s=\left\lceil-\log(6M/5)/\log\alpha\right\rceil$ seems to be too pessimistic, and the subproblems are solved to unnecessarily high accuracy. As a result, some spikes show in Figure 1, which correspond to the occasions when an inner loop for solving the subproblem completes and a new search point is obtained to update the subproblem. Early termination of the inner loops may help to remove those spikes. Reducing the number of inner iterations may not guarantee above mentioned convergence rate theoretically, but may improve the practical performance of RapGrad for this problems in our experiments. From Figure 3 and Figure 4, we can conclude that, by using smaller $s$, our randomized algorithm is able to reduce $f$ and $\|\nabla f\|^{2}$ much faster, whereas the batch version converges faster in terms of the gradient norm.

Inspired by the above experiments, we have an efficient way to tune RapGrad to yield better performance. We first run RapGrad with several different numbers of inner iterations $s^{\prime}$, for instance $s^{\prime}=s$, $s^{\prime}=s/10$, $s^{\prime}=s/100$ , for a fixed number, say $100$, of passes through the dataset, then we use the best $s^{\prime}$ corresponding to the smallest norm of gradient as the actual $s$ for the tuned RapGrad. In the following table, we compare the RapGrad without tuning, tuned RapGrad, SVRG as well as AG on testing problems of different sizes, with stoping criteria $\|\nabla f\|^{2}<${10}^{-10}$$and maximal pass$3\text{\times}{10}^{4}$. The table shows that this simple tuning technique is able to bring huge performance improvement. An interesting observation is that RapGrad without tuning is more likely to outperform SVRG when$n$is large relative to$m$.

2 Nonconvex multi-block optimization

We consider the following compressed sensing problem to test the performance of RapDual:

The numerical experiments is performed on some randomly generated data sets of size $m=1001$, $n=100$. The first $1000$ coefficient matrices $A_{i}$ , $1\leq i\leq 1000$, are of size $100\times 1$ and the last block $A_{1001}$ is an $100\times 100$ identity matrix. $A_{i}$, $1\leq i\leq 1000$ are sparse matrices with sparsity level $0.1$, and all nonzero elements of $A_{i}$ and $200$ uniformly chosen components of $\hat{x}$ are i.i.d from $N(0,1)$. The remaining variables of $\hat{x}$ are set to and $b=\sum_{i=1}^{m}A_{i}\hat{x}_{i}$. The parameters used in $p_{\lambda,\gamma,\epsilon}$ are exactly the same as the first problem, i.e., $\epsilon=${10}^{-3}$$,$\lambda=2$,$\gamma=4$. From Figure 5, we can conclude that the randomized version converges faster than its batch counterpart, in terms of the number of primal block updates required to reduce the objective value and infeasibility. We also compare our randomized algorithm with Algorithm 4 in , with$\rho=L^{2},L^{2}/10,L^{2}/20$. Note$\rho=L^{2}$only guarantees asymptotic convergence. The results in Figure 6 show that our algorithm can reduce the objective value faster than ADMM. As for the feasibility, both algorithms yield solutions that have quite tiny constraint violation. It is worth noting that the converge of ADMM can very much depend on the update order and also that RapDual requires the computation of$A_{m}^{-1}$in its current form. Similar to RapGrad, if we reduce the number of inner iterations$s$per subproblem by a factor of$10$or$20$, we obtain results in Figure 7 and Figure 8. As we can see, the total number of primal block updates needed to yield a good solution, in terms of both objective value and feasibility, can be much smaller when inner loops are terminated early.

Concluding remarks

In this paper, we propose a new randomized accelerated proximal-gradient (RapGrad) method for solving nonconvex finite-sum problems (1.1) and a new randomized primal-dual gradient (RapGrad) method for nonconvex multi-block problems (1), respectively. We demonstrate that for problem (1.1) with large condition number, our RapGrad has much better convergence rate, in terms of dependence on the large number $m$, than the state-of-art nonconvex SVRG or SAGA, as well as accelerated gradient method for nonconvex optimization. Moreover, we show that our RapDual method incorporated with randomization techniques can significantly save the number of primal block updates up to a factor of $\sqrt{m}$ than the deterministic batch methods for solving problems (1). The potential advantages of RapGrad and RapDual are also demonstrated through our preliminary numerical experiments.

We observe that the main ideas of this RapDual, i.e., using proximal point method to tranform the nonconvex multi-block problem into a series of convex subproblems, and using randomized dual method to them, can be applied for solving much more general multi-block optimization problems for which there does not exist an invertible block. In this more general case, the saddle-point subproblems will only be strongly convex in the primal space, but not in the dual space. Therefore, the complexity of solving the subproblem will only be sublinear, and as a consequence, the overall complexity will be much worse than $\mathcal{O}(1/\epsilon)$. It will be interesting to study how much benefit one can obtain by using randomized algorithms in this more general case. We leave this as an interesting topic for future research. It is also worth noting that the proposed RapGrad and RapDual implicitly assume that the parameter $\mu$ is known, or a (tight) upper bound of it can be obtained. While the values of $\mu$ for the problems considered in our numerical experiments can be tightly estimated, it will be interesting to see if we can adaptively estimate the value of $\mu$ in these algorithms applied to solve more general problems.

References

Appendix A Proof of Lemma 4

By convexity of $\psi$ and optimality of $x^{*}$, we have

For notation convenience, let $\Psi(x,z):=\psi(x)-\psi(z)-\langle\nabla\psi(z),x-z\rangle$,

Multiplying each $Q_{t}$ by a non-negative $\gamma_{t}$ and summing them up, we obtain

the first inequality follows from (A.73), (A.74), (2.38), (A) and Lemma 3, and the second inequality is implied by (2.25) and (2.26).

From the relation (2.28) and the fact $x^{-1}=x^{0}$, we have

Now we are ready to bound the last term in (A) as follows:

where (a) follows from the definition $\delta_{t}$ in (A), (b) follows relations (A) and (A) and (c) follows from the fact that $V_{\varphi}(x^{t},x^{t-1})\geq\tfrac{r}{2}\|x^{t}-x^{t-1}\|^{2}$, $\Psi({\underline{x}}^{t-1}_{i_{t}},{\underline{x}}^{t}_{i_{t}})\geq\tfrac{1}{2\hat{L}}\|\nabla\psi_{i_{t}}({\underline{x}}^{t-1}_{i_{t}})-\nabla\psi_{i_{t}}({\underline{x}}^{t}_{i_{t}})\|^{2}$.

By properly regrouping the terms on the right hand side of (LABEL:break), we have

where (a) follows from the simple relation that $b\langle u,v\rangle-a\|v\|^{2}/2\leq b^{2}\|u\|^{2}/(2a)$, $\forall a>0$ and (b) follows from (2.24), (2.27) and (2.28). By using the above inequality, (A) and (A), we obtain

where (a) follows from $\Psi({\underline{x}}^{0}_{i},x^{*})\geq\tfrac{1}{2\hat{L}}\|\nabla\psi_{i_{t}}({\underline{x}}^{0}_{i})-\nabla\psi_{i_{t}}(x^{*})\|^{2}$; (b) follows from the simple relation that $b\langle u,v\rangle-a\|v\|^{2}/2\leq b^{2}\|u\|^{2}/(2a)$, $\forall a>0$ and (c) follows from (2.29), strong convexity of $\psi_{i}$ and Lipschitz continuity of $\nabla\psi_{i}$. This completes the proof. ∎

Appendix B Proof of Lemma 11

For any $t\geq 1$, since $({\textbf{x}}^{*},y^{*})$ is a saddle point of (3.1), we have

For nonnegative $\gamma_{t}$, we further obtain

According to optimality conditions of (3.54) and (3.46) respectively, and strongly convexity of $\psi$ and $h$ we obtain

Combining the above two inequalities with relation (B.82), we have

Applying this and the results (3.55), (3.56) in Lemma 10, we further have

and the second inequality follows from (3.59) and (3.60).

where the second equality follows from (3.58) and the fact that $x^{0}=x^{-1}$.

Since $\langle\mathbf{A}({\textbf{x}}^{t-1}-{\textbf{x}}^{t-2}),y^{t-1}-y^{t}\rangle=\langle\mathbf{A}_{t-1}(x_{i_{t-1}}^{t-1}-x_{i_{t-1}}^{t-2}),y^{t-1}-y^{t}\rangle\leq\|\mathbf{A}_{i_{t-1}}\|\|x_{i_{t-1}}^{t-1}-x_{i_{t-1}}^{t-2}\|\|y^{t}-y^{t-1}\|$ and $V_{h}(y^{t},y^{t-1})\geq\tfrac{\bar{\mu}}{2}\|y^{t}-y^{t-1}\|^{2}$, from (B.84) we have

where (a) follows from regrouping the terms; (b) follows from the definition $\bar{A}=\max_{i\in[m-1]}\|\mathbf{A}_{i}\|$ and the simple relation that $b\langle u,v\rangle-a\|v\|^{2}/2\leq b^{2}\|u\|^{2}/(2a)$, $\forall a>0$; and (c) follows from (3.58) and (3.61).

By combining the relation above with (B), we obtain

In view of (B) and (B.86), we complete the proof. ∎