Efficient random coordinate descent algorithms for large-scale structured nonconvex optimization

Introduction

Coordinate descent methods are among the first algorithms used for solving general minimization problems and are some of the most successful in the large-scale optimization field Ber:99. Roughly speaking, coordinate descent methods are based on the strategy of updating one (block) coordinate of the vector of variables per iteration using some index selection procedure (e.g. cyclic, greedy, random). This often reduces drastically the complexity per iteration and memory requirements, making these methods simple and scalable. There exist numerous papers dealing with the convergence analysis of this type of methods Aus:76; LinLuc:09; Nec:13; NecPat:12; Pow:73; TseYun:09, which confirm the difficulties encountered in proving the convergence for nonconvex and nonsmooth objective functions. For instance, regarding coordinate minimization of nonconvex functions, Powell Pow:73 provided some examples of differentiable functions whose properties lead the algorithm to a closed loop. Also, proving convergence of coordinate descent methods for minimization of nondifferentiable objective functions is challenging Aus:76; FerRic:13. However, for nonconvex and nonsmooth objective functions with certain structure (e.g. composite objective functions) there are available convergence results for coordinate descent methods based on greedy index selection Bec:12; LinLuc:09; TseYun:09 or random index selection LuXia:13. Recently, Nesterov Nes:10 derived complexity results for random coordinate gradient descent methods for solving smooth and convex optimization problems. In RicTac:11 the authors generalized Nesterov’s results to convex problems with composite objective functions. An extensive complexity analysis of coordinate gradient descent methods for solving linearly constrained optimization problems with convex (composite) objective function can be found in Bec:12; Nec:13; NecPat:12; NecNes:12.

In this paper we also consider large-scale nonconvex optimization problems with the objective function consisting of a sum of two terms: one is nonconvex, smooth and given by a black-box oracle, and another is convex but simple and its structure is known. Further, we analyze unconstrained but also singly linearly constrained nonconvex problems. We also assume that the dimension of the problem is so large that traditional optimization methods cannot be directly employed since basic operations, such as the updating of the gradient, are too computationally expensive. These types of problems arise in many fields such as data analysis (classification, text mining) Bon:11; ChaSin:08, systems and control theory (optimal control, pole assignment by static output feedback) FaiMar:12; JudRay:08; NecCli:13; Par:97, machine learning ChaSin:08; NecPat:12; RicTac:13; RicTac:13_2; ShaZha:13; Vap:95 and truss topology design KocOut:06; RicTac:12. The goal of this paper is to analyze several new random coordinate gradient descent methods suited for large-scale nonconvex problems with composite objective function. Recently, after our paper came under review, a variant of random coordinate descent method for solving composite nonconvex problems was also proposed in LuXia:13. For our coordinate descent algorithm, which is designed to minimize unconstrained composite nonconvex objective functions, we prove asymptotic convergence of the generated sequence to stationary points and sublinear rate of convergence in expectation for some optimality measure. We also provide convergence analysis for a coordinate descent method designed for solving singly linearly constrained nonconvex problems and obtain similar results as in the unconstrained case. Note that our analysis is very different from the convex case Nec:13; NecPat:12; NecNes:12; Nes:10; RicTac:11 and is based on the notion of optimality measure and a supermartingale convergence theorem. Furthermore, unlike to other coordinate descent methods for nonconvex problems, our algorithms offer some important advantages, e.g. due to the randomization our algorithms are simpler and are adequate for modern computational architectures. We also present the results of preliminary computational experiments, which confirm the superiority of our methods compared with other algorithms for large-scale nonconvex optimization.

Contribution. The contribution of the paper can be summarized as follows:

For unconstrained problems we propose a 1-random coordinate descent method (1-RCD), that involves at each iteration the solution of an optimization subproblem with respect to only one (block) variable while keeping all others fixed. We show that this solution can be usually computed in closed form (Section 2.2).

For the linearly constrained case we propose a 2-random coordinate descent method (2-RCD), that involves at each iteration the solution of a subproblem depending on two (block) variables while keeping all other variables fixed. We show that in most cases this solution can be found in linear time (Section 3.1).

For each of the algorithms we introduce some optimality measure and devise a convergence analysis using this framework. In particular, for both algorithms, (1-RCD) and (2-RCD), we establish asymptotic convergence of the generated sequences to stationary points (Theorems 2.1 and 3.1) and sublinear rate of convergence for the expected values of the corresponding optimality measures (Theorems 2.2 and 3.2).

