A random coordinate descent algorithm for optimization problems with composite objective function and linear coupled constraints

Introduction

The basic problem of interest in this paper is the following convex minimization problem with composite objective function:

Linearly constrained optimization problems with composite objective function arise in many applications such as compressive sensing CanRom:06, image processing CheDon:01, truss topology design NesShp:12, distributed control NecNed:11, support vector machines TseYun:07, traffic equilibrium and network flow problems Ber:03 and many other areas. For problems of moderate size there exist many iterative algorithms such as Newton, quasi-Newton or projected gradient methods DaiFle:06; FerMun:03; LinLuc:09. However, the problems that we consider in this paper have the following features: the dimension of the optimization variables is very large such that usual methods based on full gradient computations are prohibitive. Moreover, the incomplete structure of information that may appear when the data are distributed in space and time, or when there exists lack of physical memory and enormous complexity of the gradient update can also be an obstacle for full gradient computations. In this case, it appears that a reasonable approach to solving problem (1) is to use (block) coordinate descent methods. These methods were among the first optimization methods studied in literature Ber:99. The main differences between all variants of coordinate descent methods consist of the criterion of choosing at each iteration the coordinate over which we minimize our objective function and the complexity of this choice. Two classical criteria, used often in these algorithms, are the cyclic and the greedy (e.g., Gauss-Southwell) coordinate search, which significantly differ by the amount of computations required to choose the appropriate index. The rate of convergence of cyclic coordinate search methods has been determined recently in BecTet:12; SahTew:12. Also, for coordinate descent methods based on the Gauss-Southwell rule, the convergence rate is given in TseYun:06; TseYun:07; TseYun:09. Another interesting approach is based on random coordinate descent, where the coordinate search is random. Recent complexity results on random coordinate descent methods were obtained by Nesterov in Nes:10 for smooth convex functions. The extension to composite objective functions was given in RicTak:11; RicTak:12 and for the grouped Lasso problem in QinSch:10. However, all these papers studied optimization models where the constraint set is decoupled (i.e., characterized by Cartesian product). The rate analysis of a random coordinate descent method for linearly coupled constrained optimization problems with smooth objective function was developed in NecNes:12.

The paper is organized as follows. In order to present our main results, we introduce some notations and assumptions for problem (1) in Section 1.1. In Section 2 we present the new random coordinate descent (RCD) algorithm. The main results of the paper can be found in Section 3, where we derive the rate of convergence in expectation, probability and for the strongly convex case. In Section 4 we generalize the algorithm and extend the previous results to a more general model. We also analyze its complexity and compare it with other methods from the literature, in particular the coordinate descent method of Tseng TseYun:09 in Section 5. Finally, we test the practical efficiency of our algorithm through extensive numerical experiments in Section 6.

We assume that the entire space dimension is decomposable into $N$ blocks:

We denote by $U_{i}$ the blocks of the identity matrix:

For model (1) we make the following assumptions:

The smooth and nonsmooth parts of the objective function in optimization model (1) satisfy the following properties:

Function $f$ is convex and has block-coordinate Lipschitz continuous gradient:

The nonsmooth function $h$ is convex and coordinatewise separable.

where $\lambda>0$. Often, a large $\lambda$ factor induces sparsity in the solution of optimization problem (1). Note that the function $h$ in (2) belongs to the general class of coordinatewise separable piecewise linear/quadratic functions with $\mathcal{O}(1)$ pieces. Another special case is the box indicator function, i.e.:

Adding box constraints to a quadratic objective function $f$ in (1) leads e.g., to support vector machine (SVM) problems ChaLin:11; TseYun:07. The reader can easily find many other examples of function $h$ satisfying Assumption 1 $(ii)$.

Based on Assumption 1 $(i)$, the following inequality can be derived Nes:04:

Note that these norms satisfy the Cauchy-Schwartz inequality:

Based on Assumption 1 $(ii)$ we can derive from (4) the following result:

Let function $f$ be convex and satisfy Assumption 1. Then, the function $f$ has componentwise Lipschitz continuous gradient w.r.t. every pair $(i,j)$, i.e.:

where we define $L_{ij}^{\alpha}=L_{i}^{1-\alpha}+L_{j}^{1-\alpha}$.

