Distributed Learning with Sparse Communications by Identification

TL;DR

This paper introduces an asynchronous distributed optimization algorithm that reduces communication costs through random sparsification of updates, achieving linear convergence and enabling automatic dimension reduction.

Contribution

It presents a novel asynchronous algorithm with sparse communication that converges linearly and identifies optimal sparse solutions for large-scale distributed learning.

Findings

Algorithm converges linearly for strongly convex functions.

Effectively identifies optimal sparse solutions.

Enables automatic dimension reduction in distributed settings.

Abstract

In distributed optimization for large-scale learning, a major performance limitation comes from the communications between the different entities. When computations are performed by workers on local data while a coordinator machine coordinates their updates to minimize a global loss, we present an asynchronous optimization algorithm that efficiently reduces the communications between the coordinator and workers. This reduction comes from a random sparsification of the local updates. We show that this algorithm converges linearly in the strongly convex case and also identifies optimal strongly sparse solutions. We further exploit this identification to propose an automatic dimension reduction, aptly sparsifying all exchanges between coordinator and workers.

Equations

$${\mathbf{prox}}_r(x)=\underset{y\in {\mathbb{R}}^d}{\mathrm{argmin}}\left\{r(y)+\frac{1}{2}\Vert y-x{\Vert }^2\right\}.$$

