Sparsity-Exploiting Distributed Projections onto a Simplex

Projection onto a simplex can be leveraged as a subroutine to determine projections onto more complex polyhedra. Such projections arise in numerous settings such as: image processing, e.g. labeling (Lellmann et al. 2009), or multispectral unmixing (Bioucas-Dias et al. 2012); portfolio optimization (Brodie et al. 2009); and machine learning (Blondel et al. 2014). As a particular example, projection onto a simplex can be used to project onto the parity polytope (Wasson et al. 2019):

where $\mathbb{PP}_{n}$ is a $n$-dimensional parity polytope:

Projection onto the parity polytope arises in linear programming (LP) decoding (Liu and Draper 2016, Barman et al. 2013, Zhang and Siegel 2013, Zhang et al. 2013, Wei et al. 2015), which is used for signal processing.

Another example is projection onto a $\ell_{1}$ ball:

Duchi et al. (2008) demonstrate that the solution to this problem can be easily recovered from projection onto a simplex. Furthermore, projection onto a $\ell_{1}$ ball can, in turn, be used as a subroutine in gradient-projection methods (see e.g. van den Berg (2020)) for a variety of machine learning problems that use $\ell_{1}$ penalty, such as: Lasso (Tibshirani 1996); basis-pursuit denoising (Chen et al. 1998, van den Berg and Friedlander 2009, van den Berg 2020); sparse representation in dictionaries (Donoho and Elad 2003); variable selection (Tibshirani 1997); and classification (Barlaud et al. 2017).

Finally, we note that methods for projection onto the scaled standard simplex and $\ell_{1}$ ball can be extended to projection onto the weighted simplex and weighted $\ell_{1}$ ball (Perez et al. 2020a) (see Section B.2). Projection onto the weighted simplex can, in turn, be used to solve the continuous quadratic knapsack problem (Robinson et al. 1992). Moreover, $\ell_{p}$ regularization can be handled by iteratively solving weighted $\ell_{1}$ projections (Candès et al. 2008, Chartrand and Yin 2008, Chen and Zhou 2014).

This paper presents a method to decompose the projection problem and distribute work across (up to $n$) processors. The key insight to our approach is captured by Proposition 3.4: the projection of any subvector of $d$ onto a simplex (in the corresponding space with scale factor $b$) will have zero-valued entries only if the full-dimension projection has corresponding zero-valued entries. The method can be interpreted as a sparsity-exploiting distributed preprocessing method, and thus it can be adapted to parallelize a broad range of serial projection algorithms. We furthermore provide extensive theoretical and empirical analyses of several such adaptations. We also fill in gaps in the literature on serial algorithm complexity. Our computational results demonstrate significant speedups from our distributed method compared to the state-of-the-art over a wide range of large-scale problems involving both real-world and simulated data.

Our paper contributes to the limited literature on parallel computation for large-scale projection onto a simplex. Most of the algorithms for projection onto a simplex are for serial computing. Indeed, to our knowledge, there is only one published parallel method, and one distributed method for projection problem (1). Wasson et al. (2019) parallelize a basic sort and scan (specifically prefix sum) approach–a method that we use as a benchmark in our experiments. We also develop a modest but practically significant enhancement to their approach. Iutzeler and Condat (2018) propose a gossip-based distributed ADMM algorithm for projection onto a Simplex. In this gossip-based setup, one entry of $d$ is given to each agent (e.g. processor), and communication is restricted according to a particular network topology. This differs fundamentally from our approach both in context and intended use, as we aim to solve large-scale problems and moreover our method can accommodate any number of processors up to $n$.

The remainder of the paper is organized as follows. Section 2 describes serial algorithms from the literature and develops new complexity results to fill in gaps in the literature. Section 3 develops parallel analogues of the aforementioned algorithms using our novel distributed method. Section 4 extends these parallelized algorithms to various applications of projection onto a simplex. Section 5 describes computational experiments. Section 6 concludes. Note that all appendices, mathematical proofs, as well as our code and data can be found in the online supplement.

KKT conditions characterize the unique optimal solution $v^{*}$ to problem (1):

Hence, (1) can be reduced to a univariate search problem for the optimal pivot $\tau$. Note that the nonzero (positive) entries of $v^{*}$ correspond to entries where $d_{i}>\tau$. So for a given value $t\in\mathbb{R}$ and the index set $\mathcal{I}:=\{1,...,n\}$, we denote the active index set

as the set containing all indices of active elements where $d_{i}>t$. Now consider the following function, which will be used to provide an alternate characterization of $\tau$:

The sign of $f$ only changes once, and $\tau$ is its unique root. These results can be leveraged to develop search algorithms for $\tau$, which are presented next. This corollary and the use of $f$ is our own contribution, as we have found it a convenient organizing principle for the sake of exposition We note, however, that the root finding framework has been in use in the more general constrained quadratic programming literature (see (Cominetti et al. 2014, Equation 5) and (Dai et al. 2006, Section 2, Paragraph 2)).

Observe that only the greatest $|\mathcal{I}_{\tau}|$ terms of $d$ are indexed in $\mathcal{I}_{\tau}$. Now suppose we sort $d$ in non-increasing order:

We can sequentially test these values in descending order, $f(d_{\pi_{1}}),f(d_{\pi_{2}}),$ etc. to determine $|\mathcal{I}_{\tau}|$. In particular, from Corollary 2.2 we know there exists some $\kappa:=|\mathcal{I}_{\tau}|$ such that $f(d_{\pi_{\kappa}})<0\leq f(d_{\pi_{\kappa+1}})$. Thus the projection must have $\kappa$ active elements, and since $f(\tau)=0$, we have $\tau=\frac{\sum_{i=1}^{\kappa}d_{\pi_{i}}-b}{\kappa}$. We also note that, rather than recalculating $f$ at each iteration, one can keep a running cumulative/prefix sum or scan of $\sum_{i=1}^{j}d_{\pi_{i}}$ as $j$ increments.

The bottleneck is sorting, as all other operations are linear time; for instance, QuickSort executes the sort with average complexity $O(n\log n)$ and worst-case $O(n^{2})$, while MergeSort has worst-case $O(n\log n)$ (see, e.g. (Bentley and McIlroy 1993)). Moreover, non-comparison sorting methods can achieve $O(n)$ (see, e.g. (Mahmoud 2000)), albeit typically with a high constant factor as well as dependence on the bit-size of $d$.

Sort and Scan begins by sorting all elements of $d$, but only the greatest $|\mathcal{I}_{\tau}|$ terms are actually needed to calculate $\tau$. Pivot and Partition, proposed by Duchi et al. (2008), can be interpreted as a hybrid sort-and-scan that attempts to avoid sorting all elements. We present as Algorithm 2 a variant of this method approach given by Condat (2016).

The algorithm selects a pivot $p\in[\min_{i}\{d_{i}\},\max_{i}\{d_{i}\}]$, which is intended as a candidate value for $\tau$; the corresponding value $f(p)$ is calculated. From Corollary 2.2, if $f(p)>0$, then $p<\tau$ and so $\mathcal{I}_{p}\supset\mathcal{I}_{\tau}$; consequently, a new pivot is chosen in the (tighter) interval $[\min_{i\in\mathcal{I}_{p}}\{d_{i}\},\max_{i\in\mathcal{I}_{p}}\{d_{i}\}]$, which is known to contain $\tau$. Otherwise, if $f(p)\leq 0$, then $p\geq\tau$, and so we can find a new pivot $p\in[\min_{i\in\bar{\mathcal{I}}_{p}}\{d_{i}\},\max_{i\in\bar{\mathcal{I}}_{p}}\{d_{i}\}]$, where $\bar{\mathcal{I}}_{p}:=\{1,...,n\}\setminus\mathcal{I}_{p}$ is the complement set. Repeatedly selecting new pivots and creating partitions in this manner results in a binary search to determine the correct active set $\mathcal{I}_{\tau}$, and consequently $\tau$.

Several strategies have been proposed for selecting a pivot within a given interval. Duchi et al. (2008) choose a random value in the interval, while Kiwiel (2008) uses the median value. The classical approach of Michelot (1986) can be interpreted as a special case that sets the initial pivot as $p^{(1)}=(\sum_{i\in\mathcal{I}}d_{i}-b)/|\mathcal{I}|$, and subsequently $p^{(i+1)}=f(p^{(i)})+p^{(i)}$. This ensures that $p^{(i)}\leq\tau$ which avoids extraneous re-evaluation of sums in the if condition. Note that $p^{(i)}$ generates an increasing sequence converging to $\tau$ (Condat 2016, Page 579, Paragraph 2). Michelot’s algorithm is presented separately as Algorithm 3.

Condat (2016) provided worst-case runtimes for each of the aforementioned pivot rules, as well as average case complexity (over the uniform distribution) for the random pivot rule (see Table 2). We fill in the gaps here and establish $O(n)$ runtimes for the median rule as well as Michelot’s method. We note that the median pivot method is a linear-time algorithm, but relies on a median-of-medians subroutine (Blum et al. 1973), which has a high constant factor. For Michelot’s method, we assume uniformly distributed inputs, $d_{1},\dots,d_{n}$ are $\mathrm{i.i.d}\sim U[l,u]$, and we have