Content. The structure of the paper is as follows. In Section 2 we introduce a 1-random coordinate descent algorithm for unconstrained minimization of nonconvex composite functions. Further, we analyze the convergence properties of the algorithm under standard assumptions. In Section 3 we derive a 2-random coordinate descent method for solving singly linearly constrained nonconvex problems and analyze its convergence. In Section 4 we report numerical results on large-scale eigenvalue complementarity problems, which is an important application in control theory.

Unconstrained minimization of composite objective functions

In this section we analyze a variant of random block coordinate gradient descent method, which we call the 1-random coordinate descent method (1-RCD), for solving large-scale unconstrained nonconvex problems with composite objective function. The method involves at each iteration the solution of an optimization subproblem only with respect to one (block) variable while keeping all other variables fixed. After discussing several necessary mathematical preliminaries, we introduce an optimality measure, which will be the basis for the construction and analysis of Algorithm (1-RCD). We establish asymptotic convergence of the sequence generated by Algorithm (1-RCD) to a stationary point and then we show sublinear rate of convergence in expectation for the corresponding optimality measure. For some well-known particular cases of nonconvex objective functions arising frequently in applications, the complexity per iteration of our Algorithm (1-RCD) is of order $\mathcal{O}(\max\limits_{i}n_{i})$.

The problem of interest in this section is the unconstrained nonconvex minimization problem with composite objective function:

where the function $f$ is smooth and $h$ is a convex, separable, nonsmooth function. Since $h$ is nonsmooth, then for any $x\in dom(h)$ we denote by $\partial h(x)$ the subdifferential (set of subgradients) of $h$ at $x$. The smooth and nonsmooth components in the objective function of (1) satisfy the following assumptions:

The function $f$ has block coordinate Lipschitz continuous gradient, i.e. there are constants $L_{i}>0$ such that:

The function $h$ is proper, convex, continuous and block separable:

These assumptions are typical for the coordinate descent framework and the reader can find similar variants in Nec:13; NecPat:12; Nes:10; LuXia:13; TseYun:09. An immediate consequence of Assumption 1 (i) is the following well-known inequality Nes:04:

Based on this quadratic approximation of function $f$ we get the inequality:

Any vector $x^{*}$ satisfying this relation is called a stationary point for nonconvex problem (1).

2 A 1-random coordinate descent algorithm

We analyze a variant of random coordinate descent method suitable for solving large-scale nonconvex problems of the form (1). Let $i\in\{1,\dots,N\}$ be a random variable and $p_{i_{k}}=\text{Pr}(i=i_{k})$ be its probability distribution. Given a point $x$, one block is chosen randomly with respect to the probability distribution $p_{i}$ and the quadratic model (3) derived from the composite objective function is minimized with respect to this block of coordinates (see also Nes:10; RicTac:11). Our method has the following iteration: given an initial point $x_{0}$, then for all $k\geq 0$

where the direction $d_{i_{k}}$ is computed as follows:

Note that the direction $d_{i_{k}}$ is a minimizer of the quadratic approximation model given in (3). Further, from Assumption 1 (ii) we see that $h(x^{k}+U_{i_{k}}s_{i_{k}})=h_{i_{k}}(x_{i_{k}}^{k}+s_{i_{k}})+\sum_{i\not=i_{k}}h_{i}(x_{i}^{k})$ and thus for computing $d_{i_{k}}$ we only need to know the function $h_{i_{k}}(\cdot)$. An important property of our algorithm is that for certain particular cases of function $h$, the complexity per iteration of Algorithm (1-RCD) is very low. In particular, for certain simple functions $h$, very often met in many applications from signal processing, machine learning and optimal control, the direction $d_{i_{k}}$ can be computed in closed form, e.g.:

In this case the direction $d_{i_{k}}$ has the explicit expression:

where $[x]_{[l,\ u]}$ is the orthogonal projection of vector $x$ on box set $[l,\ u]$.

In this case, considering $n=N$, the direction $d_{i_{k}}$ has the explicit expression:

where $t_{i_{k}}=x_{i_{k}}-\frac{1}{L_{i_{k}}}\nabla_{i_{k}}f(x^{k})$.

