A generic coordinate descent solver for nonsmooth convex optimization

TL;DR

This paper introduces a versatile coordinate descent solver capable of efficiently handling nonsmooth convex problems with structure, including those with linear constraints, by leveraging residual updates and a flexible modeling language.

Contribution

The paper presents a generic coordinate descent algorithm with an efficient implementation and a Python-based modeling language, enabling broad applicability to various nonsmooth convex optimization problems.

Findings

Effective for Lasso and sparse logistic regression

Handles linear and quadratic programming problems

Automatically manages dual variable duplication

Abstract

We present a generic coordinate descent solver for the minimization of a nonsmooth convex objective with structure. The method can deal in particular with problems with linear constraints. The implementation makes use of efficient residual updates and automatically determines which dual variables should be duplicated. A list of basic functional atoms is pre-compiled for efficiency and a modelling language in Python allows the user to combine them at run time. So, the algorithm can be used to solve a large variety of problems including Lasso, sparse multinomial logistic regression, linear and quadratic programs.

Tables

Table 1. Table 1: Selected papers for the minimization of smooth functions. Rate : we check whether the paper proves convergence ( × \times ), an asymptotic rate (asymp) or a finite time iteration complexity (✓). Rand : deterministic (N) or randomized (Y) selection of coordinates. Grad : exact minimization ( × \times ) or use of partial derivatives (✓). Blck : the paper considers 1D coordinates ( × \times ) or blocks of coordinates (✓). Par : Algorithm designed for parallel computing (✓). Acc : no momentum ( × \times ), accelerated but not efficient in practice (not eff), accelerated algorithm (✓)

Paper	​​Rate​​	​​Rand​​	​​Grad​​	​​Blck​​	​Par​	​Acc ​	Notable feature
Seidel ’74 [49]	$\times$	N	$\times$	$\times$	$\times$	$\times$	quadratic
Warga ’63 [62]	$\times$	N	$\times$	$\times$	$\times$	$\times$	strictly convex
Luo & Tseng ’92 [34]	​​asymp​​	N	$\times$	$\times$	$\times$	$\times$	rate for weakly convex
Leventhal & Lewis ’08 [28]	✓	Y	$\times$	$\times$	$\times$	$\times$	quadratic $f$
Tseng & Yun ’09 [60]	​​asymp​​	N	✓	✓	$\times$	$\times$	​​line search, proximal operator​​
Nesterov ’12 [39]	✓	Y	✓	$\times$	$\times$	​​not eff​​	1st acc & 1st non-uniform
Beck & Tetruashvili ’13 [4]	✓	N	✓	✓	$\times$	​​not eff​​	​​finite time analysis cyclic CD​​
Lin & Xiao ’13 [33]	✓	Y	✓	✓	$\times$	​​not eff​​	improvements on [39, 47]
Lee & Sidford ’13 [27]	✓	Y	✓	$\times$	$\times$	✓	1st efficient accelerated
Liu et al ’13 [31]	✓	Y	✓	$\times$	✓	$\times$	1st asynchronous
Glasmachers ​&​ Dogan ​’13​ [22]​​​	$\times$	Y	✓	$\times$	$\times$	$\times$	heuristic sampling
Richtárik & Takáč ’16 [46]	✓	Y	✓	$\times$	✓	$\times$	1st arbitrary sampling
Allen-Zhu et al ’16 [2]	✓	Y	✓	✓	$\times$	✓	non-uniform sampling
Sun et al ’17 [55]	✓	Y	✓	✓	✓	$\times$	better asynchrony than [31]

Table 2. Table 2: An overview of selected papers proposing and analyzing the iteration complexity of coordinate descent methods for separable non-smooth objectives. Prx : uses a proximal operator to deal with the non-smooth part of the objective. Par : updates several blocks of coordinates in parallel. Acc : uses momentum to obtain an improved rate of convergence. Dual : solves a dual problem but still proves rates in the primal (only relevant for weakly convex duals).