$$\displaystyle\left(\mathbf{prox}_{\lambda\|\cdot\|_{1}}(u)\right)_{[j]}=\left\{\begin{array}[]{ll}~{}0&\text{ if }u_{[j]}\in[-\lambda,\lambda]\\ ~{}u_{[j]}-\lambda&\text{ if }u_{[j]}>\lambda\\ ~{}u_{[j]}+\lambda&\text{ if }u_{[j]}<-\lambda\\ \end{array}\right.\qquad\text{for all $j$},$$

$$\underset{x\in {\mathbb{R}}^d}{\min }F(x)=\sum _{i=1}^M{\alpha }_i{f}_i(x)+r(x),$$

$${\alpha }_i=|{\mathcal{S}}_i|/n\text{and}\sum _{i=1}^M{\alpha }_i=1.$$

$$k_{m+1}$$

$${\Vert x^k-x^{\star }\Vert }^2\le {\Vert {\bar{x}}^k-{\bar{x}}^{\star }\Vert }^2\le {\left(1-\gamma \mu \right)}^{2m}\underset{i}{\max }{\Vert x_i^0-x_i^{\star }\Vert }^2,$$

$${\Vert x^k-x^{\star }\Vert }^2\le {\left(\frac{1-{\kappa }_{\text{eqref}eq:pb}}{1+{\kappa }_{\text{eqref}eq:pb}}\right)}^{2m}\underset{i}{\max }{\Vert x_i^0-x_i^{\star }\Vert }^2.$$

$$x_{i[j]}^k$$

$$x^k$$

$$\mathbb{P}[j\in {\mathbf{S}}^k]=p_j>0\text{for all j∈{1,…,d},}$$

$$\mathbb{E}{\Vert x^k-x^{\star }\Vert }^2\le {\left(p_{\max }(1-\gamma \mu {)}^2+1-p_{\min }\right)}^m\underset{i}{\max }{\Vert x_i^0-x_i^{\star }\Vert }^2,$$

$$\mathbb{E}{\Vert x^k-x^{\star }\Vert }^2\le {\left(p_{\max }{\left(\frac{1-{\kappa }_{\text{eqref}eq:pb}}{1+{\kappa }_{\text{eqref}eq:pb}}\right)}^2+1-p_{\min }\right)}^m\underset{i}{\max }{\Vert x_i^0-x_i^{\star }\Vert }^2.$$

$$\frac{p_{\min }}{p_{\max }}>(1-\gamma \mu {)}^2\stackrel{\gamma =\frac{2}{\mu +L}}{\ge }{\left(\frac{1-{\kappa }_{\text{eqref}eq:pb}}{1+{\kappa }_{\text{eqref}eq:pb}}\right)}^2$$

$$\displaystyle\mathbb{P}[j\in\mathbf{S}^{k}]=\left\{\begin{array}[]{cr}\pi&\text{ if }x^{k}_{[j]}=0\\ 1&\text{ if }x^{k}_{[j]}\neq 0\end{array}\right.~{}~{}~{}~{}~{}\text{for all $j\in\{1,\dots,d\}$ ~{}and~{} }\pi\in(0,1].$$

$${\kappa }_{\text{eqref}eq:pb}>{\kappa }_{\min }:=\frac{1-\sqrt{\pi }}{1+\sqrt{\pi }}.$$

$$h_{i,\ell }=f_i+\frac{\rho }{2}\Vert \cdot -x_{\ell }{\Vert }_2^2$$

$$\underset{x\in {\mathbb{R}}^d}{\min }{H}_{\ell }(x):=\sum _{i=1}^M{\alpha }_i\underset{h_{i,\ell }(x)}{\underbrace{\left(f_i(x)+\frac{\rho }{2}\Vert x-x_{\ell }{\Vert }_2^2\right)}}+r(x).$$

$$\displaystyle\kappa_{\eqref{eq:reco}}=\frac{\mu+\rho}{L+\rho}\,~{}~{}\Big{(}\geq~{}\kappa_{\eqref{eq:pb}}=\frac{\mu}{L}\Big{)}.$$

$${\mathbf{prox}}_{F/\rho }(x_{\ell })=\underset{x\in {\mathbb{R}}^d}{\mathrm{argmin}}\underset{=F(x)}{\underbrace{\sum _{i=1}^M{\alpha }_i{f}_i(x)+r(x)}}+\frac{\rho }{2}\Vert x-x_{\ell }{\Vert }_2^2.$$

$$x_{\ell +1}\approx {\mathbf{prox}}_{F/\rho }(x_{\ell })$$

$$\displaystyle p_{j,\ell}=\left\{\begin{array}[]{cr}\pi_{\ell}:=\min\left(\displaystyle\frac{c}{|\mathrm{null}(x_{\ell})|};1\right)&\text{ if }(x_{\ell})_{[j]}=0\\ 1&\text{ if }(x_{\ell})_{[j]}\neq 0\end{array}\right.~{}~{}~{}~{}\text{for all $j\in\{1,\dots,d\}$}.$$

$$\pi :=\frac{c}{d}\le \underset{\ell }{\inf }{\pi }_{\ell }.$$

$${\kappa }_{\text{eqref}eq:reco}=\frac{\mu +\rho }{L+\rho }>{\kappa }_{\min }⟺\rho >\frac{{\kappa }_{\min }L-\mu }{1-{\kappa }_{\min }}\text{ with }{\kappa }_{\min }=\frac{1-\sqrt{\pi }}{1+\sqrt{\pi }}\text{ as in eqrefeq:kappamin}.$$

$$\rho =\frac{{\kappa }_{\text{eqref}eq:reco}L-\mu }{1-{\kappa }_{\text{eqref}eq:reco}}\text{ with }{\kappa }_{\text{eqref}eq:reco}=\frac{1-\sqrt{\pi -\alpha }}{1+\sqrt{\pi -\alpha }}.$$

$$\left({\left(\frac{1-{\kappa }_{\text{eqref}eq:reco}}{1+{\kappa }_{\text{eqref}eq:reco}}\right)}^2+1-{\pi }_{\ell }\right)=\left(\pi -\alpha +1-{\pi }_{\ell }\right)=\underset{=:(1-{\alpha }_{\ell })}{\underbrace{(1-\alpha -({\pi }_{\ell }-\pi ))}}\le 1-\alpha <1.$$

$$\text{Set }\rho =\frac{\kappa L-\mu }{1-\kappa }\text{ and }\gamma \in \left(0,\frac{2}{\mu +L+2\rho }\right]\text{ with }\kappa =\frac{1-\sqrt{\pi -\alpha }}{1+\sqrt{\pi -\alpha }};\pi =\frac{c}{d}\text{ and }\alpha =\frac{c}{2d}.$$

$$\displaystyle p_{j,\ell}=\left\{\begin{array}[]{cr}\pi_{\ell}:=\min\left(\frac{c}{|\mathrm{null}(x_{\ell})|};1\right)&\text{ if }[x_{\ell}]_{j}=0\\ 1&\text{ if }[x_{\ell}]_{j}\neq 0\end{array}\right.~{}~{}~{}~{}~{}\text{for all $j\in\{1,\dots,d\}$}.$$

$$x_{\ell +1}\approx {\mathbf{prox}}_{F/\rho }(x_{\ell })=\underset{x\in {\mathbb{R}}^d}{\mathrm{argmin}}\left\{\sum _{i=1}^M{\alpha }_i\underset{h_{i,\ell }(x)}{\underbrace{\left(f_i(x)+\frac{\rho }{2}\Vert x-x_{\ell }{\Vert }_2^2\right)}}+r(x)\right\}$$

$${\mathsf{M}}_{\ell }=\lceil \frac{(1+\delta )\log (\ell )}{\log \left(\frac{1}{1-\alpha +\pi -{\pi }_{\ell }}\right)}+\frac{\log \left(\frac{2\mu +\rho }{(1-\delta )\rho }\right)}{\log \left(\frac{1}{1-\alpha +\pi -{\pi }_{\ell }}\right)}\rceil \text{ epochs.}$$

$$\Vert x_{\ell +1}-{\mathbf{prox}}_{F/\rho }(x_{\ell }){\Vert }^2\le \frac{(1-\delta )\rho }{(2\mu +\rho ){\ell }^{1+\delta }}\Vert x_{\ell }-{\mathbf{prox}}_{F/\rho }(x_{\ell }){\Vert }^2.$$

$$\Vert x_{\ell +1}-{\mathbf{prox}}_{F/\rho }(x_{\ell }){\Vert }^2\le \frac{\rho }{4(2\mu +\rho ){\ell }^{2+2\delta }}\Vert x_{\ell +1}-x_{\ell }{\Vert }^2.$$

Full text

1. Introduction

1.1. Context and contributions

We consider a distributed learning setup where observations are scattered over many machines. We assume that the machines, also referred to as workers, have a private subset of the examples, perform their computations independently, and communicate with a coordinator machine. Standard learning approaches would consider that the entire training set is moved and stored in one single machine or in a datacenter. In contrast, we consider the approach consisting in learning collaboratively a shared prediction model without moving the training data, as for instance in federated learning (see e.g. the recent review articles [22] and [29]). In this context, communications are typically the practical bottleneck of the learning process; see e.g. [28, 24, 34, 68, 37].

Recently, several approaches to mitigate this communication issue have been proposed, including different forms of gradient or model compression, which reduces representation size by quantization (reduce precision on elements) or sparsification (send only most significant gradient entries); we refer e.g. to early works [55, 61] and more recent developments [23, 59, 1, 4, 66] among many others.

In this paper, we address the question of reducing the size of communications with an approach complementary to existing compression techniques. We propose to adaptively sparsify the model updates using its progressively uncovered sparse structure. More precisely, we first propose an asynchronous distributed algorithm featuring a sparsification of upwards communications (worker-to-coordinator) by randomly zeroing the local update entries. This randomized technique is adjustable to various levels of communication costs, machines’ computational powers, and data distribution evenness. Also, an attractive and original property of this algorithm is the possibility to use a fixed learning rate that does not depend neither on communication delays nor on the number of machines. However, this random sparsification technique provably works only for i.i.d. sparsifications with either almost-uniform distributions, or well-conditioned problems. This makes aggressive sparsification or adaptation to the sparsity structure of the model impossible with such an algorithm.

To tackle this issue, we use a proximal reconditioning scheme wrapping up the previously mentioned algorithm as an inner minimization method. This allows us to perform much more aggressive sparsifications. Furthermore, we show that when using a sparsity-inducing regularizer, our reconditioned algorithm generates iterates that identify the optimal sparsity pattern of the model in finite time. This progressively uncovered pattern can be used to adaptively update the sparsification distribution of the inner method. All in all, this method only features sparse communications: the downward communications (coordinator-to-worker) consists in sending the (eventually) sparse model, and the the upwards communications (worker-to-coordinator) are adaptively and aggressively sparsified.

Finally, we show theoretically and numerically that our method has better performance than its non-sparsified version, in terms of suboptimality with respect to the quantity of information exchanged.

1.2. Related literature on sparsification

A popular way of sparsifying models is to use proximal methods: they consist in using a sparsity-inducing penalty $r$ (typically $\ell_{1}$ or $\ell_{1,2}$ norms, see e.g. [3]), and handling it with a proximity operator defined as the $\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ mapping

[TABLE]

For many popular regularizations, this operator has an explicit formulation or can be efficiently computed. For example, the proximity operator of the $\ell_{1}$-norm (with parameter $\lambda$) is given coordinate-wise111We denote by $x_{[j]}\in\mathbb{R}$ the $j$-th coordinate of vector $x$. as

[TABLE]

which has a clear sparsification effect. This kind of proximal sparsification was successfully used in machine learning (see e.g. the seminal work [12]), for algorithmic algorithmic perspectives (see e.g. [25]) and in the case of synchronous distributed systems (see e.g. [67, 56]).

Block coordinate descent algorithms, which knew a rebirth in the context of huge-scale learning [43, 51], can also be interpreted as specific sparsification of gradient updates. In the context of $\ell_{1}$-regularized problems, they can be combined with screening techniques for efficiently solving high-dimensional problems [39]. We also mention that the idea of randomly selecting some entries has been already applied in distributed algorithms: random selection is used to sparsify local gradients in the synchronous algorithm of [68], to sparsify the variance-reducing term in the stochastic algorithm of [47], or to sparsify updates in fixed-point iterations [48].

Our proposed algorithm combines these three approaches (proximal methods, coordinate descent, and random selection): the output of a proximal step is sent downwards (coordinator-to-workers) and random coordinate updates are sent upwards (workers-to-coordinator); see details in Section 3. Beyond the original combination of known ingredients, the specificity of our work is the study of the adaptive sparsification and its interplay with identification properties.

Identification of optimal structure has been extensively studied in the context of constrained convex optimization (e.g. [69]) and nonsmooth nonconvex optimization (e.g. [19]). Although many deterministic algorithms have been proved to have identification properties, the situation is less clear for randomized algorithms: in particular (proximal) stochastic gradient methods are known to be unable to identify substructure (see e.g. [26]), while it has been recently proved that (proximal) variance-reduced stochastic gradient methods do have such identification properties [50, 13]. No identification results have yet been reported for asynchronous distributed optimization algorithms.

For our algorithm, we use identification to show that the iterates produced by the coordinator become sparse when the regularization $r$ is a sparsity-inducing norm [3] and also leverage the sparsity pattern identified by the coordinator to adapt the sparsification at the workers. Finally note that, though identification has already been used to accelerate algorithms (by e.g. a better tuning of parameters [30]), it is the first time that it is used for an automatic exchange reduction and sparsification.

2. Preliminaries: problem, notations, and recalls

We consider a distributed learning setup where $n$ observations are scattered over $M$ machines; each machine $i$ having a private subset $\mathcal{S}_{i}$ of the examples. This leads to optimization problems of the form

[TABLE]

with $f_{i}(x)=\frac{1}{|\mathcal{S}_{i}|}\sum_{j\in\mathcal{S}_{i}}l_{j}(x)$ the local empirical risk at machine $i$ ($l_{j}$ standing for the smooth loss function for example $j$). The $\alpha_{i}$ are the proportion of observations locally stored in machine $i$, hence:

[TABLE]

The function $r$ is a regularization term introducing a prior on the structure of the model $x$, improving both its interpretability and its stability. One of the first and most well-known instances of such a prior is the sparsity sought in the parameters of least-squares regression with $\ell_{1}$-regularization $r=\|{\cdot}\|_{1}$), leading to the well-known lasso problem [64]. The notation related to sparsity is the following: we denote the support of a vector $x\in\mathbb{R}^{d}$ by $\mathrm{supp}(x)=\{i\in\{1,\dots,d\}:x_{[i]}\neq 0\}$, and its size by $|\mathrm{supp}(x)|$. Oppositely, we define the sparsity pattern of a vector $x\in\mathbb{R}^{d}$ by the set $\mathrm{null}(x)=\{i\in\{1,\dots,d\}:x_{[i]}=0\}$, and denote its size by $|\mathrm{null}(x)|$. Note that these two index sets are complementary; $\mathrm{supp}(x)=\overline{\mathrm{null}(x)}$.

We consider the problem ($\mathsf{P}$) in a standard convex setting, with the function $r$ convex and lower-semicontinuous (lsc), and all the functions $(f_{i})$ $L$-smooth and $\mu$-strongly convex (with $\mu\geq 0$). We then define the condition number of the smooth part of ($\mathsf{P}$) as $\kappa_{\eqref{eq:pb}}=\mu/L.$

On asynchronous distributed optimization

The general computation framework is the following: one machine gets the current model from the coordinator, improves it from its own data which produces a local model update; only this update is sent to the coordinator, which computes a new shared model. There exist many methods for solving problems of this type by iteratively sending small messages or model updates during the training process (as opposed to sending raw data); see e.g. [70, 2, 48]. The literature on distributed optimization methods without shared memory usually use on restrictive assumptions on the computing system delays, which in turn impact the obtained convergence rates. Asynchronous coordinate descent methods which are able to handle unbounded delays were recently proposed in [18, 62] but use decreasing stepsizes. In contrast, the asynchronous proximal-gradient algorithm of [41, 40] allows the use of fixed stepsizes (the same as in the vanilla proximal gradient algorithm), and their delay-independent analysis technique result an linear convergence in the strongly convex case. We present below the notation and the asynchronous algorithm of [41], which is going to be the baseline for our developments.

Notation: time, delays, and epochs

As asynchronous distributed setup allows the algorithm to carry on computations without waiting for slower machines: the machines work on outdated versions of the main variable, and the coordinator gathers the workers’ inputs on the fly to produce the updates. The following notation formalizes this framework:

• •

For the coordinator. We define the time $k$, as the number of updates the coordinator has received from any of the workers. At time $k$, the coordinator receives some input from an agent, denoted by $i^{k}$, updates its global variable, and sends it to worker $i^{k}$.

• •

For worker $i$. At time $k$, we introduce $d_{i}^{k}$ as the time elapsed since the last time the coordinator received an update from worker $i$. In particular, we have $d_{i_{k}}^{k}=0$. We also consider $D_{i}^{k}$ as the elapsed time since the penultimate update. This means that, at time $k$, the last two moments when the worker $i$ updates the coordinator are $k-d_{i}^{k}$ and $k-D_{i}^{k}$. .

We consider a totally asynchronous as per the classification of the book [7, Chap. 6.1], meaning that all workers are responsive. However the worker responses can have arbitrary delays. Following [40] we introduce the notion of epochs, as the sequence of stopping times $(k_{m})$ defined as $k_{0}=0$ and

[TABLE]

This sequence fully captures the asynchronicity and allows us to analyze algorithms independently of the computing system. All our convergence results involve this notion of epochs.

Asynchronous proximal-gradient and its rate for strongly convex objectives

As explained above, our baseline is the asynchronous proximal-gradient algorithm of [41], called DAve-PG, which solves distributed problems of the form ($\mathsf{P}$) represented by the triplet $((\alpha_{i}),(f_{i}),r)$ for (workers weights, workers functions, global regularization); see Algorithm 1. This algorithm converges linearly at a rate that only depends on the function properties but neither on the number of machines nor on the delays, directly embedded in the sequence $(k_{m})$.

Theorem 2.1** (Th. 3.2 of [40]).**

Let the functions $(f_{i})$ be $\mu$-strongly convex ($\mu>0$) and $L$-smooth. Let $r$ be convex lsc. Using $\gamma\in(0,\frac{2}{\mu+L}]$, DAve-PG converges linearly on the epoch sequence $(k_{m})$. More precisely, for all $k\in[k_{m},k_{m+1})$

[TABLE]

with the shifted local solutions $x_{i}^{\star}=x^{\star}-\gamma_{i}\nabla f_{i}(x^{\star})$ and $\bar{x}^{\star}=\sum\alpha_{i}x_{i}^{\star}$.

Furthermore, using the maximal stepsize $\gamma=\frac{2}{\mu+L}$, we obtain for all $k\in[k_{m},k_{m+1})$

[TABLE]

In this result, the stepsize $\gamma$ can be taken in the usual range for (proximal) gradient descent (in contrast with existing asynchronous algorithms, even for the works with the most realistic assumptions; see e.g. [2, 35, 62] and references therein).

3. Sparsification of upward communications

In this section, we present a first distributed algorithm for solving ($\mathsf{P}$) with sparse workers-to-coordinator communications.

3.1. Random sparsification of local updates

In the proposed method, the coordinator machine asynchronously gathers sparsified delayed updates from workers and sends back the global variable. More specifically, each worker independently compute a gradient step on its local loss for a randomly drawn subset of coordinates only. The coordinator machine keeps track of the weighted average of the most recent worker outputs, computes the proximity operator of the regularizer at this average point, and sends this result back to the updating worker $i^{k}$. At iteration $k$, the random subset of entries that worker $i^{k}$ updates is denoted by $\mathbf{S}^{k}$ (in bold, emphasizing that it is the only random variable in the algorithm). The update writes

[TABLE]

With this sparsification, the local updates correspond to a random block coordinate descent step for the workers. However, this algorithm does not boil down to an asynchronous stochastic block-coordinate descent algorithm such as [35, 62, 48, 51], since our method maintains a variable, $\bar{x}^{k}$, aggregating asynchronously all the workers’ contributions.

Assumption 1* (On the random sparsification).*

The sparsity selectors $(\mathbf{S}^{k})$ are independent and identically distributed random variables. We select a coordinate in $\mathbf{S}^{k}$ as follows:

[TABLE]

with ${p}=(p_{1},\dots,p_{d})\in(0,1]^{d}$. We denote $p_{\max}=\max_{i}p_{i}~{}\text{and}~{}p_{\min}=\min_{i}p_{i}.$

The selectors $(\mathbf{S}^{k})$ being the only random variables of the algorithm, it is natural to define the filtration $\mathcal{F}^{k}=\sigma(\{\mathbf{S}^{\ell}\}_{\ell<k})$ so that all variables at time $k$ ($x_{i}^{k}$, $\bar{x}^{k}$, $x^{k}$, $d_{i}^{k}$, $D_{i}^{k}$) are $\mathcal{F}^{k}$-measurable but $\mathbf{S}^{k}$ is not.

3.2. Distributed implementation

The algorithm introduced in the previous section can be naturally distributed. We formalize this in Algorithm 1, that we call Spy. Similarly to DAve-PG, this method is generic in the sense that none of its ingredients, including the stepsize choice, depend on the computing system. It also shares the feature that although each coordinator update involves only one agent (and thus part of the data), all the data is always implicitly involved in the coordinator variable; which allows the algorithm to cope with the heterogeneity of the computing system (data distribution, agents delays). Its presentation uses the following notation: for a vector of $x\in\mathbb{R}^{d}$ and a subset $S$ of $\{1,\dots,d\}$, $[x]_{S}$ denotes the sparse size-$d$ vector where $S$ is the set of non-null entries, for which they match those of $x$, i.e. $([x]_{S})_{[i]}=x_{[i]}$ if $i\in S$ and [math] otherwise.

The communications per iteration are (i) a blocking send/receive from a worker to the coordinator (in blue) of size $|\mathbf{S}|$, and (ii) a blocking send/receive from the coordinator to the last updating worker (in red) of the current iterate. The upward communication is thus made sparse by the algorithm and the downward communication cost depends on the structure of $x^{k}$, which is the output of a proximal operator on $r$. In the case of $\ell_{1}$-regularization, $x^{k}$ will become sparse after some iterations, leading to a two-way sparse algorithm; we will be discussed in Section 5.1.

3.3. Convergence Analysis

We study the convergence properties of Spy under standard assumptions on ($\mathsf{P}$) and no apriori assumption on the computing system or on data distribution. We emphasize here that we do not put assumptions on the delays; for instance they do not to need to be bounded or independent of the previous selectors $(\mathbf{S}^{k})$. Let us also notice that the local iterates $(x_{i})$ do not converge to a minimizer of the individual functions $(f_{i})$ but rather to local shifts of the solution $x^{\star}$ of ($\mathsf{P}$) (unique from the strong convexity assumption).

Theorem 3.1** (Reaches and Limits of Sparsification).**

Let the functions $(f_{i})$ be $\mu$-strongly convex ($\mu>0$) and $L$-smooth. Let $r$ be convex lsc. Take $\gamma\in(0,\frac{2}{\mu+L}]$. Suppose that Assumption 1 holds for the probability vector ${p}$, and that $\frac{p_{\min}}{p_{\max}}\geq(1-\gamma\mu)^{2}.$

Then, Spy on $((\alpha_{i}),(f_{i}),r~{};~{}p)$ verifies for all $k\in[k_{m},k_{m+1})$

[TABLE]

with the shifted local solutions $x_{i}^{\star}=x^{\star}-\gamma_{i}\nabla f_{i}(x^{\star})$.

Furthermore, using the maximal stepsize $\gamma=\frac{2}{\mu+L}$, we obtain for all $k\in[k_{m},k_{m+1})$

[TABLE]

Proof.

The co-existence of both deterministic and stochastic delays in the algorithm calls for an original mathematical analysis using the notation introduced in Section 2, and in particular the notion of epoch sequence (2.1). This is detailed in Supplement SM2. ∎

This result establishes bounds that lead to convergence whenever the selection probabilities are well chosen. First, if all probabilities are equal to $1$, the algorithm boils down to DAve-PG and Theorem 3.1 coincides with Theorem 2.1. In more general cases, this result has to be interpreted more carefully as developed in the following section.

3.4. On the sparsification choice

In the totally distributed setting, all machines are responsive, which means with our notation: $m\to\infty$ when $k\to\infty$. Then, Theorem 3.1 gives linear convergence of the mean squared error in terms of epochs if

[TABLE]

Thus the behavior of the algorithm depends on if the selection is performed uniformly (and thus structure-blind) or non-uniformly (to encompass some prior information). We distinguish these two cases below.

Inefficiency of uniform sparsification

If the selection is uniform, i.e. $p_{i}=\pi\in(0,1]$ for all $i$, we directly get convergence from (3.2) as the mean squared error vanishes linearly in terms of epochs with a rate $(1-\pi\gamma\mu(2-\gamma\mu))$, degraded compared to the $(1-\gamma\mu)$ rate of DAve-PG. Unfortunately, such uniform selection also results in poor performance in many cases (as illustrated in Supplement SM1). Thus, when Spy is used with uniform sampling, it has a degraded performance compared to (non-sparsified) DAve-PG, both in theory and in practice. Adaptivity is key for sparsifying efficiently.

Efficiency of adaptive sparsification

There exist several adaptation strategies for selecting coordinates. In the context of coordinate descent, we can mention greedy coordinates selection [11, 45], other heuristics [36, 15], and importance sampling [71, 52], that can sometimes be used in practice. The idea of adaptive coordinate descent methods based on the coordinate-wise Lipschitz constants and current values of the gradient is proposed in [49, 42, 60]. The idea of adaptively using the iterates structure enforced by a non-smooth regularizer was recently developed in [16].

We use here the technique of [16] specialized for the coordinate selection: when some coordinates get null, there is some hope that they will remain null for subsequent iterations, and it is thus natural to update preferentially the non-null coordinates. Mathematically, this means that we select coordinates in the active support as follows:

[TABLE]

In words, we communicate the coordinates in the support of the coordinator point $x^{k}$, together with coordinates outside the support, randomly selected with some exploration probability $\pi$.

This adaptive sampling often shows tremendous gains in practice compared to uniform sampling; however, it may not converge in some situations. This is due to two technical points related to Theorem 3.1:

• •

A good conditioning is necessary to allow for a small $\pi=p_{\min}$ (with $p_{\max}=1$). More precisely, from (3.4), we get that the minimal conditioning to allow for a probability $\pi$ of selection outside the support is:

[TABLE]

Since we aim at taking $\pi$ small to communicate little, this is a stringent condition.

• •

The sampling is not i.i.d. anymore since the probabilities depend on the points generated by the algorithm. This difficulty could be tackled with a small algorithmic fix that unfortunately degrades practical performances and a refined analysis (see the discussion to Supplement SM3), but the upcoming methods directly address this point.

While these two issues appear separate, they can both be overcome by iteratively reconditioning the problem.

4. Proximal reconditioning for adaptive sparsification

The learning problem ($\mathsf{P}$) should be well-conditioned to safely apply the random sparsification technique of Spy with reasonable exploration probability $\pi$. The idea is then not to apply Spy directly to ($\mathsf{P}$) but rather to a modified problem for which we control the condition number and the sparsification potential.

We thus propose to recondition the learning problem ($\mathsf{P}$) using the standard proximal algorithm (see e.g. [54]), as described in Section 4.1. We present in Section 4.2 our algorithmic choices and the resulting algorithm, called Reconditioned–Spy.

4.1. Proximal reconditioning

This type of methods consist in iteratively regularizing the problem with the squared distance to some center point. We call outer iteration the process of (approximatively) solving such a reconditioned problem. At outer loop222The quantities related to outer loop $\ell$ are denoted with a subscript $\ell$. $\ell$, we define the worker $i$’s regularized function as

[TABLE]

where $\rho$ is the regularization factor and $x_{\ell}$ the center point at outer loop $\ell$. The reconditioned problem for loop $\ell$ then writes

[TABLE]

For $\mu$-strongly convex $L$-smooth $(f_{i})$, the regularized functions $h_{i,\ell}$ are $(\mu+\rho)$-strongly convex and $(L+\rho)$-smooth. Hence, the condition number of the smooth part of ($\mathsf{R}_{\ell}$) writes

[TABLE]

The optimal solution of ($\mathsf{R}_{\ell}$) is exactly the proximal point (1.1) of $F/\rho$ at $x_{\ell}$

[TABLE]

Thus, for solving ($\mathsf{P}$), each outer iteration consists in an (inexact) proximal step:

[TABLE]

The proximal point algorithm is a standard regularization approach in optimization. It was presented in [6, Chap. 5] to recondition a convex quadratic objective, for which computing the proximal operator (1.1) is easy (it is the unique solution of a linear system, well-conditioned by construction). The general proximal algorithm was then popularized by the seminal works [38, 54]. The study of these algorithms, and especially their inexact variants, has attracted a lot of attention; see e.g. [17, 57, 14, 32, 33].

Practical implementations of such methods require an inner algorithm to compute the proximal point (we will use Spy here) and a rule to stop this algorithm. Several papers consider the key question of inner stopping criteria for inexact proximal methods in various contexts; see e.g. [14] in smooth optimization, [27] in nonsmooth optimization, and [58] in operator theory. In this paper, we use the standard criteria, following [54], which gave the first inexact rules for the proximal point algorithm in the context of monotone operators.

4.2. Reconditioned–Spy

We present our main algorithm, which consists in applying the inexact proximal scheme (4.1) with Spy as inner algorithm to solve ($\mathsf{P}$).

At the outer iteration $\ell$, we run Spy for solving ($\mathsf{R}_{\ell}$) with i.i.d. non-uniform sparsification probabilities given, for a fixed $0<c\leq d$, by

[TABLE]

The sparsification level over outer iterations is then bounded from below by

[TABLE]

We now choose the reconditioning parameter $\rho$ from $\pi$ so that Spy converges linearly to the solution of the reconditioned problem ($\mathsf{R}_{\ell}$). We know, from Section 3.4, that this is the case as soon as

[TABLE]

To properly handle the strict inequality above, we propose to choose a conditioning which guarantees a $(1-\alpha)$ rate for Spy on the reconditioned problems uniformly over $\ell$. Mathematically, for $0<\alpha<\pi$ (for instance $\alpha=\pi/2$), we choose

[TABLE]

Then, the contraction factor of Spy with the maximal stepsize (see (3.3)) for the reconditioned problem ($\mathsf{R}_{\ell}$) becomes

[TABLE]

This means that Spy is linearly convergent on the reconditioned problem ($\mathsf{R}_{\ell}$). Thus, it can safely be used as an inner method in the inexact proximal algorithm (4.1) to solve the original problem ($\mathsf{P}$).

The remaining part is to the choice of a stopping criterion for the inner loop. We propose to use three different criteria: epoch budget, absolute accuracy, and relative accuracy (called $\mathsf{C}_{1},\mathsf{C}_{2}$, and $\mathsf{C}_{3}$ respectively). Stopping criteria based on accuracy are usually more stringent to enforce (see e.g. [32, Sec. 2.3] and references therein), however they may bring significant performance improvement when the instantaneous rate is better than the theoretical one.

The resulting algorithm, called Reconditioned–Spy, is presented as Algorithm 1. Under any of the three stopping criteria, we recover the same convergence result, formalized in the next theorem.

Theorem 4.1**.**

Let the functions $(f_{i})$ be $\mu$-strongly convex ($\mu\geq 0$) and $L$-smooth. Let $r$ be convex lsc. If $\mu=0$ and $\mathsf{C}_{3}$ is used, we furthermore require that $F$ has a unique minimizer $x^{\star}$ and that $\liminf_{x\to x^{\star}}(F(x)-F(x^{\star}))/\|x-x^{\star}\|^{2}>0$.

*Then, the sequence generated by Reconditioned–Spy on $((\alpha_{i}),(f_{i}),r)$ with stopping criterion $\mathsf{C}_{1},\mathsf{C}_{2}$, or $\mathsf{C}_{3}$ converges almost surely to a minimizer of $F$. Furthermore, if $\mu>0$, then we have333We use the standard notation: $a_{\ell}=\tilde{\mathcal{O}}((1-r)^{\ell})$ denotes that there exists $C,p$ such that $a_{\ell}\leq C\ell^{p}(1-r)^{\ell}$. *

[TABLE]

Proof.

Even though the final result is similar for the three cases, the proof techniques are rather different. We thus present them separated in Supplement SM4. ∎

This result thus establishes that Reconditioned–Spy converges linearly to a solution of ($\mathsf{P}$). This means that Reconditioned–Spy has qualitatively the same behavior as DAve-PG, with the additional feature of having sparse local updates and therefore sparse upward communications. In other words, our algorithm is similar to the baseline in terms of iterations, but it is expected to be faster in terms of communications (more precisely in terms of quantity of information exchanged between coordinator and workers) which would result in a wallclock gain in practice, as shown in Section 6. Before this, we further investigate in the next section the theoretical gain of our sparsification technique in the case of sparse optimal solutions.

Remark 1* (Acceleration (with respect to iterations)).*

In this paper, we are interested in sparsifying communications and we primarily consider the reconditioning aspect of proximal methods, leaving aside other aspects including acceleration. As proposed in [17], the iterations of the inexact proximal algorithm can indeed be accelerated using Nesterov’s method [44]. The recent works [32, 33] also propose accelerated and quasi-Newton variants of the inexact proximal point algorithm as a meta-algorithm to improve the convergence of optimization methods (driven by machine learning applications [31]). Here, we investigate the complexity in terms of communications rather than iterations, so we do not insist much on these accelerated variants. The developments of this section could still be extended to accelerated proximal algorithm, following the meta-algorithm of [32]. We briefly study this direction in Supplement SM5. In particular, Theorem SM5.1 shows a gain of the square root on the rate for the accelerated version compared to Theorem 4.1, as usually observed with accelerated methods.

5. Identification for two-way sparse communications

In the previous sections, we present an adaptive sparsification of upward communications (worker-to-coordinator) and show that the resulting algorithm converges after proximal reconditioning. By construction, the downward communications (coordinator-to-workers) depends on the structure of $x^{k}$ (the coordinator point of the inner method), which is the output of a proximal operator on $r$. In the case of $\ell_{1}$-regularization or other sparsity-promoting regularization [3], we show in Section 5.1 that the $x^{k}$ eventually become sparse after some iterations. This automatically makes our algorithm a two-way sparse algorithm. Finally, in Section 5.2, we take a closer look to the complexity of our algorithm with respect to the communication cost.

For this study, we make an additional assumption that our problem has a strongly sparse solution. This assumption is divided into two parts: i) the regularizer $r$ should induce a stable support at the optimum (through its proximity operator); and ii) this optimal support $\mathrm{supp}(x^{\star})$ should be small with respect to the ambient dimension.

Assumption 2* (Strongly sparse optimal solution).*

Problem ($\mathsf{P}$) is $\mu$-strongly convex ($\mu>0$) and its solution $x^{\star}$ verifies

• i)