where in the third inequality we used that $\alpha a+(1-\alpha)b\leq\max\{a,b\}$ for all $\alpha\in$. Now, note that the function $g_{1}(y_{ij})=f(x+y_{ij}-x_{ij})-f(x)-\langle\nabla f(x),y_{ij}-x_{ij}\rangle$ satisfies the Assumption 1 $(i)$. If we apply the above inequality to $g_{1}(y_{ij})$ we get the following relation:

On the other hand, applying the same inequality to $g_{2}(x_{ij})=f(x)-f(x+y_{ij}-x_{ij})+\langle\nabla f(x+y_{ij}-x_{ij}),y_{ij}-x_{ij}\rangle$, which also satisfies Assumption 1 $(i)$, we have:

Further, denoting $s_{ij}=y_{ij}-x_{ij}$ and adding up the resulting inequalities we get:

It is straightforward to see that we can obtain from Lemma 1 the following inequality (see also Nes:04):

Random coordinate descent algorithm

In this section we introduce a variant of Random Coordinate Descent (RCD) method for solving problem (1) that performs a minimization step with respect to two block variables at each iteration. The coupling constraint (that is, the weighted sum constraint $a^{T}x=0$) prevents the development of an algorithm that performs a minimization with respect to only one variable at each iteration. We will therefore be interested in the restriction of the objective function $f$ on feasible directions consisting of at least two nonzero (block) components.

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 probability distribution $p_{ij}$ and a quadratic model derived from the composite objective function 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}=0$, then for all $k\geq 0$

where the directions $d_{i_{k}}$ and $d_{j_{k}}$ are chosen as follows: if we use for simplicity the notation $(i,j)$ instead of $({i_{k}},{j_{k}})$, the direction $d_{ij}=[d_{i}^{T}\;d_{j}^{T}]^{T}$ is given by

Note that for the scalar case (i.e., $N=n$) and $h$ given by (2) or (3), the direction $d_{ij}$ in (6) can be computed in closed form. For the block case (i.e., $n_{i}>1$ for all $i$) and if $h$ is a coordinatewise separable, strictly convex and piece-wise linear/quadratic function with $\mathcal{O}(1)$ pieces (e.g., $h$ given by (2)), there are algorithms for solving the above subproblem in linear-time (i.e., $\mathcal{O}(n_{i}+n_{j})$ operations) TseYun:09. Also for $h$ given by (3), there exist in the literature algorithms for solving the subproblem (6) with overall complexity $\mathcal{O}(n_{i}+n_{j})$ BerKov:93; Kiw:07.

In algorithm (RCD) we consider $(i,j)=(j,i)$ and $i\neq j$. Moreover, we know that the complexity of choosing randomly a pair $(i,j)$ with a uniform probability distribution requires $\mathcal{O}(1)$ operations. ∎

We assume that random variables $(i_{k},j_{k})_{k\geq 0}$ are i.i.d. In the sequel, we use notation $\eta^{k}$ for the entire history of random pair choices and $\phi^{k}$ for the expected value of the objective function w.r.t. $\eta^{k}$, i.e.:

We briefly review some well-known methods from the literature for solving the optimization model (1). In TseYun:06; TseYun:07; TseYun:09 Tseng studied optimization problems in the form (1) and developed a (block) coordinate gradient descent(CGD) method based on the Gauss-Southwell choice rule. The main requirement for the (CGD) iteration is the solution of the following problem: given a feasible $x$ and a working set of indexes $\mathcal{J}$, the update direction is defined by

In TseYun:09, the authors proved for the particular case when function $h$ is piece-wise linear/quadratic with $\mathcal{O}(1)$ pieces that an $\epsilon$-optimal solution is attained in $\mathcal{O}(\frac{nLR_{0}^{2}}{\epsilon})$ iterations, where $R_{0}$ denotes the Euclidean distance from the initial point to some optimal solution. Also, in TseYun:09 the authors derive estimates of order $\mathcal{O}(n)$ on the computational complexity of each iteration for this choice of $h$.