Paper	Prx	Par	Acc	Dual	Notable feature
Tseng & Yun ’09 [60]	✓	$\times$	$\times$	$\times$	1st prox, line search, deterministic
S-Shwartz & Tewari ’09 [50]	${\ell }_1$	$\times$	$\times$	$\times$	1st ${\ell }_1$-regularized w finite time bound
Bradley et al ’11 [7]	${\ell }_1$	✓	$\times$	$\times$	${\ell }_1$-regularized parallel
Richtárik & Takáč ’14 [47]	✓	$\times$	$\times$	$\times$	1st proximal with finite time bound
S-Shwartz & Zhang ’13 [51]	✓	$\times$	$\times$	✓	1st dual
Richtárik & Takáč ’15 [48]	✓	✓	$\times$	$\times$	1st general parallel
Takáč et al ’13 [56]	✓	✓	$\times$	✓	1st dual & parallel
S-Shwartz & Zhang ’14 [53]	✓	$\times$	✓	✓	acceleration in the primal
Yun ’14 [68]	✓	$\times$	$\times$	$\times$	analysis of cyclic CD
Fercoq & Richtárik ’15 [18]	✓	✓	✓	$\times$	1st proximal and accelerated
Lin, Lu & Xiao ’14 [30]	✓	$\times$	✓	$\times$	prox & accelerated on strong conv.
Richtárik & Takáč ’16 [45]	✓	✓	$\times$	$\times$	1st distributed
Fercoq et al ’14 [17]	✓	✓	✓	$\times$	distributed computation
Lu & Xiao ’15 [33]	✓	$\times$	$\times$	$\times$	improved complexity over [47, 39]
Li & Lin ’18 [29]	✓	$\times$	✓	✓	acceleration in the dual
Fercoq & Qu ’18 [16]	✓	$\times$	✓	$\times$	restart for obj with error bound

Table 3. Table 3: Selected papers for the minimization of non-smooth non-separable functions. Rate : does the paper prove rates? Const : can the method solve problems with linear equality constraints? PD-CD : the method is purely primal (P), the method updates some primal variables but all the Lagrange multipliers or some dual variables but all the primal variables ( × \times ), the method updates some primal and some dual variables at each iteration (✓).

Paper	Rate	Const	PD-CD	Notable feature
Platt ’99 [43]	$\times$	✓	P	for SVM
Tseng & Yun ’09 [61]	✓	✓	P	adapts Gauss-Southwell rule
Tao et al ’12 [57]	✓	$\times$	P	uses averages of subgradients
Necoara et al ’12 [38]	✓	✓	P	2-coordinate descent
Nesterov ’12 [40]	✓	$\times$	P	uses subgradients
Necoara & Clipici ’13 [37]	✓	✓	P	coupled constraints
Combettes & Pesquet ’14 [11]	$\times$	✓	✓	1st PD-CD, short step sizes
Bianchi et al ’14 [5]	$\times$	✓	✓	distributed optimization
Hong et al ’14 [24]	$\times$	✓	$\times$	updates all dual variables
Fercoq & Richtárik ’17 [19]	✓	$\times$	P	uses smoothing
Alacaoglu et al ’17 [1]	✓	✓	✓	1st PD-CD w rate for constraints
Xu & Zhang ’18 [66]	✓	✓	$\times$	better rate than [21]
Chambolle et al ’18 [9]	✓	✓	$\times$	updates all primal variables
Fercoq & Bianchi ’19 [15]	✓	✓	✓	1st PD-CD w long step sizes
Gao et al ’19 [21]	✓	✓	$\times$	1st primal-dual w rate for constraints
Latafat et al ’19 [26]	✓	✓	✓	linear conv w growth condition

Table 4. Table 4: Selected papers for the local minimization of nonconvex problems. Conv : study only limit points ( × \times ) or convergence of the sequence proved (✓). Smth : can deal with non-convex smooth functions. Nsmth : can deal with non-convex and non-smooth functions.

Paper	Conv	Smth	Nsmth	Notable feature
Grippo & Sciandrone ’00 [23]​​	$\times$	✓	$\times$	2-block or coordinate-wise quasiconvex
Tseng & Yun ’09 [60]	✓	✓	$\times$	convergence under error bound
Hsieh & Dhillon ’11 [25]	$\times$	✓	$\times$	non-negative matrix factorization
Breheny & Huang ’11 [8]	$\times$	$\times$	✓	regularized least squares
Mazumder et al ’11 [35]	$\times$	$\times$	✓	regularized least squares
Razaviyayn et al ’13 [44]	$\times$	✓	✓	requires uniqueness of the prox
Xu & Yin ’13 [64]	✓	✓	$\times$	multiconvex
Lu & Xiao ’13 [32]	$\times$	✓	✓	random sampling
Patrascu & Necoara ’15 [41]	✓	✓	$\times$	randomized + 1 linear constraint
Xu & Yin ’17 [65]	✓	✓	✓	convergence under KL