$\exists~{}\varepsilon>0$ such that $\mathrm{supp}(x^{\star})=\mathrm{supp}\left(\mathbf{prox}_{r}\left(x^{\star}-\sum_{i=1}^{M}\alpha_{i}\nabla f_{i}(x^{\star})+\mathbf{e}\right)\right)~{}\forall\mathbf{e}\in\mathcal{B}(0,\varepsilon)$;

• ii)

the size $s^{\star}=|\mathrm{supp}(x^{\star})|$ of the optimal support is small compared to $d$: $s^{\star}\ll d$.

While part ii) is rather explicit, part i) is quite abstract. For instance, when $r=\lambda\|\cdot\|_{1}$, using the explicit form of the proximity operator (1.5), part i) directly translates to

[TABLE]

This condition matches the nondegeneracy condition for sparse solutions commonly admitted for exact recovery in machine learning; see e.g. [46, 63]. The interest of the general assumption i) is that it accounts for a variety of sparsity-inducing regularizations, including weighted $\ell_{1}$-norms, “group” $\ell_{1}/\ell_{q}$-norms; see [3, Sec. 3.3].

5.1. Identification and consequences

The iterates of proximal algorithms usually identify the optimal structure; see e.g [65] or [21]. In the case of $\ell_{1}$-regularization, this means that proximal algorithms produce iterates that eventually have the same support as the optimal solution of ($\mathsf{P}$). Unfortunately, randomness may break this identification. For instance, it is well-known that for the proximal stochastic gradient descent, the sparse structure may not be identified with probability one; see e.g. [26] and a counter-example in [50]. We first establish that our algorithm444Note that a similar identification result holds for DAve-PG: in the same context, we have $\mathrm{supp}(x^{k})=\mathrm{supp}(x^{\star})$ for $k\geq k_{m}$ and $m$ larger than some threshold; see Supplement SM6.1. does identify the optimal support under the non-degeneracy assumption 2.