In these examples the arithmetic complexity of computing the next iterate $x^{k+1}$, once $\nabla_{i_{k}}f(x^{k})$ is known, is of order $\mathcal{O}(n_{i_{k}})$. The reader can find other favorable examples of nonsmooth functions $h$ which preserve the low iteration complexity of Algorithm (1-RCD) (see also NecPat:12; TseYun:09 for other examples). Note that most of the (coordinate descent) methods designed for solving nonconvex problems usually have complexity per iteration at least of order $\mathcal{O}(n)$ (see e.g. TseYun:09, where the authors analyze a greedy coordinate descent method). Coordinate descent methods that have similar complexity per iteration as our random method can be found e.g. in TseYun:09_2, where the index selection is made cyclically (Gauss-Seidel rule). But Algorithm (1-RCD) also offers other important advantages, e.g. due to the randomization the algorithm is adequate for modern computational architectures (e.g distributed and parallel architectures) NecCli:13; RicTac:13.

We assume that the sequence of random variables $i_{0},\dots,i_{k}$ are i.i.d. In the sequel, we use the notation $\xi^{k}$ for the entire history of random index selection

Based on this map, we now introduce an optimality measure which will be the basis for the analysis of Algorithm (1-RCD):

The map $M_{1}(x,L)$ is an optimality measure for optimization problem (1) in the sense that it is positive for all nonstationary points and zero for stationary points (see Lemma 1 below):

Note that $\psi_{L}(s;x)$ is an $1$-strongly convex function in the variable $s$ w.r.t. norm $\lVert\cdot\rVert_{L}$ and thus $d_{L}(x)$ is unique and the following inequality holds:

3 Convergence of Algorithm (1-RCD)

In this section, we analyze the convergence properties of Algorithm (1-RCD). Firstly, we prove the asymptotic convergence of the sequence generated by Algorithm (1-RCD) to stationary points. For proving the asymptotic convergence we use the following supermartingale convergence result of Robbins and Siegmund (see (Pol:87, Lemma 11 on page 50)):

Let $v_{k},u_{k}$ and $\alpha_{k}$ be three sequences of nonnegative random variables such that

where ${\cal F}_{k}$ denotes the collections $v_{0},\dots,v_{k},u_{0},\dots,u_{k}$, $\alpha_{0},\dots,\alpha_{k}$. Then, we have $\lim_{k\to\infty}v_{k}=v$ for a random variable $v\geq 0$ a.s. and $\sum_{k=0}^{\infty}u_{k}<\infty$ a.s.

In the next lemma we prove that Algorithm (1-RCD) is a descent method, i.e. the objective function is nonincreasing along its iterations:

Let $x^{k}$ be the sequence generated by Algorithm (1-RCD) under Assumption 1. Then, the following relation holds:

: From the optimality conditions of subproblem (4) we have that there exists a subgradient $g(x^{k}_{i_{k}}+d_{i_{k}})\in\partial h_{i_{k}}(x^{k}_{i_{k}}+d_{i_{k}})$ such that:

On the other hand, since the function $h_{i_{k}}$ is convex, according to Assumption 1 (ii), the following inequality holds:

Applying the previous two relations in (3) and using the separability of the function $h$, then under Assumption 1 (ii) we have that

Using Lemma 3, we state the following result regarding the asymptotic convergence of Algorithm (1-RCD).

If Assumption 1 holds for the composite objective function $F$ of problem (1) and the sequence $x^{k}$ is generated by Algorithm (1-RCD) using the uniform distribution, then the following statements are valid:

The sequence of random variables $M_{1}(x^{k},L)$ converges to 0 a.s. and the sequence $F(x^{k})$ converges to a random variable $\bar{F}$ a.s.

Any accumulation point of the sequence $x^{k}$ is a stationary point for optimization problem (1).

We now take the expectation conditioned on $\xi^{k-1}$ and note that $i_{k}$ is independent on the past $\xi^{k-1}$, while $x^{k}$ is fully determined by $\xi^{k-1}$. We thus obtain:

Using the supermartingale convergence theorem given in Lemma 2 in the previous inequality, we can ensure that

for a random variable $\theta\geq 0$ and thus $\bar{F}=\theta+F^{*}$. Further, due to almost sure convergence of sequence $F(x^{k})$, it can be easily seen that $\lim\limits_{k\to\infty}F(x^{k})-F(x^{k+1})=0$ a.s. From $x^{k+1}-x^{k}=U_{i_{k}}d_{i_{k}}$ and Lemma 3 we have:

As $\lVert d_{i_{k}}\rVert\to 0$ a.s., we can conclude that the random variable $E[\lVert d_{i_{k}}\rVert|\xi^{k-1}]\to 0$ a.s. or equivalently $M_{1}(x^{k},L)\to 0$ a.s.