Table 5. Table 5: Convergence speed of the algorithms implemented

	$h=0$	$h\ne 0$
PD-CD Alg. 1	$O\left(1/k\right)$	$O\left(1/\sqrt{k}\right)$
APPROX / SMART-CD Alg. 2	$O\left(1/k^2\right)$	$O\left(1/k\right)$

Table 6. Table 6: Comparison of our code with a pure Python code and reference implementations for performing 100 N 𝑁 N coordinate descent iterations for the Lasso problem on the Leukemia dataset with regularization parameter λ = 0.1 ​ ∥ ( A f ) ⊤ ​ b f ∥ ∞ 𝜆 0.1 subscript delimited-∥∥ superscript superscript 𝐴 𝑓 top superscript 𝑏 𝑓 \lambda=0.1\lVert(A^{f})^{\top}b^{f}\rVert_{\infty} , and for 10 N 𝑁 N coordinate descent iterations for the dual SVM problem on the RCV1 dataset with penalty parameter C = 10 𝐶 10 C=10 .

Lasso
Pure Python	308.76s
cd_solver	0.43s
Scikit learn Lasso	0.11s

Table 7. Table 7: Time to reach precision ϵ italic-ϵ \epsilon . scr = with gap safe screening. For SVM, we scale the data so that in each column of the data matrix A 𝐴 A , A i , j ∈ [ − 1 , 1 ] subscript 𝐴 𝑖 𝑗 1 1 A_{i,j}\in[-1,1] .

Problem	Dataset	Alternative solver 1	Alternative solver 2	Algorithm 1	Algorithm 2
Lasso ($\lambda =0.1{\parallel X^{\top }y\parallel }_{\infty }$)	triazines ($J$=186, $I$=60)	Scikit-learn: 0.005s	OSQP: 0.033s	0.107s; scr: 0.101s	0.066s; scr: 0.079s
$ϵ={10}^{-6}$	scm1d ($J$=9,803, $I$=280)	Scikit-learn: 24.40s	OSQP: 33.21s	91.97s; scr: 8.73s	24.13; scr: 3.21s
	news20.binary ($J$=19,996, $I$=1,355,191)	Scikit-learn: 64.6s	OSQP: $>$2,000s	267.3s; scr: 114.4s	169.8s; scr: 130.7s
TV-regularized regression	fMRI [58, 13] $\left({\alpha }_1={\alpha }_2={5.10}^{-3}\right)$	inexact FISTA: 24,341s	LBFGS+homotopy: 6,893s	$>$25,000s	2,734s
$ϵ=1$	($J$=768, $I$=65,280, $L$=195,840)
Dual SVM with intercept	ionosphere ($J$=14, $I$=351)	libsvm: 0.04s	OSQP: 0.12s	3.23s	0.42s
$ϵ={10}^{-3}$	leukemia ($J$=7,129, $I$=72)	libsvm: 0.1s	OSQP: 3.2s	25.5s	0.8s
	madelon ($J$=500, $I$=2,600)	libsvm: 50s	OSQP: 37s	3842s	170s
	gisette ($J$=5,000, $I$=6,000)	libsvm: 70s	OSQP: 936s	170s	901s
	rcv1 ($J$=47,236, $I$=20,242)	libsvm: 195s	OSQP: Memory error	$>$5000s	63s

Equations

$$\underset{x\in {\mathbb{R}}^N}{\min }\sum _{j=1}^{J}f(A_{j}x-b_j)+\sum _{i=1}^{I}g(x^{(i)})$$

$$\underset{x\in {\mathbb{R}}^N}{\min }\frac{1}{2}{x}^{\top }Qx+\sum _{j=1}^J{c}_j^f{f}_j(A_j^{f}x-b_j^f)+\sum _{i=1}^I{c}_i^g{g}_i(D_i^g{x}^{(i)}-b_i^g)+\sum _{l=1}^L{c}_l^h{h}_l(A_l^{h}x-b_l^h).$$

$$\underset{x\in {\mathbb{R}}^N}{\min }f(x)$$