Theorem 5.1** (Identification).**

Let the functions $(f_{i})$ be $\mu$-strongly convex and $L$-smooth. Let $r$ be convex lsc. Under Assumption 2, the outer and inner iterates of Reconditioned–Spy identify the optimal structure in finite time: with probability one there exists a finite time $\Lambda<\infty$ such that

[TABLE]

where $x^{k}_{\ell}$ denotes the $k$-th iterate produced by Spy during the $\ell$-th outer loop.

Proof.

We proved in the previous section that Reconditioned–Spy converges almost surely, but it is not enough to guarantee identification in general555Take $d=1$, $F(x)=|x|$, and $x_{\ell+1}=\mathbf{prox}_{|\cdot|}(x_{\ell})+1/\ell^{2}$. The minimum of $F$ is [math] but we have $x_{\ell}=1/(\ell-1)^{2}>0$ for all $\ell>1$.. Here it is the fact that the inner algorithm Spy features a proximity operator with a non-vanishing stepsize that yields the following identification property. This proof is detailed in Supplement SM6.1. ∎

This identification has two consequences on communications in our distributed setting. First, identification implies that the variables communicated by the coordinator to the workers will eventually be sparse. Second, this sparsity is also leveraged in the sparsification strategies of Reconditioned–Spy where only the coordinates in $\hbox{}(x_{\ell})$ are randomly zeroed. Thus, for sparsity inducing problems such as $\ell_{1}$-regularized learning problems, our distributed algorithm has, structurally, two-way sparse communications.