(ii) For brevity we assume that the entire sequence $x^{k}$ generated by Algorithm (1-RCD) is convergent. Let $\bar{x}$ be the limit point of the sequence $x^{k}$. In the first part of the theorem we proved that the sequence of random variables $d_{L}(x^{k})$ converges to $0$ a.s. Using the definition of $d_{L}(x^{k})$ we have:

and taking the limit $k\to\infty$ and using Assumption 1 (ii) we get:

This shows that $d_{L}(\bar{x})=0$ is the minimum in subproblem (7) for $x=\bar{x}$ and thus $M_{1}(\bar{x},L)=0$. From Lemma 1 we conclude that $\bar{x}$ is a stationary point for optimization problem (1). ∎

The next theorem proves the convergence rate of the optimality measure $M_{1}(x^{k},L)$ towards $0$ in expectation.

Let $F$ satisfy Assumption 1. Then, the Algorithm (1-RCD) based on the uniform distribution generates a sequence $x^{k}$ satisfying the following convergence rate for the expected values of the optimality measure:

: For simplicity of the exposition we use the following notation: given the current iterate $x$, denote the next iterate $x^{+}=x+U_{i}d_{i}$, where direction $d_{i}$ is given by (4) for some random chosen index $i$ w.r.t. uniform distribution. For brevity, we also adapt the notation of expectation upon the entire history, i.e. $(\phi,\phi^{+},\xi)$ instead of $(\phi^{k},\phi^{k+1},\xi^{k-1})$. From Assumption 1 and inequality (3) we have:

We now take the expectation conditioned on $\xi$:

After rearranging the above expression we get:

Now, by taking the expectation in (10) w.r.t. $\xi$ we obtain:

and then using the $1-$strong convexity property of $\psi_{L}$ we get:

Now coming back to the notation dependent on $k$ and summing w.r.t. the entire history we have:

which leads to the statement of the theorem. ∎

It is important to note that the convergence rate for the Algorithm (1-RCD) given in Theorem 2.2 is typical for the class of first order methods designed for solving nonconvex and nonsmooth optimization problems (see e.g. Nes:07 for more details). Recently, after our paper came under review, a variant of 1-random coordinate descent method for solving composite nonconvex problems was also proposed in LuXia:13. However, the authors in LuXia:13 do not provide complexity results for their algorithm, but only asymptotic convergence in expectation. Note also that our convergence results are different from the convex case Nes:10; RicTac:11, since here we introduce another optimality measure and we use the supermartingale convergence theorem in the analysis.

Furthermore, when the objective function $F$ is smooth and nonconvex, i.e. $h=0$, the first order necessary conditions of optimality become $\nabla f(x^{*})=0$. Also, note that in this case, the optimality measure $M_{1}(x,L)$ is given by: $M_{1}(x,L)=\lVert\nabla f(x)\rVert^{*}_{L}$. An immediate consequence of Theorem 2.2 in this case is the following result:

Let $h=0$ and $f$ satisfy Assumption 1 (i). Then, in this case, the Algorithm (1-RCD) based on the uniform distribution generates a sequence $x^{k}$ satisfying the following convergence rate for the expected values of the norm of the gradients:

Constrained minimization of composite objective functions

In this section we present a variant of random block coordinate gradient descent method for solving large-scale nonconvex optimization problems with composite objective function and a single linear equality constraint:

The function $f$ has 2-block coordinate Lipschitz continuous gradient, i.e. there are constants $L_{ij}>0$ such that:

The function $h$ is proper, convex, continuous and coordinatewise separable:

Note that these assumptions are frequently used in the area of coordinate descent methods for convex minimization, e.g. Bec:12; Nec:13; NecPat:12; NecNes:12; TseYun:09. Based on this assumption the first order necessary optimality conditions become: if $x^{*}$ is a local minimum of (13), then there exists a scalar $\lambda^{*}$ such that:

and then we can bound the function $F$ with the following quadratic expression:

Let $(i,j)$ be a two dimensional random variable, where $i,j\in\{1,\dots,N\}$ with $i\neq j$ and $p_{i_{k}j_{k}}=\text{Pr}((i,j)=(i_{k},j_{k}))$ be its probability distribution. Given a feasible $x$, two blocks are chosen randomly with respect to a given probability distribution $p_{ij}$ and the quadratic model (15) is minimized with respect to these coordinates. Our method has the following iteration: given a feasible initial point $x^{0}$, that is $a^{T}x^{0}=b$, then for all $k\geq 0$

