A generic coordinate descent solver for nonsmooth convex optimization

Smooth objectives are a natural starting point for algorithmic innovations. The optimization problem at stake writes

where $f$ is a convex differentiable function with Lipschitz-continuous partial derivatives.

In Table 1, we compare papers that introduced important improvements to coordinate descent methods. We shall in particular stress the seminal paper by Tseng and Yun [60]. It features coordinate gradient steps instead of exact minimization. Indeed a coordinate gradient steps gives similar iteration complexity both in theory and in practice for a much cheaper iteration cost. Moreover, this opened the door for many innovations: blocks of coordinates and the use of proximal operators were developed in the same paper. Another crucial step was made in [39]: Nesterov showed how randomization can help finding finite-time complexity bounds and proposed an accelerated version of coordinate descent. He also proposed to use a non-uniform sampling of coordinates depending on the coordinate-wise Lipschitz constants.

A large literature has been devoted to composite optimization problems with separable non-smooth functions:

where $f$ is a convex differentiable function with Lipschitz-continuous partial derivatives and for all $i$, $g_{i}$ is a convex function whose proximal operator is easy to compute. Indeed, regularized expected risk minimization problems often fit into this framework and this made the success of coordinate descent methods for machine learning applications.

Some papers study $\min_{x\in\mathbb{R}^{p}}f(x)+\sum_{i=1}^{n}g_{i}((Ax)^{(i)})$, where $f$ is strongly convex and apply coordinate descent to a dual problem written as

One of the challenges of these works is to show that even though we are solving the dual problem, one can still recover rates for a sequence minimizing the primal objective.

We present our selection of papers devoted to this type of problems in Table 2.

Non-separable non-smooth objective functions are much more challenging to coordinate descent methods. One wishes to solve

where $f$ is a convex differentiable function with Lipschitz-continuous partial derivatives, $g$ and $h$ are convex functions whose proximal operator are easy to compute and $A$ is a linear operator. Indeed, the linear operator introduces a coupling between the coordinates and a naive approach leads to a method that does not converge to a minimizer [3]. When $h=\iota_{\{b\}}$, the convex indicator function of the set $\{b\}$, we have equality constraints.

We present our selection of papers devoted to this type of problems in Table 3.

The goal here is to find a local minimum to the problem

without any assumption on the convexity of $f$ nor $g_{i}$. The function $f$ should be continuously differentiable and the proximal operator of each function $g_{i}$ should be easily computable. Note that in the non-convex setting, the proximal operator may be set-valued.

The algorithm we implemented is a coordinate descent primal-dual method developed in [15]. Let introduce the notation $F(x)=\frac{1}{2}x^{\top}Qx+\sum_{j=1}^{J}c^{f}_{j}f_{j}(A^{f}_{j}x-b^{f}_{j})$, $G(x)=\sum_{i=1}^{I}c^{g}_{i}g_{i}(D^{g}_{i}x^{(i)}-b^{g}_{i})$, $H(z)=\sum_{l=1}^{L}c^{h}_{l}h_{l}(z^{(l)}-b^{h}_{l})$, $\mathcal{J}(i)=\{j\;:\;A^{h}_{j,i}\neq 0\}$, $\mathcal{I}(j)=\{i\;:\;A^{h}_{j,i}\neq 0\}$, $m_{j}=|\mathcal{I}(j)|$ and $\rho(A)$ the spectral radius of matrix $A$. We shall also denote $\mathcal{J}^{f}(i)=\{j\;:\;A^{f}_{j,i}\neq 0\}$, $\mathcal{J}^{Q}(i)=\{j\;:\;Q_{j,i}\neq 0\}$, $A^{f}\in\mathbb{R}^{\sum_{j}M^{f}_{j}\times N}$ the matrix which stacks the matrices $(A^{f}_{j})_{1\leq j\leq J}$ and $A^{h}$ the matrix which stacks the matrices $(A^{h}_{l})_{1\leq l\leq L}$. The algorithm writes then as Algorithm 1.