The same argument holds for the median pivot rule, as half the elements are guaranteed to be removed each iteration, and the operations per iteration are within a constant factor of Michelot’s; we omit a formal proof of its $O(n)$ average runtime for brevity.

Condat’s method (Condat 2016), presented as Algorithm 5, can be seen as a modification of Michelot’s method in two ways. First, Condat replaces the initial scan with a Filter to find an initial pivot, presented as Algorithm 4. Lemma 2 (see Appendix. A) establishes that the Filter provides a greater (or equal) initial starting pivot compared to Michelot’s initialization; furthermore, since Michelot approaches $\tau$ from below, this results in fewer iterations (see proof of Proposition 2.4). Second, Condat’s method dynamically updates the pivot value whenever an inactive entry is removed from $\mathcal{I}_{p}$, whereas Michelot’s method updates the pivot every iteration by summing over all entries.

Condat (2016) supplies a worst-case complexity of $O(n^{2})$. We supplement this with average-case analysis under uniformly distributed inputs, e.g.$d_{1},\dots,d_{n}$ are $\mathrm{i.i.d}\sim U[l,u]$.

Table 1 shows that all presented algorithms attain $O(n)$ performance on average given uniformly i.i.d. inputs. The methods are ordered by publication date, starting from the oldest result. As described in Section 2.2, Sort and Scan can be implemented with non-comparison sorting to achieve $O(n)$ worst-case performance. However, as with the linear-time median pivot rule, there are tradeoffs: increased memory, overhead, dependence on factors such as input bit-size, etc.

Both sorting and scanning are (separately) well-studied in parallel algorithm design, so the Sort and Scan idea lends itself to a natural decomposition for parallelism (discussed in Section 3.1). The other methods integrate sorting and scanning in each iterate, and it is no longer clear how best to exploit parallelism directly. We develop in the next section a distributed preprocessing scheme that works around this issue in the case of sparse projections. Note that the table includes the Bucket Method; details on the algorithm are provided in Appendix B.1.

Wasson et al. (2019) parallelize Sort and Scan in a natural way: first applying a parallel merge sort (see e.g. (Cormen et al. 2009, p. 797)) and then a parallel scan (Ladner and Fischer 1980) on the input vector. However, their scan calculates $\sum_{i=1}^{j}d_{\pi_{i}}$ for all $j\in\mathcal{I}$, but only $\sum_{i=1}^{\kappa}d_{\pi_{\kappa}}$ is needed to calculate $\tau$. We modify the algorithm accordingly, presented as Algorithm 6: checks are added (lines 7 and 14) in the for-loops to allow for possible early termination of scans. As we are adding constant operations per loop, Algorithm 6 has the same complexity as the original Parallel Sort and Scan. We combine this with parallel mergesort in Algorithm 7 and empirically benchmark this method with the original (parallel) version in Section 5.

Our main idea is motivated by the following two theorems that establish that projections with i.i.d. inputs and fixed right-hand-side $b$ become increasingly sparse as the problem size $n$ increases.

Theorem 3.1 establishes that, for i.i.d. uniformly distributed inputs, the projection has $O(\sqrt{n})$ active entries in expectation and thus has considerable sparsity as $n$ grows; we also show this in the computational experiments of Appendix E.1.

Theorem 3.2 establishes arbitrarily sparse projections over arbitrary i.i.d. distributions given fixed $b$. Note that if $b$ is sufficiently large with respect to $n$ (rather than fixed), then the resulting projection could be too dense to attain the theorem result. However, sparsity can be assured provided $b$ does not grow too quickly with respect to $n$, namely:

We apply Theorem 3.2 to example distributions in Appendix C, and test the bounds empirically in Appendix E.1.

Proposition 3.4 tells us that if we project a subvector of some length $m\leq n$ onto the same $b$-scaled simplex in the corresponding $\mathbb{R}^{m}$ space, the zero entries in the projected subvector must also be zero entries in the projected full vector.

Our idea is to partition and distribute the vector $d$ across cores (broadcast); have each core find the projection of its subvector (local projection); and combine the nonzero entries from all local projections to form a vector $\hat{v}$ (reduce), and apply a final (global) projection to $\hat{v}$. The method is outlined in Figure 1. Provided the projection $v^{*}$ is sufficiently sparse, which (for instance) we have established is the case for i.i.d. distributed large-scale problems, we can expect $\hat{v}$ to have been far less than $n$ entries. We demonstrate the practical advantages of this procedure with various computational experiments in Section 5.