Furthermore, for a quadratic function $f$ and a box indicator function $h$ (e.g., support vector machine (SVM) applications) one of the first decomposition approaches developed similar to (RCD) is Sequential Minimal Optimization (SMO) Pla:99. SMO consists of choosing at each iteration two scalar coordinates with respect to some heuristic rule based on KKT conditions and solving the small QP subproblem obtained through the decomposition process. However, the rate of convergence is not provided for the SMO algorithm. But the numerical experiments show that the method is very efficient in practice due to the closed form solution of the QP subproblem. List and Simon LisSim:05 proposed a variant of block coordinate descent method for which an arithmetic complexity of order $\mathcal{O}(\frac{n^{2}LR_{0}^{2}}{\epsilon})$ is proved on a quadratic model with a box indicator function and generalized linear constraints. Later, Hush et al. HusKel:06 presented a more practical decomposition method which attains the same complexity as the previous methods.

A random coordinate descent algorithm for model (1) with $a=0$ and $h$ being the indicator function for a Cartesian product of sets was analyzed by Nesterov in Nes:10. The generalization of this algorithm to composite objective functions has been studied in QinSch:10; RicTak:11. However, none of these papers studied the application of coordinate descent algorithms to linearly coupled constrained optimization models. A similar random coordinate descent algorithm as the (RCD) method described in the present paper, for optimization problems with smooth objective and linearly coupled constraints, has been developed and analyzed by Necoara et al. in NecNes:12. We further extend these results to linearly constrained composite objective function optimization and provide in the sequel the convergence rate analysis for the previously presented variant of the (RCD) method (see Algorithm 1 (RCD)).

Convergence results

In the following subsections we derive the convergence rate of Algorithm 1 (RCD) for composite optimization model (1) in expectation, probability and for strongly convex functions.

In this section we study the rate of convergence in expectation of algorithm (RCD). We consider uniform probability distribution, i.e., the event of choosing a pair $(i,j)$ can occur with probability:

since we assume that $(i,j)=(j,i)$ and $i\neq j\in\{1,\dots,N\}$ (see Remark 1 (ii)). In order to provide the convergence rate of our algorithm, first we have to define the conformal realization of a vector introduced in Roc:67; Roc:84.

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

Based on Exercise 10.6 in Roc:84 we state the following lemma:

Roc:84 Given $d\in Null(A)$, if $d$ is an elementary vector, then $\left|supp(d)\right|\leq rank(A)+1\leq m+1$. Otherwise, $d$ has a conformal realization:

where $s\geq 1$ and $d^{t}\in\text{Null}(A)$ are elementary vectors conformal to $d$ for all $t=1,\dots,s$.

For the scalar case, i.e., $N=n$ and $m=1$, the method provided in TseYun:09 finds a conformal realization with dimension $s\leq|\text{supp}(d)|-1$ within $\mathcal{O}(n)$ operations. We observe that elementary vectors $d^{t}$ in Lemma 2 for the case $m=1$ (i.e., $A=a^{T}$) have at most $2$ nonzero components.

Our convergence analysis is based on the following lemma, whose proof can be found in (TseYun:09, Lemma 6.1):

For the simplicity of the analysis we introduce the following linear subspaces:

A simplified update rule of algorithm (RCD) is expressed as:

We denote by $F^{*}$ and $X^{*}$ the optimal value and the optimal solution set for problem (1), respectively. We also introduce the maximal residual defined in terms of the norm $\|\cdot\|_{\alpha}$:

which measures the size of the level set of $F$ given by $x^{0}$. We assume that this distance is finite for the initial iterate $x^{0}$.

Now, we prove the main result of this section:

Let $F$ satisfy Assumption 1. Then, the random coordinate descent algorithm (RCD) based on the uniform distribution generates a sequence $x^{k}$ satisfying the following convergence rate for the expected values of the objective function:

For simplicity, we drop the index $k$ and use instead of $(i_{k},j_{k})$ and $x^{k}$ the notation $(i,j)$ and $x$, respectively. Based on (5) we derive:

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

for all possible $s_{ij}\in S_{ij}$ and pairs $(i,j)$ with $i\neq j\in\{1,\dots,N\}$.