Even better, once this identification occurs, the rate of the inner algorithm Spy dramatically improves to match the rate of its non-sparsified version DAve-PG. To get this improved rate, a small additional property is needed on the regularizer: $r$ has to be separable with respect to $\mathrm{supp}(x^{\star})$ i.e. $r(x)=r_{1}([x]_{\mathrm{supp}(x^{\star})})+r_{2}([x]_{\mathrm{null}(x^{\star})})$ which holds true for almost all sparsity inducing regularizations [3, Sec. 3.3].

Theorem 5.2** (Improved rate).**

Let the functions $(f_{i})$ be $\mu$-strongly convex and $L$-smooth. Let $r$ be convex lsc, separable with respect to $\mathrm{supp}(x^{\star})$. Under Assumption 2, the inner iterates of Reconditioned–Spy benefit from an improved rate after identification. There is $\Lambda<\infty$ such that for all $\ell>\Lambda$ and $k\in[k_{m},k_{m+1})$, using the maximal stepsize $\gamma=\frac{2}{\mu+L+2\rho}$,

[TABLE]

where $x^{k}_{\ell}$ denotes the $k$-th iterate produced by Spy during the $\ell$-th outer loop.

Proof.

The proof is reported to Supplement SM6.2. ∎

Thus our algorithm eventually has the practical interest of having sparse two-way communication, at almost no additional cost.

5.2. Communication complexity

We study in this section the asymptotic communication complexity of our method in terms of number of coordinates (real numbers) exchanged between the coordinator and the workers. To do so, we combine the number of outer iterations to reach an accuracy of $\varepsilon$ with the asymptotic communication cost of the inner iterations needed to reach the stopping criteria. Mathematically, we define

[TABLE]

where

• •

$\mathbf{c^{\mathrm{up}}}$ (resp. $\mathbf{c^{\mathrm{down}}}$) is the (expected) number of coordinates communicated from the coordinator to the active worker (resp. from the active worker to the coordinator) during one iteration and $K$ is the average number of iterations per epoch;

• •

$\mathbf{M}^{\mathsf{C}}$ is the (expected) number of inner epochs to reach stopping criterion $\mathsf{C}$;

• •

$\mathbf{L}(\varepsilon)$ is the number of outer loops to reach accuracy $\varepsilon$.

Focusing on the final regime of the algorithm when identification has taken place (as per Section 5.1), we get for Reconditioned–Spy

[TABLE]

which leads to the following communication complexity for Reconditioned–Spy.

Theorem 5.3** (Communication complexity of Reconditioned–Spy).**