We will denote $U_{1},\ldots,U_{I}$ the columns of the identity matrix corresponding to the blocks of $x=(x^{(1)},\ldots,x^{(I)})$, so that $U_{i}x^{(i)}\in\mathbb{R}^{N}$ and $V_{1},\ldots,V_{J}$ the columns of the identity matrix corresponding to the blocks of $A^{f}x-b^{f}=(A^{f}_{1}x-b^{f}_{1},\ldots,A^{f}_{J}x-b^{f}_{J})$.

For simplicity of implementation, we are assuming that $G$ is separable and the blocks of variable will follow the block structure of $G$. This implies in particular that at each iteration, only $\nabla_{i}F(x_{k})$ needs to be computed. This partial derivative needs to be calculated efficiently because it needs to be performed at each iteration of the algorithm. We now describe the efficient residual update method, which is classically used in coordinate descent implementations [39].

Denote $r^{f,x}_{k}=A^{f}x_{k}-b^{f}$ and $r^{Q,x}_{k}=Qx_{k}$. By the chain rule, we have

If $r^{f,x}_{k}$ and $r^{Q,x}_{k}$ are pre-computed, only $O(|\mathcal{J}^{f}(i)|+|\mathcal{J}^{Q}(i)|)$ operations are needed.

For an efficient implementation, we will update the residuals $r^{f,x}_{k}$ as follows, using the fact that only the coordinate block $i_{k+1}$ is updated:

Hence, updating $r^{f,x}_{k+1}$ also requires only $O(|\mathcal{J}^{f}(i_{k+1})|)$ iterations.

Similarly, updating the residuals $r^{Q,x}_{k}$, $r^{h,x}_{k}=A^{h}x_{k}-b^{h}$, $w_{k}$ and $z_{k}$ can be done in $O(|\mathcal{J}^{Q}(i_{k+1})|)$ and $O(|\mathcal{J}(i_{k+1})|)$ operations.

Although this technique is well known, it is not trivial how to write it in a generic fashion, since residual updates are needed at each iteration and should be written in a compiled language. We coded the residual update using abstract atom functions in order to achieve this goal.

Another major step in the method is the computation of the $i^{\text{th}}$ coordinate of $\operatorname{prox}_{\tau,G}(x)$ for a given $x\in\mathbb{R}^{N}$.

As $D^{g}$ is assumed to be diagonal, $G$ is separable. Hence, by the change of variable $\bar{z}=D^{g}_{i}\bar{x}-b^{g}_{i}$,

where we used the abuse of notation that $D_{i}^{g}$ is either the scaled identity matrix or any of its diagonal elements. This derivation shows that to compute $(\operatorname{prox}_{\tau,G}(x))_{i}$ we only need linear algebra and the proximal operator of the atom function $g_{i}$.

We can similarly compute $\operatorname{prox}H$. To compute $\operatorname{prox}_{\sigma,H^{\star}}$, we use Moreau’s formula:

Algorithm 1 maintains duplicated dual variables $\boldsymbol{y}_{k}\in\mathbb{R}^{\mathrm{nnz}(A^{h})}$ as well as averaged dual variables $z_{k}\in\mathbb{R}^{M^{h}}$ where $M^{h}=\sum_{l=1}^{L}M^{h}_{l}$ and $A^{h}_{l,i}$ is of size $M_{l}^{h}\times N_{i}$. The sets $\mathcal{J}(i)$ for all $i$ are given by the sparse column format representation of $A^{h}$. Yet, for all $i$, we need to construct the set of indices of $\boldsymbol{y}_{k+1}$ that need to be updated. This is the table dual_vars_to_update in the code. Moreover, as $H$ is not separable in general, in order to compute $\bar{y}_{k+1}^{j}$, for $j\in\mathcal{J}(i_{k+1})$, we need to determine the set of dual indices $j^{\prime}$ that belong to the same block as $j$ with respect to the block decomposition of $H$. This is the purpose of the tables inv_blocks_h and blocks_h.

The procedure allows us to only compute the entries of $\bar{y}_{k+1}$ that are required for the update of $\boldsymbol{y}_{k}$.