$$\underset{x\in {\mathbb{R}}^N}{\min }f(x)+\sum _{i=1}^n{g}_i(x^{(i)})$$

$$\underset{y\in {\mathbb{R}}^N}{\min }{f}^{\ast }(-A^{\top }y)+\sum _{i=1}^n{g}_i^{\ast }(y^{(i)}).$$

$$\underset{x\in {\mathbb{R}}^N}{\min }f(x)+g(x)+h(Ax)$$

$$\underset{x\in {\mathbb{R}}^N}{loc-min}f(x)+\sum _{i=1}^n{g}_i(x^{(i)})$$

$${\tau }_i<\frac{1}{{\beta }_i+\rho \left({\sum }_{j\in \mathcal{J}(i)}{m}_j{\sigma }_j(A^h{)}_{j,i}^{\top }{A}_{j,i}^h\right)}.$$

$${\overline{y}}_{k+1}$$

$${\overline{x}}_{k+1}$$

$$x_{k+1}^{(i)}={\overline{x}}_{k+1}^{(i)}$$

$${\mathbf{y}}_{k+1}^{(j)}(i)={\overline{y}}_{k+1}^{(j)}$$

$$w_{k+1}^{(i)}=w_k^{(i)}+\sum _{j\in J(i)}(A^h{)}_{j,i}^{\top }({\mathbf{y}}_{k+1}^{(j)}(i)-{\mathbf{y}}_k^{(j)}(i))$$

$$z_{k+1}^{(j)}=z_k^{(j)}+\frac{1}{m_j}({\mathbf{y}}_{k+1}^{(j)}(i)-{\mathbf{y}}_k^{(j)}(i)).$$

$${\nabla }_{i}F(x_k)=\sum _{j=1}^J{c}_j^f(A^f{)}_{j,i}^{\top }\nabla f_j(A_j^f{x}_k-b_j^f)+Q{x}_k=\sum _{j\in {\mathcal{J}}^f(i)}{c}_j^f(A^f{)}_{j,i}^{\top }\nabla f_j((r_k^{f,x}{)}_j)+\sum _{j\in {\mathcal{J}}^Q(i)}{U}_i{r}_k^{Q,x}$$

$$r^{f,x}_{k+1}=A^{f}x_{k+1}-b^{f}=A^{f}\big{(}x_{k}+U_{i_{k+1}}(x_{k+1}^{(i_{k+1})}-x_{k}^{(i_{k+1})})\big{)}-b^{f}=r^{f,x}_{k}+A^{f}U_{i_{k+1}}(x_{k+1}^{(i_{k+1})}-x_{k}^{(i_{k+1})})\\ =r^{f,x}_{k}+\sum_{j\in\mathcal{J}^{f}(i_{k+1})}V_{j}A^{f}_{j,i_{k+1}}(x_{k+1}^{(i_{k+1})}-x_{k}^{(i_{k+1})})$$

$$({\mathrm{prox}}_{\tau ,G}(x){)}_i$$

$$\displaystyle=(D_{i}^{g})^{-1}\Big{(}b^{g}_{i}+\arg\min_{\bar{z}\in\mathbb{R}^{N_{i}}}c_{i}^{g}g_{i}(\bar{z})+\frac{1}{2\tau_{i}}\lVert(D_{i}^{g})^{-1}(b^{g}_{i}+\bar{z})-x^{(i)}\rVert^{2}\Big{)}$$

$$\displaystyle=(D_{i}^{g})^{-1}\Big{(}b^{g}_{i}+\arg\min_{\bar{z}\in\mathbb{R}^{N_{i}}}g_{i}(\bar{z})+\frac{1}{2c_{i}^{g}(D_{i}^{g})^{2}\tau_{i}}\lVert\bar{z}-(D^{g}_{i}x^{(i)}-b^{g}_{i})\rVert^{2}\Big{)}$$

$$\displaystyle=(D_{i}^{g})^{-1}\Big{(}b^{g}_{i}+\operatorname{prox}_{c_{i}^{g}(D_{i}^{g})^{2}\tau_{i}g_{i}}(D_{i}^{g}x^{(i)}-b_{i}^{g})\Big{)}$$

$${\mathrm{prox}}_{\sigma ,H^{\star }}(z)=z-D(\sigma ){\mathrm{prox}}_{{\sigma }^{-1},H}(D(\sigma {)}^{-1}z)$$