Let the functions $(f_{i})$ be $\mu$-strongly convex and $L$-smooth ($\mu>0$), $r$ be convex lsc, separable with respect to $\mathrm{supp}(x^{\star})$. Let Assumption 2 hold. If the parameter $c$ is of the same order as $s^{\star}$ compared to $d$ ($c\approx s^{\star}\ll d$), then the communication complexity (5.1) of Reconditioned–Spy with criteria $\mathsf{C}_{2}$ or $\mathsf{C}_{3}$ is:

[TABLE]

Proof.

The proof consists in evaluating the terms in (5.1), one by one, in the right regime. It is given in Supplement SM6.3. ∎

In order to show the benefit in terms of communications of the proposed method, we can compare this result with our baseline DAve-PG. For DAve-PG, there is no inner loop so $\mathbf{M}^{\mathsf{C}}=1$ and Theorem 2.1 gives us $\mathbf{L}(\varepsilon)=\mathcal{O}((\mu+L)/\mu\log(1/\varepsilon))$. However, even if DAve-PG identifies which implies that $\mathbf{c^{\mathrm{up}}}=s^{\star}$ as previously, the cost of an upward communication is $\mathbf{c^{\mathrm{down}}}=d$ since there is no sparsification. This yields the following gain in communication complexity of our algorithm over DAve-PG ($1$ meaning similar performances; the greater, the better):

[TABLE]

The gain shows a product of three terms. The first one is greater than $1$ and depends on the conditioning; the second one is in $(0,1]$ but should be not far from $1$, provided that the final sparsity is not too poorly estimated. Finally, the last term fully exhibits the merits of adaptive sparsification with a term in $d+s^{\star}$ for DAve-PG which is much greater than the $\sqrt{ds^{\star}}$ for Reconditioned–Spy. This last term thus shows a nice dependence in the dimension of the problem and optimal solution for the proposed method. This theoretical gain of Reconditioned–Spy compared to DAve-PG is confirmed in the next numerical illustrations.

6. Numerical illustrations

In this section, we illustrate the communication gain provided by adaptive sparsification, on two classic $\ell_{1}$-regularized empirical risk minimization problems.

Problems. We consider two types of problems: lasso with randomly generated data and $\ell_{1}$-regularized logistic regression on popular datasets. We first consider a lasso problem

[TABLE]

with $d=1000$ features and an example set of size $m=500$. We take $A$ randomly generated from the standard normal distribution, $b=Ax_{0}+e$ where $x_{0}$ is a $99\%$ sparse vector and $e$ is taken from the normal distribution with standard deviation $0.01$. We take $\lambda_{1}$ to reach an optimal solution of size $|\mathrm{supp}(x^{\star})|=12$ (i.e. 1.2% of $d=1000$).

We also examine the regularized logistic regression with elastic net

[TABLE]

on two standard datasets from the LibSVM repository: madelon ($d=500$ $m=2000$) rcv1_train($d=47236$ $m=20242$). We take hyperparameters as follows: for madelon, $\lambda_{2}=0.001$ and $\lambda_{1}=0.03$ chosen to reach $99\%$ sparsity (i.e. $|\mathrm{supp}(x^{\star})|=5$); for rcv1, $\lambda_{2}=0.0001$ and $\lambda_{1}=0.001$ to reach $99.7\%$ sparsity (i.e. $|\mathrm{supp}(x^{\star})|=61$).

We distribute these problems on a machine with $32$ cores and $256$ Gb of RAM: one core plays the role the coordinator while $M$ cores act as workers ($M=5$ for LASSO problems, $M=10$ for madelon, and $M=20$ for rcv1). The data sets are split evenly between the $M$ workers.

Algorithms. We illustrate our sparsified algorithm Reconditioned–Spy (Algorithm 1) for different the amount of randomly chosen coordinates $c$. We take a simplified stopping criteria $\mathsf{C}_{1}$ with $\mathsf{M}_{\ell}=1$; we thus stop the inner iterations after two passes over the data, following the practical guidelines of Catalyst [32]. We observe that this simple stopping rule gives similar empirical convergence as $\mathsf{C}_{3}$ with respect to both iterations and scalars exchanged (without the additional computational cost of the test); for a numerical illustration, see Supplement SM7.

We display the performances of the algorithms in three ways:

• •

size of support vs number666Identification is illustrated plotting the size of $\mathrm{supp}(x^{k})$. For Reconditioned–Spy, $x^{k}$ refers to $x^{j}_{\ell}$ where there have been $k$ inner iterations during the $\ell-1$ first outer loops plus $j$ in the $\ell$-th inner loop. of inner iterations, showing the identification

• •

functional suboptimality vs number of inner iterations,

• •

functional suboptimality vs communication cost, modelled as the number of couples (coordinate, value) sent from and to the coordinator.

Empirical results. We make the following observations from experiments with both problems (reported in Figures 1 and 2 respectively). When the support of iterates is far from the optimal one, sparsification is generally bad for convergence in terms of iterations (as shown by the slopes in the plots “suboptimality vs iterations”), but even in terms of communications (see the beginning of the curves “suboptimality vs exchanges’). On the other hand, when the iterates begin to be closer to the optimal support, adaptive sparsification becomes highly beneficial as illustrated by the final slopes of the plots “suboptimality vs exchanges”.

Since there is no guarantee that the currently identified support is the optimal one, it is impossible to restrict ourselves to a subset of the coordinates; here comes the need for our adaptively sparsified method, that keeps exploring dimensions, additionally to those in the current support. The quantity of randomly chosen dimensions, controlled by $c$ has an slight impact: we see on that that (relatively) small and large values of $c$ (yellow and black curves) lead to slightly worse slopes on the convergence plots, compared to $c$ being in the range of 1 to 3 times the optimal support (green and pink curves) which is our recommendation both theoretically (see Theorem 5.3) and in practice.

Finally, in order to mitigate the above-mentioned negative effects of sparsification in the first iterations for large problems, we propose to use a warmstart strategy: in the first iterations, we use a non-sparsified (eg. DAve-PG) or a moderately sparsified method to allow for a sharp initial functional decrease, leading to some partial identification; after this warmstart we switch to our sparsified method to fully benefit from identification. This strategy is illustrated in Figure 3, with warmstarted algorithms on the right-hand-side vs. the non-warmstarted ones on the left-hand-side. We see a drastic improvement in terms of communication offered by the quick identification, for all versions of the sparsified method.

SM1 Inefficiency of uniform sparsification illustrated

In Section 3.4, we discussed the impact of the random selection in our sparsification technique. The uniform sparsification is proved to have a degraded convergence rate in terms of epochs; we illustrate in Figure SM1 that moreover it shows a degraded empirical performance in terms of exchanges.

SM2 Proof of Theorem 3.1

From the solution $x^{\star}$ of ($\mathsf{P}$) (unique from strong convexity), we define for each worker $i$, the local shift $x_{i}^{\star}=x^{\star}-\gamma\nabla f_{i}(x^{\star})$. From those, one can define $\bar{x}^{\star}=\sum_{i=1}^{M}\alpha_{i}x_{i}^{\star}$. First-order optimality conditions

[TABLE]

imply that

[TABLE]

which directly leads to $\mathbf{prox}_{\gamma r}(\bar{x}^{\star})=x^{\star}$ (see Chap. 16 of [5]).

We can now lay down the proof using these above-defined variables. For a time $k$ and a worker $i$, we have that $x_{i}^{k}=x_{i}^{k-d_{i}^{k}}$ depends on i) $x^{k-D_{i}^{k}}$ which is $\mathcal{F}^{k-D_{i}^{k}}$-measurable; and ii) $\mathbf{S}^{k-D_{i}^{k}}$ which is i.i.d. . First, we are going to control the term $\|x_{i}^{k}-x_{i}^{\star}\|^{2}$.

Let us define $\|x\|_{p}^{2}=\sum_{i=1}^{d}p_{i}x_{[i]}^{2}$ where $(p_{1},\dots,p_{d})$ is the vector of probabilities of Assumption 1. The conditional expectation can be developed as follows:

[TABLE]

Let us now bound both terms of this sum above using $p_{\max}=\max_{i}p_{i}$ and $p_{\min}=\min_{i}p_{i}$.

[TABLE]

We now use the $\mu$-strong convexity and $L$-smoothness of $f_{i}$ to write (see e.g. [9, Chap. 3.4.2]),

[TABLE]

Thus, for any $\gamma\in(0,2/(\mu+L)]$,

[TABLE]