Based on Lemma 2 for $m=1$, it follows that any $d\in S$ has a conformal realization defined by $d=\sum\limits_{t=1}^{s}d^{t}$, where the vectors $d^{t}\in S$ are conformal to $d$ and have only two nonzero components. Thus, for any $t=1,\dots,s$ there is a pair $(i,j)$ such that $d^{t}\in S_{ij}$. Therefore, for any $d\in S$ we can choose an appropriate set of pairs $(i,j)$ and vectors $s_{ij}^{d}\in S_{ij}$ conformal to $d$ such that $d=\sum\limits_{i,j}s_{ij}^{d}$. As we have seen, the above chain of relations in (8) holds for any set of pairs $(i,j)$ and vectors $s_{ij}\in S_{ij}$. Therefore, it also holds for the set of pairs $(i,j)$ and vectors $s_{ij}^{d}$ such that $d=\sum\limits_{i,j}s_{ij}^{d}$. In conclusion, we have from (8) that:

for all $d\in S$. Moreover, observing that $L_{ij}^{\alpha}\leq 2L^{1-\alpha}$ and applying Lemma 3 in the previous inequality for coordinatewise separable functions $\left\|\cdot\right\|^{2}_{\alpha}$ and $h(\cdot)$, we obtain:

Based on this choice and using similar reasoning as in Nes:07; RicTak:11 for proving the convergence rate of gradient type methods for composite objective functions, we derive the following:

where in the first inequality we used the convexity of $f$ while in the second and third inequalities we used basic optimization arguments. Therefore, at each iteration $k$ the following inequality holds:

Taking expectation with respect to $\eta_{k}$ and using convexity properties we get:

Further, if we denote $\Delta^{k}=\phi^{k}-F^{*}$ and $\gamma=N(N-1)L^{1-\alpha}R^{2}_{\alpha}$ we get:

Dividing both sides with $\Delta^{k}\Delta^{k+1}>0$ and using the fact that $\Delta^{k+1}\leq\Delta^{k}$ we get:

Finally, summing up from $0,\dots,k$ we easily get the above convergence rate. ∎

Let us analyze the convergence rate of our method for the two most common cases of the extended norm introduced in this section: w.r.t. extended Euclidean norm $\left\|\cdot\right\|_{0}$ (i.e., $\alpha=0$) and norm $\left\|\cdot\right\|_{1}$ (i.e., $\alpha=1$). Recall that the norm $\left\|\cdot\right\|_{1}$ is defined by:

Under the same assumptions of Theorem 3.1, the algorithm (RCD) generates a sequence $x^{k}$ such that the expected values of the objective function satisfy the following convergence rates for $\alpha=0$ and $\alpha=1$:

We usually have $R_{1}^{2}\leq LR_{0}^{2}$ and this shows the advantages that the general norm $\left\|\cdot\right\|_{\alpha}$ has over the Euclidean norm. Indeed, if we denote by $r_{i}^{2}=\max_{x}\{\max_{x^{*}\in X^{*}}\left\|x_{i}-x^{*}_{i}\right\|^{2}:\;F(x)\leq F(x^{0})\}$, then we can provide upper bounds on $R_{1}^{2}\leq\sum_{i=1}^{N}L_{i}r_{i}^{2}$ and $R_{0}^{2}\leq\sum_{i=1}^{N}r_{i}^{2}$. Clearly, the following inequality is valid:

and the inequality holds with equality only for $L_{i}=L$ for all $i=1,\dots,N$. We recall that $L=\max_{i}L_{i}$. Therefore, in the majority of cases the estimate for the rate of convergence based on norm $\left\|\cdot\right\|_{1}$ is much better than that based on the Euclidean norm $\left\|\cdot\right\|_{0}$.

2 Convergence for strongly convex functions

Now, we assume that the objective function in (1) is $\sigma_{\alpha}$-strongly convex with respect to norm $\left\|\cdot\right\|_{\alpha}$, i.e.:

where $F^{\prime}(y)$ denotes some subgradient of $F$ at $y$. Note that if the function $f$ is $\sigma$-strongly convex w.r.t. extended Euclidean norm, then we can remark that it is also $\sigma_{\alpha}$-strongly convex function w.r.t. norm $\left\|\cdot\right\|_{\alpha}$ and the following relation between the strong convexity constants holds:

Taking $y=x^{*}$ in (11) and from optimality conditions $\langle F^{\prime}(x^{*}),x-x^{*}\rangle\geq 0$ for all $x\in S$ we obtain:

Next, we state the convergence result of our algorithm (RCD) for solving the problem (1) with $\sigma_{\alpha}$-strongly convex objective w.r.t. norm $\left\|\cdot\right\|_{\alpha}$.

Under the assumptions of Theorem 3.1, let $F$ be also $\sigma_{\alpha}$-strongly convex w.r.t. $\left\|\cdot\right\|_{\alpha}$. For the sequence $x^{k}$ generated by algorithm (RCD) we have the following rate of convergence of the expected values of the objective function:

Then, using similar derivation as in Theorem 1 we have:

where the last inequality results from (12). The statement of the theorem is obtained by noting that $\beta^{*}=\min\{1,\frac{\sigma_{\alpha}}{4L^{1-\alpha}}\}$ and the following subcases:

If $\beta^{*}=\frac{\sigma_{\alpha}}{4L^{1-\alpha}}$ and we take the expectation w.r.t. $\eta^{k}$ we get:

if $\beta^{*}=1$ and we take the expectation w.r.t. $\eta^{k}$ we get:

3 Convergence in probability

Further, we establish some bounds on the required number of iterations for which the generated sequence $x^{k}$ attains $\epsilon$-accuracy with prespecified probability. In order to prove this result we use Theorem 1 from RicTak:11 and for a clear understanding we present it bellow.

RicTak:11 Let $\xi^{0}>0$ be a constant, $0<\epsilon<\xi^{0}$ and consider a nonnegative nonincreasing sequence of (discrete) random variables $\{\xi^{k}\}_{k\geq 0}$ with one of the following properties:

$E[\xi^{k+1}|\xi^{k}]\leq\xi^{k}-\frac{(\xi^{k})^{2}}{c}$ for all $k$, where $c>0$ is a constant,

$E[\xi^{k+1}|\xi^{k}]\leq\left(1-\frac{1}{c}\right)\xi^{k}$ for all $k$ such that $\xi^{k}\geq\epsilon$, where $c>1$ is a constant.

Then, for some confidence level $\rho\in(0,1)$ we have in probability that:

for a number $K$ of iterations which satisfies

Based on this lemma we can state the following rate of convergence in probability:

Let $F$ be a $\sigma_{\alpha}$-strongly convex function satisfying Assumption 1 and $\rho>0$ be the confidence level. Then, the sequence $x^{k}$ generated by algorithm (RCD) using uniform distribution satisfies the following rate of convergence in probability of the expected values of the objective function:

where $\gamma=\begin{cases}1-\frac{\sigma_{\alpha}}{8L^{1-\alpha}},&if\ \sigma_{\alpha}\leq 4L^{1-\alpha}\\ \frac{2L^{1-\alpha}}{\sigma_{\alpha}},&\ \text{otherwise}.\end{cases}$

Based on relation (10), we note that taking $\xi^{k}$ as $\xi^{k}=\phi^{k}-F^{*}$, the property $(1)$ of Lemma 4 holds and thus we get the first part of our result. Relations (13) and (14) in the strongly convex case are similar instances of property $(2)$ in Theorem 4 from which we get the second part of the result. ∎

Generalization

In this section we study the optimization problem (1), but with general linearly coupling constraints:

where the direction $d_{\mathcal{N}_{k}}$ is chosen as follows:

We can easily see that the linearly coupling constraints $Ax=0$ prevent the development of an algorithm that performs at each iteration a minimization with respect to less than $m+1$ coordinates. Therefore we are interested in the class of iteration updates which restricts the objective function on feasible directions that consist of at least $m+1$ (block) components.

Let $F$ satisfy Assumption 1. Then, the random coordinate descent algorithm (RCD)N that chooses uniformly at each iteration $m+1$ blocks generates a sequence $x^{k}$ satisfying the following rate of convergence for the expected values of the objective function:

The proof is similar to that of Theorem 3.1 and we omit it here for brevity.

Complexity analysis

In this section we analyze the total complexity (arithmetic complexity Nes:04) of algorithm (RCD) based on extended Euclidean norm for optimization problem (1) and compare it with other complexity estimates. Tseng presented in TseYun:09 the first complexity bounds for the (CGD) method applied to our optimization problem (1). Up to our knowledge there are no other complexity results for coordinate descent methods on the general optimization model (1).