$$\mathbb{P}(i_{k+1}=i)=\{\begin{array}{ll}\frac{1}{n}& \text{ if }|I_{\mathrm{kink}}|=n\\ \frac{1}{2n}& \text{ if }|I_{\mathrm{kink}}|<n\text{ and }i\in I_{\mathrm{kink}}\\ \frac{1}{2n}+\frac{1}{2(n-|I_{\mathrm{kink}}|)}& \text{ if }|I_{\mathrm{kink}}|<n\text{ and }i\notin I_{\mathrm{kink}}\end{array}$$

$${\overline{y}}_{k+1}$$

$${\overline{x}}_{k+1}$$

$${\tilde{x}}_{k+1}^{(i)}={\overline{x}}_{k+1}^{(i)}$$

$${\hat{x}}_{k+1}^{(i)}={\hat{x}}_k^{(i)}-\frac{1-{\theta }_k/{\theta }_0}{c_k}({\tilde{x}}_{k+1}^{(i)}-{\tilde{x}}_k^{(i)})$$

$$\{\begin{array}{ll}{\theta }^3+{\theta }^2+{\theta }_k^2\theta -{\theta }_k^2=0& \text{ if }h\ne 0\\ {\theta }^2+{\theta }_k^2\theta -{\theta }_k^2=0& \text{ if }h=0\end{array}$$

$${\sigma }_C(x)=\underset{y\in C}{\sup }⟨y,x⟩.$$

$${\sigma }_C^{\circ }(\bar{x})=\underset{x:{\sigma }_C(x)\le 1}{\sup }⟨\bar{x},x⟩.$$

$$\underset{(\zeta ,\omega )\in \mathcal{R}}{\max }{\sigma }_{\partial G_i({\hat{x}}^{(i)})}^{\circ }((A^f{)}_i^{\top }\zeta +{\omega }^{(i)})<1\Rightarrow x_{\star }^{(i)}={\hat{x}}^{(i)}.$$

$$(\bar{\zeta},\bar{\omega})=\frac{(z,w)}{\max(1,\max_{1\leq i\leq I}\big{(}\sigma_{\operatorname{dom}G_{i}^{*}}^{\circ}((A_{i}^{f})^{\top}z+w^{(i)})\big{)})}$$

$$r=\sqrt{\frac{2\mathrm{Gap}(x,\bar{\zeta },\bar{\omega })}{L_{f,Q}}}$$

$$L_{f,Q}=\max (\underset{1\le j\le J}{\max }L(\nabla f_j),\underset{1\le i\le I}{\max }\rho (Q_{i,i}))$$

$$\mathrm{Gap}(x,\bar{\zeta },\bar{\omega })=\frac{1}{2}{x}^{\top }Qx+f(A^{f}x)+G(x)+G^{\ast }(-(A^f{)}^{\top }\bar{\zeta }-\bar{\omega })+\frac{1}{2}{\bar{\omega }}^{\top }{Q}^{†}\bar{\omega }+f^{\ast }(\bar{\zeta }).$$

$${\sigma }_{\partial G_i({\hat{x}}^{(i)})}^{\circ }$$

Full text

1 Introduction

Coordinate descent methods decompose a large optimization problem into a sequence of one-dimensional optimization problems. The algorithm was first described for the minimization of quadratic functions by Gauss and Seidel in [49]. Coordinate descent methods have become unavoidable in machine learning because they are very efficient for key problems, namely Lasso [20], logistic regression [67] and support vector machines [43, 52]. Moreover, the decomposition into small subproblems means that only a small part of the data is processed at each iteration and this makes coordinate descent easily scalable to high dimensions.

One of the main ingredients of an efficient coordinate descent solver is its ability to compute efficiently partial derivatives of the objective function [39]. In the case of least squares for instance, this involves the definition of a vector of residuals that will be updated during the run of the algorithm. As this operation needs to be performed at each iteration, and millions of iterations are usually needed, the residual update and directional derivative computation must be coded in a compiled programming language.

Many coordinate descent solvers have been written in order to solve a large variety of problems. However, most of the existing solvers can only solve problems of the type

[TABLE]

where $x^{(i)}\in\mathbb{R}^{N_{i}}$ is the $i$th block of $x$, $\sum_{i=1}^{I}N_{i}=N$, $A_{j}\in\mathbb{R}^{M_{j}\times N}$ is a matrix and $b_{j}\in\mathbb{R}^{M_{j}}$ is a vector, and where $f$ is a convex differentiable function and $g$ is a convex lower-semicontinuous function whose proximal operator is easy to compute (a.k.a. a proximal-friendly convex function). Each piece of code usually covers only one type of function [14, 42]. Moreover, even when the user has a choice of objective function, the same function is used for every block [6].

In this work, we propose a generic coordinate descent method for the resolution of the convex optimization problem

[TABLE]

We shall call $f_{j}$, $g_{i}$ and $h_{l}$ atom functions. Each of them may be different. We will assume that $f_{j}$’s are differentiable and convex, $g_{i}$’s and $h_{l}$’s are proximal-friendly convex functions. As before $A^{f}_{j}\in\mathbb{R}^{M^{f}_{j}\times N}$ and $A^{h}_{l}\in\mathbb{R}^{M^{h}_{l}\times N}$ are matrices, $D_{i}^{g}$ is a multiple of the identity matrix of size $N_{i}$, $b^{f}_{j}\in\mathbb{R}^{M^{f}_{j}}$, $b^{g}_{i}\in\mathbb{R}^{N_{i}}$ and $b^{h}_{l}\in\mathbb{R}^{M^{h}_{l}}$ are vectors, $c^{f}_{j}$, $c^{g}_{i}$ and $c^{h}_{l}$ are positive real numbers, $Q$ is a $N\times N$ positive semi-definite matrix.

The algorithm we implemented is described in [15] and can be downloaded on https://bitbucket.org/ofercoq/cd_solver. The present paper focuses on important implementation details about residual updates and dual variable duplication. The novelty of our code is that it allows a generic treatment of these algorithmic steps and includes a modelling interface in Python for the definition of the optimization problem. Note that unlike most coordinate descent implementations, it can deal with nonseparable nonsmooth objectives and linear constraints.

2 Literature review on coordinate descent methods

A thorough review on coordinate descent is beyond the scope of this paper. We shall refer the interested reader to the review papers [63] and [54]. Instead, for selected papers dealing with smooth functions, separable non-smooth functions or non-separable non-smooth function, we list their main features. We also quickly review what has been done for non-convex functions. We sort papers in order of publication except when there are an explicit dependency between a paper and a follow-up.

2.1 Smooth functions

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

[TABLE]

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.

2.2 Separable non-smooth functions

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

[TABLE]

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

[TABLE]

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.

2.3 Non-separable non-smooth functions

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

[TABLE]

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.

2.4 Non-convex functions

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

[TABLE]

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.

3 Description of the Algorithm

3.1 General scheme

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})$.