The distributed method outlined in Figure 1 can be applied directly to Pivot and Partition, as described in Algorithm 8. Note that, as presented, $v^{*}$ is a sparse vector: entries not processed in the final Pivot and Project iteration are set to zero (recall Proposition 3.4).

We assume uniformly distributed inputs, $d_{1},\dots,d_{n}$ are $\mathrm{i.i.d}\sim U[l,u]$, and we have

In the worst case we may assume the distributed projections are ineffective, $\mbox{dim}(\hat{v})\in O(n)$, and so the final projection is bounded above by $O(n^{2})$ with random pivots and Michelot’s method, and $O(n)$ with the median pivot rule.

We could apply the distributed sparsity idea as a preprocessing step for Condat’s method. However, due to Proposition 3.6 (and confirmed via computational experiments) we have found that Filter tends to discard many non-active elements. Therefore, we propose instead to apply our distributed method to parallelize the Filter itself. Our Distributed Filter is presented in Algorithm 9: we partition $d$ and broadcast it to the cores, and in each core we apply (serial) Filter on its subvector. Condat’s method with the distributed Filter is presented as Algorithm 10.

We assume uniformly distributed inputs, $d_{1},\dots,d_{n}$ are $\mathrm{i.i.d}\sim U[l,u]$, and we have

Under the same assumption (uniformly distributed inputs) to Proposition 3.6, we have

In the worst-case we can assume Distributed Filter is ineffective ($|\mathcal{I}_{p}|\in O(n)$), and so the complexity of Parallel Condat’s method is $O(n^{2})$, same as the serial method.

Complexity results for parallel algorithms developed throughout this section, as well as their serial counterparts, are presented in Table 2. Parallelized Sort and Scan has a dependence on $\frac{1}{k}$; indeed, sorting (and scanning) are very well-studied problems from the perspective of parallel computing. Parallel (Bitonic) Merge Sort has an average-case (and worst case) complexity in $O((n\log n)/k)$ (Greb and Zachmann 2006), and Parallel Scan has an average-case (and worst case) complexity in $\Theta(n/k+\log k)$ (Wang et al. 2020); thus, running parallel sort followed by parallel scan is $O(\frac{n}{k}\log n)$. Now, Michelot’s and Condat’s serial methods are observed to have favorable practical performance in our computational experiments; this is expected as these more modern approaches were explicitly developed to gain practical advantages in e.g. the constant runtime factor. Conversely, our distributed method does not improve upon the worst-case complexity of Michelot’s and Condat’s methods, but is able to attain a $\frac{1}{k}$ factor for average complexity for $n\gg k$, which is the case on large-scale instances, i.e. for all practical purposes. The average case analyses were conducted under the admittedly limited (typical) assumption of uniform i.i.d. entries, but our computational experiments over other distributions and real-world data confirm favorable practical speedups from our parallel algorithms.

Consider projection onto an $\ell_{1}$ ball:

where $\mathcal{B}_{b}$ is given by Equation (3). Duchi et al. (2008) show Problem (6) is linear-time reducible to projection onto simplex (see (Duchi et al. 2008, Section 4)). Hence, any parallel method for projection onto a simplex can be applied to Problem (6).

As mentioned in Section 1, Problem (6) can itself be used as a subroutine in solving the Lasso problem, via (e.g.) Projected Gradient Descent (PGD) (see e.g. (Boyd and Vandenberghe 2004, Exercise 10.2)). To handle large-scale datasets, we instead use the mini-batch gradient descent method (Zhang et al. 2023, Sec. 12.5) in PGD in Sec 5.

Leveraging the solution to Problem (1), Wasson et al. (Wasson et al. 2019, Algorithm 2) develop a method to project a vector onto the centered parity polytope, $\mathbb{PP}_{n}-\frac{1}{2}$ (recall Problem 2); we present a slightly modified version as Algorithm 2 in Appendix D. The modification is on line $11$, where we determine whether a simplex projection is required to avoid unnecessary operations; the original method executes line $14$ before line $11$.

In this subsection, we compare runtime results for serial methods and their parallel versions with two measures: absolute speedup and relative speedup. Absolute speedup is the parallel method’s runtime vs the fastest serial method’s runtime, e.g. serial Condat time divided by parallel Sort and Scan time; relative speedup is the parallel method’s runtime vs its serial equivalent’s runtime, e.g. serial Sort and Scan time divided by parallel Sort and Scan time. We test parallel implementations using $8,16,24,...,80$ cores. Note that we have verified that each parallel algorithm run on a single core is slower than its serial equivalent (see Appendix E.4).