Note that the algorithm (RCD) has an overall complexity w.r.t. extended Euclidean norm given by:

where $\mathcal{O}(i_{RCD})$ is the complexity per iteration of algorithm (RCD). On the other hand, algorithm (CGD) has the following complexity estimate:

where $\mathcal{O}(i_{CGD})$ is the iteration complexity of algorithm (CGD). Based on the particularities and computational effort of each method, we will show in the sequel that for some optimization models arising in real-world applications the arithmetic complexity of (RCD) method is lower than that of (CGD) method. For certain instances of problem (1) we have that the computation of the coordinate directional derivative of the smooth component of the objective function is much more simpler than the function evaluation or directional derivative along an arbitrary direction. Note that the iteration of algorithm (RCD) uses only a small number of coordinate directional derivatives of the smooth part of the objective, in contrast with the (CGD) iteration which requires the full gradient. Thus, we estimate the arithmetic complexity of these two methods applied to a class of optimization problems containing instances for which the directional derivative of objective function can be computed cheaply. We recall that the process of choosing a uniformly random pair $(i,j)$ in our method requires $\mathcal{O}(1)$ operations.

Let us structure a general coordinate descent iteration in two phases: Phase 1: Gather first-order information to form a quadratic approximation of the original optimization problem. Phase 2: Solve a quadratic optimization problem using data acquired at Phase 1 and update the current vector. Both algorithms (RCD) and (CGD) share this structure but, as we will see, there is a gap between computational complexities. We analyze the following example:

Further, we estimate the iteration complexity of the algorithms (RCD) and (CGD). Given a feasible $x$, from the expression

we note that if the residual $r(x)=Zx$ is already known, then the computation of $\nabla_{i}f(x)$ requires $\mathcal{O}(p)$ operations. We consider that the dimension $n_{i}$ of each block is of order $\mathcal{O}(\frac{n}{N})$. Thus, the (RCD) method updates the current point $x$ on $\mathcal{O}(\frac{n}{N})$ coordinates and summing up with the computation of the new residual $r(x^{+})=Zx^{+}$, which in this case requires $\mathcal{O}(\frac{pn}{N})$ operations, we conclude that up to this stage, the iteration of (RCD) method has numerical complexity $\mathcal{O}(\frac{pn}{N})$. However, the (CGD) method requires the computation of the full gradient for which are necessary $\mathcal{O}(np)$ operations. As a preliminary conclusion, Phase 1 has the following complexity regarding the two algorithms:

Suppose now that for a given $x$, the blocks $(\nabla_{i}f(x),\nabla_{j}f(x))$ are known for (RCD) method or the entire gradient vector $\nabla f(x)$ is available for (CGD) method within previous computed complexities, then the second phase requires the finding of an update direction with respect to each method. For the general linearly constrained model (1), evaluating the iteration complexity of both algorithms can be a difficult task. Since in TseYun:09 Tseng provided an explicit total computational complexity for the cases when the nonsmooth part of the objective function $h$ is separable and piece-wise linear/quadratic with $\mathcal{O}(1)$ pieces, for clarity of the comparison we also analyze the particular setting when $h$ is a box indicator function as given in equation (3). For algorithm (RCD) with $\alpha=0$, at each iteration, we require the solution of the following problem (see (3)):

It is shown in Kiw:07 that problem (18) can be solved in $\mathcal{O}(n_{i}+n_{j})$ operations. However, in the scalar case (i.e., $N=n$) problem (18) can solved in closed form. Therefore, Phase 2 of algorithm (RCD) requires $\mathcal{O}(\frac{n}{N})$ operations. Finally, we estimate for algorithm (RCD) the total arithmetic complexity in terms of the number of blocks $N$ as:

On the other hand, due to the Gauss-Southwell rule, the (CGD) method requires at each iteration the solution of a quadratic knapsack problem of dimension $n$. It is argued in Kiw:07 that for solving the quadratic knapsack problem we need $\mathcal{O}(n)$ operations. In conclusion, the Gauss-Southwell procedure in algorithm (CGD) requires the conformal realization of the solution of a continuous knapsack problem and the selection of a “good” set of blocks $\mathcal{J}$. This last process has a different cost depending on $m$. Overall, we estimate the total complexity of algorithm (CGD) for one equality constraint, $m=1$, as:

First, we note that in the case $m=1$ and $N<<n$ (i.e., the block case) algorithm (RCD) has better arithmetic complexity than algorithm (CGD) and previously mentioned block-coordinate methods HusKel:06; LisSim:05 (see Table 1). When $m=1$ and $N=n$ (i.e., the scalar case), by substitution in the above expressions from Table 1, we have a total complexity for algorithm (RCD) comparable to the complexity of algorithm (CGD) and the algorithms from HusKel:06; LisSim:05.

On the other hand, the complexity of choosing a random pair $(i,j)$ in algorithm (RCD) is very low, i.e., we need $\mathcal{O}(1)$ operations. Thus, choosing the working pair $(i,j)$ in our algorithm (RCD) is much simpler than choosing the working set $\mathcal{J}$ within the Gauss-Southwell rule for algorithm (CGD) which assumes the following steps: first, compute the projected gradient direction and second, find the conformal realization of computed direction; the overall complexity of these two steps being $\mathcal{O}(n)$. In conclusion, the algorithm (RCD) has a very simple implementation due to simplicity of the random choice for the working pair and a low complexity per iteration.

For the case $m=2$ the algorithm (RCD) needs in Phase 1 to compute coordinate directional derivatives with complexity $\mathcal{O}(\frac{pn}{N})$ and in Phase 2 to find the solution of a 3-block dimensional problem of the same structure as (18) with complexity $\mathcal{O}(\frac{n}{N})$. Therefore, the iteration complexity of the (RCD) method in this case is still $\mathcal{O}(\frac{pn}{N})$. On the other hand, the iteration complexity of the algorithm (CGD) for $m=2$ is given by $\mathcal{O}(pn+n\log n)$ TseYun:09.

For $m>2$, the complexity of Phase 1 at each iteration of our method still requires $\mathcal{O}(\frac{pn}{N})$ operations and the complexity of Phase 2 is $\mathcal{O}(\frac{mn}{N})$, while in the (CGD) method the iteration complexity is $\mathcal{O}(m^{3}n^{2})$ TseYun:09.

For the case $m>1$, a comparison between arithmetic complexities of algorithms (RCD) and (CGD) is provided in Table 2. We see from this table that depending on the values of $n,m$ and $N$, the arithmetic complexity of (RCD) method can be better or worse than that of the (CGD) method.

We conclude from the rate of convergence and the previous complexity analysis that algorithm (RCD) is easier to be implemented and analyzed due to the randomization and the typically very simple iteration. Moreover, on certain classes of problems with sparsity structure, that appear frequently in many large-scale real applications, the arithmetic complexity of (RCD) method is better than that of some well-known methods from the literature. All these arguments make the algorithm (RCD) to be competitive in the composite optimization framework. Moreover, the (RCD) method is suited for recently developed computational architectures (e.g., distributed or parallel architectures).

Numerical Experiments

We have implemented all the algorithms in C-code and the experiments were run on a PC with an Intel Xeon E5410 CPU and 8 GB RAM memory. In all algorithms we considered the scalar case, i.e., $N=n$ and we worked with the extended Euclidean norm ($\alpha=0$). In our applications the smooth part $f$ of the composite objective function is of the form (17). The coordinate directional derivative at the current point for algorithm (RCD) $\nabla_{i}f(x)=\langle z_{i},Zx\rangle+q_{i}$ is computed efficiently by knowing at each iteration the residual $r(x)=Zx$. For the (CGD) method, the working set is chosen accordingly to Section 6 in TseYun:07. Therefore, the entire gradient at the current point, $\nabla f(x)=Z^{T}Zx+q$, is required, which is computed efficiently using the residual $r(x)=Zx$. For gradient and residual computations we used an efficient sparse matrix-vector multiplication procedure. We coded the standard (CGD) method presented in TseYun:09 and we have not used any heuristics recommended by Tseng in TseYun:07, e.g., the “3-pair” heuristic technique. The direction $d_{ij}$ at the current point from subproblem (6) for algorithm (RCD) is computed in closed form for all three applications considered in this section. For computing the direction $d_{H}(x;\mathcal{J})$ at the current point from subproblem (7) in the (CGD) method for the first two applications we coded the algorithm from Kiw:07 for solving quadratic knapsack problems of the form (18) that has linear time complexity. For the second application, the direction at the current point for algorithm (GM) is computed using a linear time simplex projection algorithm introduced in JudRay:08. For the third application, we used the equivalent formulation of the subproblem (7) given in TseYun:09, obtaining for both algorithms (CGD) and (GM) an iteration which requires the solution of some double size quadratic knapsack problem of the form (18).