3.2 Computation of partial derivatives

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

[TABLE]

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:

[TABLE]

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.

3.3 Computation of proximal operators using atom functions

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}$,

[TABLE]

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:

[TABLE]

3.4 Duplication of dual variables

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}$.

4 Code structure

The code is organized in nine files. The main file is cd_solver.pyx. It contains the Python callable and the data structure for the problem definition. The other files are atoms.pyx/pxd, algorithm.pyx/pxd, helpers.pyx/pxd and screening.pyx/pxd. They contain the definition of the atom functions, the algorithms and the functions for computing the objective value. In Figure 1, we show for each subfunction, in which function it is used. The user needs to call the Python class Problem and the Python function coordinate_descent. Atom functions can be added by the user without modifying the main algorithm.

All tables are defined using Numpy’s array constructor in the coordinate_descent function. The main loop of coordinate descent and the atom functions are pre-compiled for efficiency.

5 Atom functions

The code allows us to define atom functions independently of the coordinate descent algorithm. As an example, we provide in Figure 2 the code for the square function atom.

As inputs, it gets x (an array of numbers which is the point where the operation takes place), buff (the buffer for vectorial outputs), nb_coord (is the size of x), mode, prox_param and prox_param2 (numbers which are needed when computing the proximal operator). The input mode can be:

• •

GRAD in order to compute the gradient.

• •

PROX to compute the proximal operator.

• •

PROX_CONJ uses Moreau’s formula to compute the proximal operator of the conjugate function.

• •

LIPSCHITZ to return the Lipschitz constant of the gradient.

• •