where we used that $\|x^{k-D_{i}^{k}}-x^{\star}\|^{2}=\|\mathbf{prox}_{\gamma r}(\bar{x}^{k-D_{i}^{k}})-\mathbf{prox}_{\gamma r}(\bar{x}^{\star})\|^{2}\leq\|\bar{x}^{k-D_{i}^{k}}-\bar{x}^{\star}\|^{2}$ by definition and non-expansiveness of the proximity operator of $r$.

Taking full expectation on both sides, we get

[TABLE]

Then, using that $\bar{x}^{k-D_{i}^{k}}-\bar{x}^{\star}=\sum_{i=1}^{M}\alpha_{i}({x}_{i}^{k-D_{i}^{k}}-\bar{x}_{i}^{\star})$ and the convexity of $\|\cdot\|^{2}$, we get

[TABLE]

Let $c_{k}=\max_{i=1,\dots,M}\mathbb{E}\left\|x_{j}^{k}-x_{j}^{\star}\right\|^{2}$ and $\beta=\left(p_{\max}\left(1-\gamma\mu\right)^{2}+1-p_{\min}\right)$ (note that the assumptions imply that $\beta\leq 1$), then the above result implies that

[TABLE]

and using the definition of the sequence $(k_{m})$, we get

[TABLE]

Thus for all $k\geq k_{m}$, $c_{k}\leq\beta~{}\max_{\ell\in[k_{m-1},k_{m})}c_{\ell}$. This implies that the sequence $\widetilde{c}_{m}$ defined by $\widetilde{c}_{m}=\max_{\ell\in[k_{m},k_{m+1})}c_{\ell}$ has an exponential bound:

[TABLE]

Finally, we to use once again the non-expansivity of the proximity operator of $r$ and the definitions to get that for all $k\in[k_{m},k_{m+1}),$

[TABLE]

which concludes the proof.

SM3 Discussion on non-i.i.d. selections

Imagine that one wants to build an algorithm similar to Spy but where the update probabilities are equal to $1$ over an adaptively chosen set $\mathsf{A}^{k}$ (depending e.g. on the support of $x^{k}$) and $\pi<1$ for the others. As mentioned in the main text, this poses a problem since the selection is not i.i.d. anymore due to the adaptive selection.

This drawback can actually be fixed by “slowing down” the updates in the adaptive set $\mathsf{A}^{k}$ by a factor $\pi$. This gives Algorithm 1. A proof very similar to the one of Supplement SM2 enables to prove a $(1-\pi\gamma\mu(2-\gamma\mu))$ linear rate (the same as for uniform selection with probability $\pi$). Unfortunately, the artificial slowdown of the iterates is harmful in practice and causes this algorithm to be rather slow, which is why we chose not to consider this path.

SM4 Proofs of convergence for reconditioned algorithms

SM4.1. Basic lemmas for proofs of reconditioned algorithms

We state here two simple lemmas that we have not found as such in the literature and that are required in the proofs of Section 4.2. They use the fact that the unique minimum $x^{\star}$ of a strongly convex function $F$ is a fixed point of $\mathbf{prox}_{F/{\rho}}$ for any $\rho>0$; see [5, Prop. 12.28]).

Lemma SM4.1**.**

Let $F:\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$ be $\mu$-strongly convex lsc and $\rho>0$. Then, for any $x\in\mathbb{R}^{d}$,

[TABLE]

Proof.

The proof simply consists in developing norms as follows:

[TABLE]

then a reordering concludes the proof. The inequality uses the the fact that the resolvent of the $\mu/{\rho}$ strongly monotone operator $\partial F/{\rho}$ is $(1+\mu/{\rho})$-cocoercive; see [5], particularly Proposition 23.11. ∎

Lemma SM4.2**.**

Let $F:\mathbb{R}^{d}\to\mathbb{R}\cup\{+\infty\}$ be $\mu$-strongly convex lsc, $\rho>0$. Then, for any $x,x^{\prime}\in\mathbb{R}^{d}$ such that

[TABLE]

for some $\nu>0$, we have for any $\varepsilon\in(0,1)$

[TABLE]

Proof.

Using Young’s inequality, we get that for any $\varepsilon\in(0,1)$,

[TABLE]

where we used Lemma SM4.1. ∎

SM4.2. Proof of Theorem 4.1

Theorem 4.1 gives very similar results for the three inner stopping criteria. However, the proof techniques are quite different for each criteria.

Proof for criterion $\mathsf{C}_{1}$ (epoch budget)

We start by noticing that, at outer loop $\ell$, Spy solves the reconditioned problem (4.11) over which it has a contraction factor of $(1-\alpha_{\ell})$ with $1-\alpha_{\ell}:=1-\alpha+\pi-\pi_{\ell}$; see (4.6) and (3.3). This means that Spy initialized with $x_{\ell}$ verifies after $m$ epochs with the maximal777The initialization means that $x^{0}=x_{\ell}$ and $x_{i}^{0}=x_{\ell}-\gamma\nabla h_{i,\ell}(x_{\ell})$. The use of the maximal stepsize is required here to guarantee the optimal rate over the reconditioned problem, see (3.3). stepsize

[TABLE]

where $x_{\ell}^{\star}$ is the unique solution of ($\mathsf{R}_{\ell}$), and $x_{i,\ell}^{\star}=x_{\ell}^{\star}-\gamma\nabla h_{i,\ell}(x^{\star}_{\ell})$ are its local shifts.

We now apply Lemma SM4.2 with the following input: $x^{\prime}=x_{\ell+1}=x_{\ell}^{k_{\mathsf{M}_{\ell}}}$; $x=x_{\ell}$; $F=F$ (which is $\mu$ strongly convex); $x^{\star}=x^{\star}$ (minimizer of $F$); and $\nu=\left(1-\alpha_{\ell}\right)^{\mathsf{M}_{\ell}}$. We get for $\varepsilon=\frac{1}{\ell^{1+\delta}}\in(0,1)$

[TABLE]

Choosing $\mathsf{M}_{\ell}$ as per $\mathsf{C}_{1}$ guarantees that

[TABLE]

Thus, for any $\mu\geq 0$, we have

[TABLE]

Convergence. For any $\mu\geq 0$, (SM4.2) tells us that

[TABLE]

By Robbins-Siegmund theorem [53, Th. 1] (see also [20, 8, 10] for applications to optimization), we have that i) $(\left\|x_{\ell}-x^{\star}\right\|^{2})$ converges almost surely to a random variable with finite support; and ii) $\sum_{\ell=1}^{\infty}\|x_{\ell}-\mathbf{prox}_{F/{\rho}}(x_{\ell})\|^{2}<\infty$. This means that we can extract a subsequence $(x_{\ell^{n}})$ that converges almost surely to some $y$ which is necessarily a minimizer from ii). Using (SM4.3) again with $x^{\star}=y$, we see that $(x_{\ell})$ converges to $x^{\star}$ almost surely.

Rate. Now, if $\mu>0$, we get by dropping the last term in (SM4.2) and successively taking expectations that

[TABLE]

Proof for criterion $\mathsf{C}_{2}$ (absolute accuracy)

We apply Lemma SM4.2 with the following input: $x^{\prime}\!=x_{\ell+1}$; $x=x_{\ell}$; $F=F$ (which is $\mu$ strongly convex); $x^{\star}=x^{\star}$ (minimizer of $F$); and $\nu=(1-\delta)\rho/((2\mu+\rho)\ell^{1+\delta})$ (noting that the condition on $x_{\ell+1}$ is almost sure). We get for $\varepsilon=\frac{1}{\ell^{1+\delta}}\in(0,1)$

[TABLE]

This directly gives the rate of convergence when $\mu>0$. When $\mu=0$, the inequality can be simplified to

[TABLE]

In this case, the same arguments as for criterion $\mathsf{C}_{1}$ enable to get almost sure convergence.

Proof for criterion $\mathsf{C}_{3}$ (relative accuracy)

Denoting $\beta:=\sqrt{\rho/(2\mu+\rho)}\in(0,1]$, the stopping criterion $\mathsf{C}_{3}$ writes

[TABLE]

Convergence. The condition (SM4.5) matches condition (B) of [54, Th. 2]. We also have clearly $\sum_{\ell}\varepsilon_{\ell}<+\infty$ and the regularity assumption of the operator is verified, by our extra assumption and [54, Prop. 7]. Thus [54, Th. 2] directly gives us that $(x_{\ell})$ converges to a minimizer of $F$, that we denote by $x^{\star}$.

Rate. When $F$ is $\mu$-strongly convex, we can furthermore develop:

[TABLE]

where the first inequality used both condition $\mathsf{C}_{3}$ and Lemma SM4.1. This implies that

[TABLE]