We added the following feature for an improved efficiency. Under the argument sampling=’kink_half’, the algorithms periodically detects the set of blocks $I_{\mathrm{kink}}$ such that $i\in I_{\mathrm{kink}}$ if $x^{(i)}$ is at a kink of $g_{i}$. Then, block $i$ is selected with probability law

The rationale for this probability law is that blocks at kinks are likely to incur no move when we try to update them. We thus put more computational power for non-kinks. On the other hand, we still keep an update probability weight of at least $\frac{1}{2n}$ for each block, so even in unfavourable cases, we should not observe too much degradation in the performance as compared to the uniform law.

We also coded accelerated coordinate descent [18], as well as its restarted [16] and primal-dual [1] variants. The algorithm is given in Algorithm 2. As before, $\bar{y}_{k+1}$ and $\bar{x}_{k+1}$ should not be computed: only the relevant coordinates should be computed.

The accelerated algorithms improve the worst case guarantee as explained in Table 5:

However, accelerated algorithms do not take profit of regularity properties of the objective like strong convexity. Hence, they are not guaranteed to be faster, even though restart may help.

The code includes the Gap Safe screening method presented in [36]. Note that the method has been studied only for the case where $h=0$. Given a non-differentiability point $\hat{x}^{(i)}$ of the function $g_{i}$ where the subdifferential $\partial g_{i}(\hat{x}^{(i)})$ has a non-empty interior, a test is derived to check whether $\hat{x}^{(i)}$ is the $i^{\text{th}}$ variable of an optimal solution. If this is the case, one can set $x^{(i)}=\hat{x}^{(i)}$ and stop updating this variable. This may lead to a huge speed up in some cases. As the test relies on the computation of the duality gap, which has a nonnegligible cost, it is only performed from time to time.

In order to state Gap Safe screening in a general setting, we need the concept of polar of a support function. Let $C$ be a convex set. The support function of $C$ is $\sigma_{C}$ defined by

The polar to $\sigma_{C}$ is

In particular, if $\bar{x}\in C$, then $\sigma_{C}^{\circ}(\bar{x})\leq 1$. Denote $f(z)=\sum_{j=1}^{J}c^{f}_{j}f_{j}(z^{(j)}-b^{f}_{j})$, so that $(\nabla f(A^{f}x))_{j}=c^{f}_{j}\nabla f_{j}(A^{f}_{j}x-b_{j}^{f})$, $G_{i}(x)=c^{g}_{i}g_{i}(D_{i}^{g}x^{(i)}-b_{i}^{g})$, $x_{\star}$ a solution to the optimization problem (1) and suppose we have a set $\mathcal{R}$ such that $(\nabla f(A^{f}x_{\star}-b^{f}),Qx_{\star})\in\mathcal{R}$. Gap Safe screening states that

Denote $z=\nabla f(A^{f}x)$ and $w=Qx$¹¹1We reuse the notation $w$ and $z$ here for the purpose of this section.. We choose $\mathcal{R}$ as a sphere centered at

and with radius

where

and

Note that as $\bar{\omega}$ is a rescaled version of $w=Qx$, we do not need to know $Q^{\dagger}$ in order to compute $\bar{\omega}^{\top}Q^{\dagger}\bar{\omega}$. It is proved in [36] that this set $\mathcal{R}$ contains $(\nabla f(A^{f}x_{\star}-b^{f}),Qx_{\star})$ for any optimal solutions $x_{\star}$ to the primal problem. In the case where $G_{i}$ is a norm and $x^{(i)}=0$, these expressions simplify since the $\sigma^{\circ}_{\partial G_{i}(0)}=\sigma^{\circ}_{\operatorname{dom}G_{i}^{*}}$ is nothing else than the dual norm associated to $G_{i}$.

For the estimation of $\max_{(\zeta,\omega)\in\mathcal{R}}\sigma^{\circ}_{\partial G_{i}(\hat{x}^{(i)})}((A^{f})_{i}^{\top}\zeta+\omega^{(i)})$, we use the fact that the polar of a support function is sublinear and positively homogeneous. Indeed, we have