VAL_CONJ to return the value of the conjugate function. As this mode is used only by compute_smoothed_gap for printing purposes, its implementation is optional and can be approximated using the helper function val_conj_not_implemented. Indeed, for a small $\epsilon>0$, $h^{*}(y)=\sup_{z}\langle z,y\rangle-h(z)\approx\sup_{z}\langle z,y\rangle-h(z)-\frac{\epsilon}{2}\lVert z\rVert^{2}=\langle p,y\rangle-h(p)-\frac{\epsilon}{2}\lVert p\rVert^{2}$, where $p=\operatorname{prox}_{h/\epsilon}(y/\epsilon)$.

• •

VAL to return the value of the function.

Some functions naturally require multi-dimensional inputs, like $\lVert\cdot\rVert_{2}$ or the log-sum-exp function. For consistency, we define all the atoms with multi-dimensional inputs: for an atom function $f_{0}:\mathbb{R}\to\mathbb{R}$, we extend it to an atom function $f:\mathbb{R}^{N_{i}}\to\mathbb{R}$ by $f(x)=\sum_{l=1}^{N_{i}}f_{0}(x_{l})$.

For efficiency purposes, we are bypassing the square atom function when computing a gradient and implemented it directly in the algorithm.

6 Modelling language

In order to use the code in all its generality, we defined a modelling language that can be used to define the optimization problem we want to solve (1).

The user defines a problem using the class Problem. Its arguments can be:

• •

N the number of variables, blocks the blocks of coordinates coded in the same fashion as the indptr index of sparse matrices (default [0, 1, …, N]), x_init the initial primal point (default 0) and y_init the initial duplicated dual variable (default 0)

• •

Lists of strings f, g and h that code for the atom functions used. The function string_to_func is responsible for linking the atom function that corresponds to the string. Our convention is that the string code is exactly the name of the function in atoms.pyx. The size of the input of each atom function is defined in blocks_f, blocks and blocks_h. The function strings f, g or h may be absent, which means that the function does not appear in the problem to solve.

• •

Arrays and matrices cf, Af, bf, cg, Dg, bg, ch, Ah, bh, Q. The class initiator transforms matrices into the sparse column format and checks whether Dg is diagonal.

For instance, in order to solve the Lasso problem, $\min_{x}\frac{1}{2}\lVert Ax-b\rVert^{2}_{2}+\lambda\lVert x\rVert_{1}$, one can type

pb_lasso = cd_solver.Problem(N=A.shape[1], f=["square"]*A.shape[0], Af=A,

bf=b, cf=[0.5]*A.shape[0], g=["abs"]*A.shape[1], cg=[lambda]*A.shape[1])

cd_solver.coordinate_descent(pb_lasso)

7 Extensions

7.1 Non-uniform probabilities

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

[TABLE]

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.

7.2 Acceleration

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.

7.3 Variable screening

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

[TABLE]

The polar to $\sigma_{C}$ is

[TABLE]

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

[TABLE]

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

[TABLE]

and with radius

[TABLE]

where

[TABLE]

and

[TABLE]

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

[TABLE]

Here $\hat{x}_{g_{i}}=D^{g}_{i}\hat{x}^{(i)}-b^{g}_{i}$ is a point where $\partial g_{i}(\hat{x}_{g_{i}})$ has a nonempty interior. Some care should be taken when $\partial G_{i}(\hat{x}^{(i)})$ is unbounded, so that we first check whether $(A^{f})_{i}^{\top}\zeta+\omega^{(i)}\in\operatorname{dom}\sigma^{\circ}_{\partial g_{i}(\hat{x}_{g_{i}})}$ for all $(\zeta,\omega)\in\mathcal{R}$.

Here also, the novelty lies in the genericity of the implementation.

8 Numerical validation

8.1 Performance

In order to evaluate the performance of the implementation, we compare our implementation with a pure Python coordinate descent solver and code written for specific problems: Scikit learn’s Lasso solver and Liblinear’s SVM solver. We run the code on an Intel Xeon CPU at 3.07GHz.

We can see on Table 6 that our code is hundreds of times faster than the pure Python code. This is due to the compiled nature of our code, that does not suffer from the huge number of iterations required by coordinate descent. On the other hand, our code is about 4 times slower than state-of-the-art coordinate descent implementations designed for a specific problem. We can see it in both examples we chose. This overhead is the price of genericity.