In the following tables we present for each algorithm the final objective function value (obj), the number of iterations (iter) and the necessary CPU time for our computer to execute all the iterations. As the algorithms (CGD), LIBSVM and (GM) use the whole gradient information to obtain the working set and to find the direction at the current point, we also report for the algorithm (RCD) the equivalent number of full-iterations which means the total number of iterations divided by $\frac{n}{2}$ (i.e., the number of iterations groups $x^{0},x^{n/2},\dots,x^{kn/2}$).

In order to better understand the practical performance of our method, we have tested the algorithms (RCD), (CGD) and LIBSVM on two-class data classification problems with linear kernel, which is a well-known real-world application that can be posed as a large-scale optimization problem in the form (1) with a sparsity structure. In this section, we describe our implementation of algorithms (RCD), (CGD) TseYun:07 and LIBSVM ChaLin:11 and report the numerical results on different test problems. Note that linear SVM is a technique mainly used for text classification, which can be formulated as the following optimization problem:

We report in Table $3$ the results for algorithms (RCD), (CGD) and LIBSVM implemented in the scalar case, i.e., $N=n$. The data used for the experiments can be found on the LIBSVM webpage (http://www.csie.ntu.edu.tw/cjlin/libsvmtools/ datasets/). For problems with very large dimensions, we generated the data randomly (see “test1” and “test2”) such that the nonzero elements of $Z$ fit into the available memory of our computer. For each algorithm we present the final objective function value (obj), the number of iterations (iter) and the necessary CPU time (in minutes) for our computer to execute all the iterations. For the algorithm (RCD) we report the equivalent number of full-iterations, that is the number of iterations groups $x^{0},x^{n/2},\dots,x^{kn/2}$. On small test problems we observe that LIBSVM outperforms algorithms (RCD) and (CGD), but we still have that the CPU time for algorithm (RCD) does not exceed $30$ min, while algorithm (CGD) performs much worse. On the other hand, on large-scale problems the algorithm (RCD) has the best behavior among the three tested algorithms (within a factor of $10$). For very large problems ($n\geq 10^{6}$), LIBSVM has not returned any result within $10$ hours.

For the block case (i.e., $N\leq n$), we have plotted for algorithm (RCD) on the test problem “a7a” the CPU time and total time (in minutes) to solve knapsack problems (left) and the number of full-iterations (right) for different dimensions of the blocks $n_{i}$. We see that the number of iterations decreases with the increasing dimension of the blocks, while the CPU time increases w.r.t. the scalar case due to the fact that for $n_{i}>1$ the direction $d_{ij}$ cannot be computed in closed form as in the scalar case (i.e., $n_{i}=1$), but requires solving a quadratic knapsack problem (18) whose solution can be computed in $\mathcal{O}(n_{i}+n_{j})$ operations Kiw:07.

2 Chebyshev center of a set of points

where $r$ is the radius and $z_{c}$ is the center of the enclosing ball. It can be immediately seen that the dual formulation of this problem is a particular case of our linearly constrained optimization model (1):

where $Z$ is the matrix containing the given points $z_{i}$ as columns. Once an optimal solution $x^{*}$ for the dual formulation is found, a primal solution can be recovered as follows:

The direction $d_{ij}$ at the current point in the algorithm (RCD) is computed in closed form. For computing the direction in the (CGD) method we need to solve a quadratic knapsack problem that has linear time complexity Kiw:07. The direction at the current point for algorithm (GM) is computed using a linear time simplex projection algorithm introduced in JudRay:08. We compare algorithms (RCD), (CGD) and (GM) for a set of large-scale problem instances generated randomly with a uniform distribution. We recover a suboptimal radius and Chebyshev center using the same set of relations (21) evaluated at the final iteration point $x^{k}$ for all three algorithms.

References