where directions $d_{i_{k}j_{k}}=[d_{i_{k}}^{T}\;d_{j_{k}}^{T}]^{T}$ minimize the quadratic model (15):

In order to analyze the convergence of Algorithm (2-RCD), we introduce an optimality measure:

2 Convergence of Algorithm (2-RCD)

We introduce the notion of elementary vectors for the linear subspace $S=Null(a^{T})$.

An elementary vector $d$ of $S$ is a vector $d\in S$ for which there is no nonzero $d^{\prime}\in S$ conformal to $d$ and $\text{supp}(d^{\prime})\neq\text{supp}(d)$.

We now present some results for elementary vectors and conformal realization, whose proofs can be found in Roc:69; Roc:84; TseYun:09. A particular case of Exercise 10.6 in Roc:84 and an interesting result in Roc:69 provide us the following lemma:

Roc:69; Roc:84 Given $d\in S$, if $d$ is an elementary vector, then $\lvert\text{supp}(d)\rvert\leq 2$. Otherwise, $d$ has a conformal realization $d=d^{1}+\dots+d^{s}$, where $s\geq 2$ and $d^{t}\in S$ are elementary vectors conformal to $d$ for all $t=1,\dots,s$.

An important property of convex and separable functions is given by the following lemma:

where $d^{t}\in S$ are elementary vectors conformal to $d$ for all $t=1,\dots,s$.

If Assumption 2 holds and sequence $x^{k}$ is generated by Algorithm (2-RCD) using the uniform distribution, then the following inequality is valid:

: As in the previous sections, for a simple exposition we drop $k$ from our derivations: e.g. the current point is denoted $x$, next iterate $x^{+}=x+U_{i}d_{i}+U_{j}d_{j}$, where direction $d_{ij}$ is given by Algorithm (2-RCD) for some random selection of pair $(i,j)$ and $\xi$ instead of $\xi^{k-1}$. From the relation (17) and the property of minimizer $d_{ij}$ we have:

Taking expectation in both sides w.r.t. random variable $(i,j)$ conditioned on $\xi$ and recalling that $p_{ij}=\frac{2}{N(N-1)}$, we get:

for all $s_{ij}\in S_{ij}$. We can apply Lemma 6 for coordinatewise separable functions $\lVert\cdot\rVert^{2}$ and $h(\cdot)$ and we obtain:

for all $s_{ij}\in S_{ij}$. From Lemma 5 it follows that any $s\in S$ has a conformal realization defined by $s=\sum_{t}s^{t}$, where the vectors $s^{t}\in S$ are elementary vectors conformal to $s$. Therefore, observing that every elementary vector $s^{t}$ has at most two nonzero blocks, then any vector $s\in S$ can be generated by $s=\sum_{i,j}s_{ij}$, where $s_{ij}\in S$ are conformal to $s$ and have at most two nonzero blocks, i.e. $s_{ij}\in S_{ij}$ for some pair $(i,j)$. Due to conformal property of the vectors $s_{ij}$, the expression $\lVert\sum_{i,j}\sqrt{L_{ij}}s_{ij}\rVert^{2}$ is nondecreasing in the weights $L_{ij}$ and taking in account that $L_{ij}\leq\min\{N\Gamma_{i},N\Gamma_{j}\}$, the previous inequality leads to:

for all $s\in S$. As the last inequality holds for any vector $s\in S$, it also holds for the particular vector $d_{N\Gamma}(x)\in S$:

The main convergence properties of Algorithm (2-RCD) are given in the following theorem:

If Assumption 2 holds for the composite objective function F of problem (13) and the sequence $x^{k}$ is generated by Algorithm (2-RCD) using the uniform distribution, then the following statements are valid:

The sequence of random variables $M_{2}(x^{k},\Gamma)$ converges to 0 a.s. and the sequence $F(x^{k})$ converges to a random variable $\bar{F}$ a.s.

Any accumulation point of the sequence $x^{k}$ is a stationary point for optimization problem (13).

: (i) Using a similar reasoning as in Lemma 3 but for the inequality (15) we can show the following decrease in the objective function for Algorithm (2-RCD) (i.e. Algorithm (2-RCD) is also a descent method):