We believe that, except for critical applications like Lasso or SVM, a 4 times speed-up does not justify writing a new code from scratch, since a separate piece of code for each problem makes it difficult to maintain and to improve with future algorithmic advances.

8.2 Genericity

We tested our algorithm on the following problems:

• •

Lasso problem

[TABLE]

• •

Binomial logistic regression

[TABLE]

where $b_{i}\in\{-1,1\}$ for all $i$.

• •

Sparse binomial logistic regression

[TABLE]

• •

Dual SVM without intercept

[TABLE]

where $\iota_{[0,1]^{n}}$ is the convex indicator function of the set $[0,1]^{n}$ and encodes the constraint $x\in[0,1]^{n}$.

• •

Dual SVM with intercept

[TABLE]

• •

Linearly constrained quadratic program

[TABLE]

• •

Linear program

[TABLE]

• •

TV+$\ell_{1}$-regularized regression

[TABLE]

where $D$ is the discrete gradient operator and $\lVert y\rVert_{2,1}=\sum_{i,j,k}\sqrt{\sum_{l=1}^{3}y_{i,j,k,l}^{2}}$.

• •

Sparse multinomial logistic regression

[TABLE]

where $b_{i,j}\in\{0,1\}$ for all $i,j$.

This list demonstrates that the method is able to deal with differentiable functions, separable or nonseparable nondifferentiable functions, as well as use several types of atom function in a single problem.

8.3 Benchmarking

In Table 7, we compare the performance of Algorithm 1 with and without screening, Algorithm 2 with and without screening as well as 2 alternative solvers for 3 problems exhibiting various situations:

• •

Lasso: the nonsmooth function in the Lasso problem is separable;

• •

the TV-regularized regression problem has a nonsmooth, nonseparable regularizer whose matrix $A^{h}$ is sparse;

• •

the dual SVM with intercept has a single linear nonseparable constraint.

For Algorithm 2, we set the restart with a variable sequence as in [16]. We did not tune the algorithmic parameters for each instance. We evaluate the precision of a primal-dual pair $(x,y)$ as follows. We define the smoothed gap [59] as

[TABLE]

and we choose the positive parameters $\beta$ and $\gamma$ as

[TABLE]

It is shown in [59] that when $\mathcal{G}_{\beta,\gamma}(x,y,\zeta,\omega)$, $\beta$ and $\gamma$ are small then the objective value and feasibility gaps are small.

For the Lasso problems, we compare our implementations of the coordinate descent method with scikit-learn’s coordinate descent [42] and CVXPY’s augmented Lagrangian method [12] called OSQP. As in Table 6, we have a factor 4 between Alg. 1 without screening and scikit-learn. Acceleration and screening allows us to reduce this gap without sacrificing generality. OSQP is efficient on small problems but is not competitive on larger instances.

For the TV-regularized regression problems, we compare ourself with FISTA where the proximal operator of the total variation is computed inexactly and with LBFGS where the total variation is smoothed with a decreasing smoothing parameter. Those two methods have been implemented for [13]. They manage to solve the problem to an acceptable accuracy in a few hours. As the problem is rather large, we did not run OSQP on it. For our methods, as $h$ is nonzero, we cannot use variable screening with the current theory. Alg. 1 quickly reduces the objective value but fails to get a high precision in a reasonable amount of time. On the other hand, Alg. 2 is the quickest among the four solvers tested here.

The third problem we tried is dual support vector machine with intercept. A very famous solver is libsvm [10], which implements SMO [43], a 2-coordinate descent method that ensures the feasibility of the constraints at each iteration. The conclusions are similar to what we have seen above. The specialized solver remains a natural choice. OSQP can only solve small instances. Alg. 1 has trouble finding high accuracy solutions. Alg. 2 is competitive with respect to the specialized solver.

9 Conclusion

This paper introduces a generic coordinate descent solver. The technical challenge behind the implementation is the fundamental need for a compiled language in the low-level operations that are partial derivative and proximal operator computations. We solved it using pre-defined atom functions that are combined at run time using a python interface.

We show how genericity allows us to decouple algorithm development from a particular application problem. As an example, our software can solve at least 12 types of optimization problems on large instances using primal-dual coordinate descent, momentum acceleration, restart and variable screening.

As future works, apart from keeping the algorithm up to date with the state of the art, we plan to bind our solver with CVXPY in order to simplify further the user experience.