Finally, denoting $d_{\ell}:=\ell^{1+\delta}/(\ell^{1+\delta}-1)>1$, we have888Note that $\varepsilon_{\ell}^{2}=\frac{\rho}{4(2\mu+\rho)\ell^{2+2\delta}}\leq\left(\frac{\beta}{2}\right)^{2}\frac{(d_{\ell}-1)^{2}}{d_{\ell}^{2}}\leq\beta^{2}\frac{(d_{\ell}-1)^{2}}{(1+\beta d_{\ell})^{2}}$.

[TABLE]

This yields

[TABLE]

which gives the result.

SM5 Acceleration à la Catalyst

Catalyst is a popular accelerated inexact proximal point meta-algorithm for solving machine learning objectives [31, 32]. It consists in adding an acceleration step for each outer iteration to reach a faster rate. The parameters of Catalyst are i) the choice of the inertial sequence $(\beta_{\ell})$; and ii) the choice of the stopping criteria that comes in the same flavor as for Reconditioned–Spy. Applied to our context, this gives Algorithm 1 where the main difference with Reconditioned–Spy is the final line with the acceleration.

Under the parameters given in [32], the fast outer rate is stated in the following result which shows in the square root compared to Theorem 4.1 for the direct proximal reconditioning.

Theorem SM5.1**.**

Let the functions $(f_{i})$ be $\mu$-strongly convex ($\mu\geq 0$) and $L$-smooth. Let $r$ be convex lsc. Then, Catalyst–Spy on $((\alpha_{i}),(f_{i}),r)$ with stopping criterion $\mathsf{C}^{\prime}_{2}$ or $\mathsf{C}^{\prime}_{3}$ converges in terms of suboptimality at rate $1/\ell^{2}$.

Furthermore, if $\mu>0$, then

[TABLE]

Proof.

The result directly follows from [32], in particular Propositions 5,6,8, and 9. ∎

Note that $\mathsf{C}_{1}$ is not supported in this result, since it is not covered by theory in [32]; however, a criterion close to $\mathsf{C}_{1}$ is used in the computational experiments of [32].

As for the other two criteria, $\mathsf{C}^{\prime}_{2}$ and $\mathsf{C}^{\prime}_{3}$ imply rules similar to $\mathsf{C}_{2}$ and $\mathsf{C}_{3}$ since

[TABLE]

by using the $(\mu+\rho)$-strong convexity of the inner problem. We note that $\mathsf{C}^{\prime}_{2}$ is then much more stringent that $\mathsf{C}_{2}$. For $\mathsf{C}^{\prime}_{3}$, the conditions are similar to $\mathsf{C}_{3}$ in the non-strongly convex case, however they appear to be better in the strongly convex case. This is an artefact due to our use of one criterion both strongly convex and non-strongly convex cases. Indeed, in the strongly convex case, one can see from (SM4.6) that a fixed choice of $\varepsilon_{\ell}$ is possible, leading to a similar strategy and rate.

SM6 Proofs related to identification

SM6.1. Proof of Theorem 5.1

We notice that DAve-PG, as many proximal algorithms (see e.g. [21, Sec. 4]), has the following convergence property

[TABLE]

from Theorem 2.1. This property is the key to get the identification result, as we will develop at the end of the proof. We start by proving that Reconditioned–Spy (as well as Catalyst–Spy) shares a similar convergence property.

Both Reconditioned–Spy and Catalyst–Spy are of the form

[TABLE]

(with $\beta_{\ell}\equiv 0$ for Reconditioned–Spy and $\beta_{\ell}\in(0,1)$ for Catalyst–Spy); and verify for all $\ell,k$

[TABLE]

for some $C,C^{\prime}>0$ and $\rho\in(0,1)$, where

[TABLE]

As in Supplement SM2, we also consider $\bar{x}^{\star}=x^{\star}-\gamma\sum_{i=1}^{M}\alpha_{i}\nabla f_{i}(x^{\star})$ given by (SM2.1). Then, we have

[TABLE]

for some $D,D^{\prime}>0$. For any $\ell,k$, we then have

[TABLE]

Hence, by Markov’s inequality and Borel-Cantelli’s lemma, $\overline{x}_{\ell}^{k}\to\overline{x}^{\star}$ almost surely. As a direct result, we get for our two random algorithms, the same convergence as (SM6.1) for DAve-PG,

[TABLE]

This convergence implies identification of optimal support (see e.g. the recent survey [21, Cor. 1]). Recalling Assumption 2i, there exists $\varepsilon>0$ such that

[TABLE]

for all $\mathbf{e}\in\mathcal{B}(0,\gamma\varepsilon)$. Hence, since $\bar{x}_{\ell}^{k}\to_{\ell\to\infty}\bar{x}^{\star}$ almost surely, $\bar{x}_{\ell}^{k}$ will belong to the ball of radius $\gamma\varepsilon$ centered on $\bar{x}^{\star}$ in a finite number of outer iterations, say $\Lambda<\infty$, with probability one. The equation above then directly implies that for all $\ell\geq\Lambda$, $\mathrm{supp}(x^{\star})=\mathrm{supp}(\mathbf{prox}_{\gamma r}(\bar{x}_{\ell}^{k}))=\mathrm{supp}(x_{\ell}^{k})$.

Finally, it suffices to notice that $x_{\ell+1}=x_{\ell}^{k}$ for some $k$ to conclude the proof.

SM6.2. Proof of Theorem 5.2

From Theorem 5.1 we know that identification takes place i.e. that we have $\mathrm{null}(x^{k}_{\ell})=\mathrm{null}(x^{\star}):=n^{\star}$ for all $k,\ell$ ($\ell\geq\Lambda$). In this case,

[TABLE]

where we denote by $\overline{n^{\star}}$ the complementary of $n^{\star}$. In addition, if $x=\mathbf{prox}_{\gamma r}(\bar{x})$ is such that $\mathrm{null}(x)=n^{\star}$, then $x=[x]_{\overline{n^{\star}}}=\mathbf{prox}_{\gamma r}([\bar{x}]_{\overline{n^{\star}}})$ by separability of $r$. Then, recalling that all coordinates in $\mathrm{supp}(x_{\ell}^{k})=\overline{n^{\star}}$ are updated with the choice of (4.10), we have for all $k,\ell\geq\Lambda$

[TABLE]

This exactly coincides with a non-sparsified update on the restriction of $f_{i}$ to the subspace of vectors with null coordinates in ${\overline{n^{\star}}}$. More specifically, let $S^{\star}=\{x\in\mathbb{R}^{d}:\mathrm{null}(x)=n^{\star}\}$ be the subspace of vectors with null coordinates in ${\overline{n^{\star}}}$, and $f_{i|{\overline{n^{\star}}}}$ be the restriction of $f_{i}$ to $S^{\star}$. Then the above iteration coincides with non-sparsified update on $((\alpha_{i}),(f_{i|{\overline{n^{\star}}}}),r)$. In other words, after identification, Spy is no longer random and has the same iterates as DAve-PG on $S^{\star}$ (while in $S^{\star\perp}$, the algorithm has converged to [math]). Theorem 2.1 therefore guarantees that Spy benefits from a $(1-\gamma(\mu+\rho))^{2}$ rate in terms of epochs (since $(\mu+\rho)$ is the modulus of strong convexity of each $f_{i|{\overline{n^{\star}}}}$).

SM6.3. Proof of Theorem 5.3

Following the choice of (4.10),

[TABLE]

for both $\mathsf{C}_{2}$ and $\mathsf{C}_{3}$ from Theorem 5.2 (see also [32, Lemma 11]). Finally,

[TABLE]

from Theorem 4.1. Put together, this gives the following complexity:

[TABLE]

Since $c\approx s^{\star}\ll d$, we have by definition (in (4.5)) $\rho\gg L\geq\mu$ (we will use $\mu+L\leq 2(\sqrt{2}-1)\rho$ in the $\lesssim$ step of the upcoming equation). Taking the maximal $\gamma=2/(\mu+L+2\rho)$, we have

[TABLE]

SM7 Numerical illustration of the stopping tests

In Algorithm 1, the stopping criteria $\mathsf{C}_{1}$ with prescribed number of inner loops is simple and natural. In contrast, enforcing criteria $\mathsf{C}_{3}$ needs a full gradient evaluation (see e.g. [54]), which requires some partial synchronization and then breaks the asynchronous nature of the method.

In our numerical experiments, we use a simplification of the stopping test $\mathsf{C}_{1}$ with only one epoch, following the principle of “1 pass over the data = 1 restart” used in Catalyst [32]. We illustrate here that this simple stopping rule gives similar empirical convergence as $\mathsf{C}_{3}$ without its computational issues.