Further, subtracting $F^{*}$ from both sides, applying expectation conditioned on $\xi^{k-1}$ and then using supermartingale convergence theorem given in Lemma 2 we obtain that $F(x^{k})$ converges to a random variable $\bar{F}$ a.s. for $k\to\infty$. Due to almost sure convergence of sequence $F(x^{k})$, it can be easily seen that $\lim\limits_{k\to\infty}F(x^{k})-F(x^{k+1})=0$ a.s. Moreover, from (19) we have:

As in the previous section, for a simple exposition we drop $k$ from our derivations: e.g. the current point is denoted $x$, next iterate $x^{+}=x+U_{i}d_{i}+U_{j}d_{j}$, where direction $d_{ij}$ is given by Algorithm (2-RCD) for some random selection of pair $(i,j)$ and $\xi$ stands for $\xi^{k-1}$. From Lemma 7, we obtain a sequence which bounds from below $\psi_{N\Gamma}(d_{N\Gamma}(x);x)$ as follows:

On the other hand, from Lemma 5 it follows that any $s\in S$ has a conformal realization defined by $s=\sum_{i,j}s_{ij}$, where $s_{ij}\in S$ are conformal to $s$ and have at most two nonzero blocks, i.e. $s_{ij}\in S_{ij}$ for some pair $(i,j)$. Using now Jensen inequality we derive another sequence which bounds $\psi_{N\Gamma}(d_{N\Gamma}(x);x)$ from above:

Note that $\psi_{N\Gamma}(0;x^{k})=F(x^{k})$ and since both sequences $\psi_{N\Gamma}(0;x^{k})$ and $\psi_{N\Gamma}(d_{N\Gamma}(x^{k});x^{k})$ converge to $\bar{F}$ a.s. for $k\to\infty$, from the above strong convexity relation it follows that the sequence $M_{2}(x^{k};\Gamma)=\lVert d_{N\Gamma}(x^{k})\rVert_{\Gamma}$ converges to $0$ a.s. for $k\to\infty$.

(ii) The proof follows the same ideas as in the proof of Theorem 1 (ii). ∎

We now present the convergence rate for Algorithm (2-RCD).

Let $F$ satisfy Assumption 2. Then, the Algorithm (2-RCD) based on the uniform distribution generates a sequence $x^{k}$ satisfying the following convergence rate for the expected values of the optimality measure:

: Given the current feasible point $x$, denote $x^{+}=x+U_{i}d_{i}+U_{j}d_{j}$ as the next iterate, where direction $(d_{i},d_{j})$ is given by Algorithm (2-RCD) for some random chosen pair $(i,j)$ and we use the notation $(\phi,\phi^{+},\xi)$ instead of $(\phi^{k},\phi^{k+1},\xi^{k-1})$. Based on Lipschitz inequality (15) we derive:

Taking expectation conditioned on $\xi$ in both sides and using Lemma 7 we get:

Taking now expectation w.r.t. $\xi$, we can derive:

where we used the strong convexity property of function $\psi_{N\Gamma}(s;x)$. Now, considering iteration $k$ and summing up with respect to entire history we get:

This inequality leads us to the above result. ∎

3 Constrained minimization of smooth objective functions

where $\nabla f(x)_{\perp}$ is the projection of the gradient vector $\nabla f(x)$ onto the subspace $S$ orthogonal to the vector $a$. Since $\nabla f(x)_{\perp}=\nabla f(x)+\lambda a$, we defined a particular optimality measure:

It is straightforward to see that $Q_{ij}$ is positive semidefinite (notation $Q_{ij}\succeq 0$) and $Q_{ij}a=0$ for all pairs $(i,j)$ with $i\neq j$. Given a probability distribution $p_{ij}$, let us define the matrix:

that is also symmetric and positive semidefinite, since $L_{ij},p_{ij}>0$ for all $(i,j)$. Furthermore, since we consider all possible pairs $(i,j)$, with $i\not=j\in\{1,\dots,N\}$, it can be shown that the matrix $Q$ has an eigenvalue $\nu_{1}(Q)=0$ (which is a simple eigenvalue) with the associated eigenvector $a$. It follows that $\nu_{2}(Q)$ (the second smallest eigenvalue of $Q$) is positive. Since $h=0$, we have $F=f$. Using the same reasoning as in the previous sections we can easily show that the sequence $f(x^{k})$ satisfies the following decrease:

We now give the convergence rate of Algorithm (2-RCD) for this particular case:

Let $h=0$ and $f$ satisfy Assumption 2 (i). Then, Algorithm (2-RCD) based on a general probability distribution $p_{ij}$ generates a sequence $x^{k}$ satisfying the following convergence rate for the expected values of the norm of the projected gradients onto subspace $S$:

As in the previous section, for a simple exposition we drop $k$ from our derivations: e.g. the current point is denoted $x$, and $x^{+}=x+U_{i}d_{i}+U_{j}d_{j}$, where direction $d_{ij}$ is given by Algorithm (2-RCD) for some random selection of pair $(i,j)$. Since $h=0$, we have $F=f$. From (21) we have the following decrease: $f(x^{+})\leq f(x)-\frac{1}{2L_{ij}}\nabla f(x)^{T}Q_{ij}\nabla f(x)$. Taking now expectation conditioned in $\xi$ in this inequality we have:

From the above decomposition of the gradient $\nabla f(x)=\nabla f(x)_{\perp}-\lambda a$ and the observation that $Qa=0$, we conclude that the previous inequality does not change if we replace $\nabla f(x)$ with $\nabla f(x)_{\perp}$:

Note that $\nabla f(x)_{\perp}$ is included in the orthogonal complement of the span of vector $a$, so that the above inequality can be relaxed to:

Coming back to the notation dependent on $k$ and taking expectation in both sides of inequality (22) w.r.t. $\xi^{k-1}$, we have:

Summing w.r.t. the entire history, we obtain the above result. ∎

Note that our convergence proofs given in this section (Theorems 4, 5 and 6) are different from the convex case Nec:13; NecPat:12; NecNes:12, since here we introduce another optimality measure and we use supermartingale convergence theorem in the analysis. It is important to see that the convergence rates for the Algorithm (2-RCD) given in Theorems 3.2 and 3.3 are typical for the class of first order methods designed for solving nonconvex and nonsmotth optimization problems, e.g. in Bec:12; Nes:07 similar results are obtained for other gradient based methods designed to solve nonconvex problems.

Numerical Experiments

In this section we analyze the practical performance of the random coordinate descent methods derived in this paper and compare our algorithms with some recently developed state-of-the-art algorithms from the literature. Coordinate descent methods are one of the most efficient classes of algorithms for large-scale optimization problems. Therefore, we present extensive numerical simulation for large-scale nonconvex problems with dimension ranging from $n=10^{3}$ to $n=10^{7}$. For numerical experiments, we implemented all the algorithms in C code and we performed our tests on a PC with Intel Xeon E5410 CPU and 8 Gb RAM memory.

If matrices A and B are symmetric, then we have symmetric (EiCP). It has been shown in ThiMoe:10 that symmetric (EiCP) is equivalent with finding a stationary point of a generalized Rayleigh quotient on the simplex:

where an upper bound on Lipschitz constant $L_{ij}^{A}$ is given by

: The Hessian of the function $g_{A}(x)$ is given by

Note that $\nabla^{2}_{ij}g_{A}(x)=\frac{2A_{ij}}{x^{T}Ax}-\frac{4(Ax)_{ij}(Ax)_{ij}^{T}}{(x^{T}Ax)^{2}}$. With the same arguments as in ThiMoe:10 we have that: $\lVert\nabla^{2}_{ij}g_{A}(x)\rVert\leq\lVert\frac{2A_{ij}}{x^{T}Ax}\rVert$. From the mean value theorem we obtain:

for any $x,x+s_{ij}\in\Delta_{n}$. Taking norm in both sides of the equality results in:

Note that $\min\limits_{x\in\Delta_{n}}x^{T}Ax>0$ since we have:

and the above result can be easily derived. ∎

Based on the previous notation, the objective function of the logarithmic formulation (23) is given by:

Therefore, the local Lipschitz constants $L_{ij}$ of function $f$ are estimated very easily and numerically cheap as:

where $\mu$ is a parameter chosen in a preliminary stage of the algorithm such that the function $x\mapsto\frac{1}{2}\mu\lVert x\rVert^{2}+\ln(x^{T}Ax)$ is convex. In both algorithms we use the following stopping criterion: $|f(x^{k})-f(x^{k+1})|\leq\epsilon$, where $\epsilon$ is some chosen accuracy. Note that Algorithm (DC) is based on full gradient information and in the application (EiCP) the most computations consists of matrix vector multiplication and a projection onto simplex. When at least one matrix $A$ and $B$ is dense, the computation of the sequence $y^{k}$ is involved, typically $\mathcal{O}(n^{2})$ operations. However, when these matrices are sparse the computation can be reduced to $\mathcal{O}(pn)$ operations, where $p$ is the average number of nonzeros in each row of the matrix $A$ and $B$. Further, there are efficient algorithms for computing the projection onto simplex, e.g. block pivotal principal pivoting algorithm described in JudRay:08, whose arithmetic complexity is of order $\mathcal{O}(n)$. As it appears in practice, the value of parameter $\mu$ is crucial in the rate of convergence of Algorithm (DC). The authors in ThiMoe:10 provide an approximation of $\mu$ that can be computed easily when the matrix $A$ from (23) is positive definite. However, for general copositive matrices (as the case of nonnegative irreducible matrices considered in this paper) one requires the solution of certain NP-hard problem to obtain a good approximation of parameter $\mu$. On the other hand, for our Algorithm (2-RCD) the computation of the Lipschitz constants $L_{ij}$ is very simple and numerically cheap (see previous lemma). Further, for the scalar case (i.e. $n=N$) the complexity per iteration of our method applied to (EiCP) problem is $\mathcal{O}(p)$ in the sparse case.

In Table 1 we compare the two algorithms: (2-CRD) and (DC). We generated random sparse symmetric nonnegative and irreducible matrices of dimension ranging from $n=10^{3}$ to $n=10^{7}$ using the uniform distribution. Each row of the matrices has only $p=10$ nonzero entries. In both algorithms we start from random initial points. In the table we present for each algorithm the final objective function value ($F^{*}$), the number of iterations (iter) and the necessary CPU time (in seconds) for our computer to execute all the iterations. As Algorithm (DC) uses the whole gradient information to obtain the next iterate, we also report for Algorithm (2-RCD) the equivalent number of full-iterations which means the total number of iterations divided by $n/2$ (i.e. the number of iterations groups $x^{0},x^{n/2},...,x^{kn/2}$). Since computing $\mu$ is very difficult for this type of matrices, we try to tune $\mu$ in Algorithm (DC). We have tried four values for $\mu$ ranging from $0.01n$ to $50n$. We have noticed that if $\mu$ is not carefully tuned Algorithm (DC) cannot find the optimal value $f^{*}$ in a reasonable time. Then, after extensive simulations we find an appropriate value for $\mu$ such that Algorithm (DC) produces an accurate approximation of the optimal value. From the table we see that our Algorithm (2-RCD) provides better performance in terms of objective function values and CPU time (in seconds) than Algorithm (DC). We also observe that our algorithm is not sensitive w.r.t. the Lipschitz constants $L_{ij}$ and also w.r.t. the initial point, while Algorithm (DC) is very sensitive to the choice of $\mu$ and the initial point.

Further, in Fig. 1 we plot the evolution of the objective function w.r.t. time for Algorithms (2-RCD) and (DC), in logarithmic scale, on a random (EiCP) problem with dimension $n=5\cdot 10^{5}$ (Algorithm (DC) with parameter left: $\mu=1.42\cdot n$; right: $\mu=50\cdot n$). For a good choice of $\mu$ we see that in the initial phase of Algorithm (DC) the reduction in the objective function is very fast, but while approaching the optimum it slows down. On the other hand, due to the sparsity and randomization our proposed algorithm is faster in numerical implementation than the (DC) scheme.

In Fig. 2 we plot the evolution of CPU time, in logarithmic scale, required for solving the problem w.r.t. the average number of nonzeros entries $p$ in each row of the matrix $A$. We see that for very sparse matrices (i.e. for matrices with relatively small number of nonzeros per row $p\ll n$), our Algorithm (2-RCD) performs faster in terms of CPU time than (DC) method. The main reason is that our method has a simple implementation, does not require the use of other algorithms at each iteration and the arithmetic complexity of an iteration is of order $\mathcal{O}(p)$. On the other hand, Algorithm (DC) is using the block pivotal principal pivoting algorithm described in JudRay:08 at each iteration for projection on simplex and the arithmetic complexity of an iteration is of order $\mathcal{O}(pn)$.

We conclude from the theoretical rate of convergence and the previous numerical results that Algorithms (1-RCD) and (2-RCD) are easier to be implemented and analyzed due to the randomization and the typically very simple iteration. Furthermore, on certain classes of problems with sparsity structure, that appear frequently in many large-scale real applications, the practical complexity of our methods is better than that of some well-known methods from the literature. All these arguments make our algorithms to be competitive in the large-scale nonconvex optimization framework. Moreover, our methods are suited for recently developed computational architectures (e.g., distributed or parallel architectures NecCli:13; RicTac:13).

References