Fast Differentiable Sorting and Ranking

Introduction

Modern deep learning architectures are built by composing parameterized functional blocks (including loops and conditionals) and are trained end-to-end using gradient backpropagation. This has motivated the term differentiable programming, recently popularized, among others, by LeCun (2018). Despite great empirical successes, many operations commonly used in computer programming remain poorly differentiable or downright pathological, limiting the set of architectures for which a gradient can be computed.

We focus in this paper on two such operations: sorting and ranking. Sorting returns the given input vector with its values re-arranged in monotonic order. It plays a key role to handle outliers in robust statistics, as in least-quantile (Rousseeuw, 1984) or trimmed (Rousseeuw & Leroy, 2005) regression. As a piecewise linear function, however, the sorted vector contains many kinks where it is non-differentiable. In addition, when used in composition with other functions, sorting often induces non-convexity, thus rendering model parameter optimization difficult.

In this paper, we propose the first differentiable sorting and ranking operators with $O(n\log n)$ time and $O(n)$ memory complexity. Our proposals enjoy exact computation and differentiation (i.e., they do not involve differentiating through the iterates of an approximate algorithm). We achieve this feat by casting differentiable sorting and ranking as projections onto the permutahedron, the convex hull of all permutations, and using a reduction to isotonic optimization. While the permutahedron had been used for learning before (Yasutake et al., 2011; Ailon et al., 2016; Blondel, 2019), it had not been used to define fast differentiable operators. The rest of the paper is organized as follows.

We review the necessary background (§2) and show how to cast sorting and ranking as linear programs over the permutahedron, the convex hull of all permutations (§3).

We introduce regularization in these linear programs, which turns them into projections onto the permutahedron and allows us to define differentiable sorting and ranking operators. We analyze the properties of these operators, such as their asymptotic behavior (§4).

Using a reduction to isotonic optimization, we achieve $O(n\log n)$ computation and $O(n)$ differentiation of our operators, a key technical contribution of this paper (§5).

We show that our approach is an order of magnitude faster than existing approaches and showcase two novel applications: differentiable Spearman’s rank coefficient and soft least trimmed squares (§6).

Preliminaries

We define the argsort of $\bm{\theta}$ as the indices sorting $\bm{\theta}$, i.e.,

where $\theta_{\sigma_{1}(\bm{\theta})}\geq\dots\geq\theta_{\sigma_{n}(\bm{\theta})}$. If some of the coordinates of $\bm{\theta}$ are equal, we break ties arbitrarily. We define the sort of $\bm{\theta}$ as the values of $\bm{\theta}$ in descending order, i.e.,

We define the rank of $\bm{\theta}$ as the function evaluating at coordinate $j$ to the position of $\theta_{j}$ in the descending sort (smaller rank $r_{j}(\bm{\theta})$ means that $\theta_{j}$ has higher value). It is formally equal to the argsort’s inverse permutation, i.e.,

For instance, if $\theta_{3}\geq\theta_{1}\geq\theta_{2}$, then $\sigma(\bm{\theta})=(3,1,2)$, $s(\bm{\theta})=(\theta_{3},\theta_{1},\theta_{2})$ and $r(\bm{\theta})=(2,3,1)$. All three operations can be computed in $O(n\log n)$ time. Note that throughout this paper, we use descending order for convenience. The ascending order counterparts are easily obtained by $\sigma(-\bm{\theta})$, $-s(-\bm{\theta})$ and $r(-\bm{\theta})$, respectively.

Sorting and ranking as linear programs

We show in this section how to cast sorting and ranking operations as linear programs over the permutahedron. To that end, we first formulate the argsort and ranking operations as optimization problems over the set of permutations $\Sigma$. {lemma}Discrete optimization formulations

A well-known object in combinatorics (Bowman, 1972; Ziegler, 2012), the permutahedron of ${\bm{w}}$ is a convex polytope, whose vertices correspond to permutations of ${\bm{w}}$. It is illustrated in Figure 1. In particular, when ${\bm{w}}=\bm{\rho}$, $\mathcal{P}({\bm{w}})=\operatorname*{conv}(\Sigma)$. With this defined, we can now derive linear programming formulations of sort and ranks. {proposition}Linear programming formulations

A proof is provided in §B.2. The key idea is to perform a change of variable to “absorb” the permutation in (4) and (5) into a permutahedron. From the fundamental theorem of linear programming (Dantzig et al., 1955, Theorem 6), an optimal solution of a linear program is almost surely achieved at a vertex of the convex polytope, a permutation in the case of the permutahedron. Interestingly, $\bm{\theta}$ appears in the constraints and $\bm{\rho}$ appears in the objective for sorting, while this is the opposite for ranking.

For $s(\bm{\theta})$, the fact that $\bm{\theta}$ appears in the linear program constraints makes $s(\bm{\theta})$ piecewise linear and thus differentiable almost everywhere. When $\sigma(\bm{\theta})$ is unique at $\bm{\theta}$, $s(\bm{\theta})=\bm{\theta}_{\sigma(\bm{\theta})}$ is differentiable at $\bm{\theta}$ and its Jacobian is the permutation matrix associated with $\sigma(\bm{\theta})$. When $\sigma(\bm{\theta})$ is not unique, we can choose any matrix in Clarke’s generalized Jacobian, i.e., any convex combination of the permutation matrices associated with $\sigma(\bm{\theta})$.

Lack of useful Jacobian of ranking.

On the other hand, for $r(\bm{\theta})$, since $\bm{\theta}$ appears in the objective, a small perturbation to $\bm{\theta}$ may cause the solution of the linear program to jump to another permutation of $\bm{\rho}$. This makes $r(\bm{\theta})$ a discontinuous, piecewise constant function. This means that $r(\bm{\theta})$ has null or undefined partial derivatives, preventing its use within a neural network trained with backpropagation.

Differentiable sorting and ranking

As we have already motivated, our primary goal is the design of efficiently computable approximations to the sorting and ranking operators, that would smoothen the numerous kinks of the former, and provide useful derivatives for the latter. We achieve this by introducing strongly convex regularization in our linear programming formulations. This turns them into efficiently computable projection operators, which are differentiable and amenable to formal analysis.

i.e., the Euclidean projection of $\bm{z}$ onto $\mathcal{P}({\bm{w}})$. We also consider entropic regularization $E(\bm{\mu})\coloneqq\langle\bm{\mu},\log\bm{\mu}-\mathbf{1}\rangle$, popularized in the optimal transport literature (Cuturi, 2013; Peyré & Cuturi, 2017). Subtly, we define

More generally, we can use any strongly convex regularization $\Psi$ under mild conditions. For concreteness, we focus our exposition in the main text on $\Psi\in\{Q,E\}$. We state all our propositions for these two cases and postpone a more general treatment to the appendix.

Soft operators.

We now build upon these projections to define soft sorting and ranking operators. To control the regularization strength, we introduce a parameter $\varepsilon>0$ which we multiply $\Psi$ by (equivalently, divide $\bm{z}$ by).

For sorting, we choose $(\bm{z},{\bm{w}})=(\bm{\rho},\bm{\theta})$ and therefore define the $\Psi$-regularized soft sort as

For ranking, we choose $(\bm{z},{\bm{w}})=(-\bm{\theta},\bm{\rho})$ and therefore define the $\Psi$-regularized soft rank as

We illustrate the behavior of both of these soft operations as we vary $\varepsilon$ in Figures 2 and 3. As for the hard versions, the ascending-order soft sorting and ranking are obtained by negating the input as $-s_{\varepsilon\Psi}(-\bm{\theta})$ and $r_{\varepsilon\Psi}(-\bm{\theta})$, respectively.

Properties.

We can further characterize these approximations. Namely, as we now formalize, they are differentiable a.e., and not only converge to the their “hard” counterparts, but also satisfy some of their properties for all $\varepsilon$. {proposition}Properties of $s_{\varepsilon\Psi}(\bm{\theta})$ and $r_{\varepsilon\Psi}(\bm{\theta})$

Differentiability. For all $\varepsilon>0$, $s_{\varepsilon\Psi}(\bm{\theta})$ and $r_{\varepsilon\Psi}(\bm{\theta})$ are differentiable (a.e.) w.r.t. $\bm{\theta}$.

where $f_{Q}(\bm{u})\coloneqq\text{mean}(\bm{u})$, $f_{E}(\bm{u})\coloneqq\log f_{Q}(\bm{u})$.

The last property describes the behavior as $\varepsilon\to 0$ and $\varepsilon\to\infty$. Together with the proof of Section 4, we include in §B.3 a slightly stronger result. Namely, we derive an explicit value of $\varepsilon$ below which our operators are exactly equal to their hard counterpart, and a value of $\varepsilon$ above which our operators can be computed in closed form.

Convexification effect.

Proposition 4 shows that $[s_{\varepsilon\Psi}(\bm{\theta})]_{i}$ and $[r_{\varepsilon\Psi}(\bm{\theta})]_{i}$ for all $i\in[n]$ converge to convex functions of $\bm{\theta}$ as $\varepsilon\to\infty$. This suggests that larger $\varepsilon$ make the objective function increasingly easy to optimize (at the cost of departing from “hard” sorting or ranking). This behavior is also visible in Figure 3, where $[s_{\varepsilon Q}(\bm{\theta})]_{2}$ converges towards the mean $f_{Q}$, depicted by a straight line.

On tuning ε𝜀\varepsilon (or not).

The parameter $\varepsilon>0$ controls the trade-off between approximation of the original operator and “smoothness”. When the model $g(\bm{x})$ producing the scores or “logits” $\bm{\theta}$ to be sorted/ranked is a homogeneous function, from (12) and (13), $\varepsilon$ can be absorbed into the model. In our label ranking experiment, we find that indeed tuning $\varepsilon$ is not necessary to achieve excellent accuracy. On the other hand, for top-$k$ classification, we find that applying a logistic map to squash $\bm{\theta}$ to $^{n}$ and tuning $\varepsilon$ is important, confirming the empirical finding of Cuturi et al. (2019).

Relation to linear assignment formulation.

Similarly, we can rewrite (7) as $s(\bm{\theta})=\bm{P}(\bm{\theta})^{\top}\bm{\theta}$. To obtain a differentiable operator, Cuturi et al. (2019) (see also (Adams & Zemel, 2011)) propose to replace the permutation matrix $\bm{P}(\bm{\theta})$ by a doubly stochastic matrix $\bm{P}_{\varepsilon E}(\bm{\theta})\coloneqq\operatorname*{argmin}_{\bm{P}\in\mathcal{B}}~{}\langle\bm{P},D(-\bm{\theta},\bm{\rho})\rangle+\varepsilon E(\bm{P})$, which is computed approximately in $O(Tn^{2})$ using Sinkhorn (1967). In comparison, our approach is based on regularizing ${\bm{y}}=\bm{P}\bm{\rho}$ with $\Psi\in\{Q,E\}$ directly, the key to achieve $O(n\log n)$ time and $O(n)$ space complexity, as we now show.

Fast computation and differentiation

As shown in the previous section, computing our soft sorting and ranking operators boils down to projecting onto a permutahedron. Our key contribution in this section is the derivation of an $O(n\log n)$ forward pass and an $O(n)$ backward pass (multiplication with the Jacobian) for these projections. Beyond soft sorting and ranking, this is an important sensitivity analysis question in its own right.

We now show how to reduce the projections to isotonic optimization, i.e., with simple chain constraints, which is the key to fast computation and differentiation. We will w.l.o.g. assume that ${\bm{w}}$ is sorted in descending order (if not the case, we sort it first). {proposition}Reduction to isotonic optimization

The function $\bm{v}_{Q}$ is classically known as isotonic regression. The fact that it can be used to solve the Euclidean projection onto $\mathcal{P}({\bm{w}})$ has been noted several times (Negrinho & Martins, 2014; Zeng & Figueiredo, 2015). The reduction of Bregman projections, which we use here, to isotonic optimization was shown by Lim & Wright (2016). Unlike that study, we use the KL projection of $e^{\bm{z}}$ onto $\mathcal{P}(e^{\bm{w}})$, and not of $\bm{z}$ onto $\mathcal{P}({\bm{w}})$, which simplifies many expressions. We include in §B.4 a simple unified proof of Section 5 based on Fenchel duality and tools from submodular optimization. We also discuss an interpretation of adding regularization to the primal linear program as relaxing the equality constraints of the dual linear program in §B.5.

Computation.

Hence, PAV returns an exact solution of both $\bm{v}_{Q}(\bm{s},{\bm{w}})$ and $\bm{v}_{E}(\bm{s},{\bm{w}})$ in $O(n)$ time (Best et al., 2000). This means that we do not need to choose a number of iterations or a level of precision, unlike with Sinkhorn. Since computing $P_{Q}(\bm{z},{\bm{w}})$ and $P_{E}(\bm{z},{\bm{w}})$ requires obtaining $\bm{s}=\bm{z}_{\sigma(\bm{\theta})}$ beforehand, the total computational complexity is $O(n\log n)$.

Differentiating isotonic optimization.

The block-wise structure of the solution also makes its derivatives easy to analyze, despite the fact that we are differentiating the solution of an optimization problem. Since the coordinates of the solution in block $\mathcal{B}_{j}$ are all equal to $\gamma(\mathcal{B}_{j})$, which in turn depends only on a subset of the parameters, the Jacobian has a simple block-wise form, which we now formalize. {lemma}Jacobian of isotonic optimization

Let $\mathcal{B}_{1},\dots,\mathcal{B}_{m}$ be the ordered partition of $[n]$ induced by $\bm{v}_{\Psi}(\bm{s},{\bm{w}})$ from Proposition 5. Then,

There are interesting differences between the two forms of regularization. For quadratic regularization, the Jacobian only depends on the partition $\mathcal{B}_{1},\dots,\mathcal{B}_{m}$ (not on $\bm{s}$) and the blocks have constant value. For entropic regularization, the Jacobian does depend on $\bm{s}$ and the blocks are constant column by column. Both formulations are averaging the incoming gradients, one uniformly and the other weighted.

Differentiating the projections.

We now combine Section 5 with Section 5 to characterize the Jacobians of the projections onto the permutahedron and show how to multiply arbitrary vectors with them in linear time. {proposition}Jacobian of the projections

Let $P_{\Psi}(\bm{z},{\bm{w}})$ be defined in Proposition 5. Then,

where $\bm{J}_{\pi}$ is the matrix obtained by permuting the rows and columns of $\bm{J}$ according to $\pi$, and where

Again, the Jacobian w.r.t. ${\bm{w}}$ is entirely symmetric. Unlike the Jacobian of isotonic optimization, the Jacobian of the projection is not block diagonal, as we need to permute its rows and columns. We can nonetheless multiply with it in linear time by using the simple identity $(\bm{J}_{\pi})\bm{z}=(\bm{J}\bm{z}_{\pi^{-1}})_{\pi}$, which allows us to reuse the $O(n)$ multiplication with the Jacobian of isotonic optimization.

With the Jacobian of $P_{\Psi}(\bm{z},{\bm{w}})$ w.r.t. $\bm{z}$ and ${\bm{w}}$ at hand, differentiating $s_{\varepsilon\Psi}$ and $r_{\varepsilon\Psi}$ boils down to a mere application of the chain rule to (12) and (13). To summarize, we can multiply with the Jacobians of our soft operators in $O(n)$ time and space.

Experiments

We present in this section our empirical findings. NumPy, JAX, PyTorch and Tensorflow versions of our sorting and ranking operators are available at https://github.com/google-research/fast-soft-sort/.

OT (Cuturi et al., 2019): optimal transport formulation.

All-pairs (Qin et al., 2010): noting that $[r(\bm{\theta})]_{i}$ is equivalent to $\sum_{j\neq i}\mathbf{1}[\theta_{i}<\theta_{j}]+1$, one can obtain soft ranks in $O(n^{2})$ by replacing the indicator function with a sigmoid.

Proposed: our $O(n\log n)$ soft ranks $r_{Q}$ and $r_{E}$. Although not used in this experiment, for top-$k$ ranking, the complexity can be reduced to $O(n\log k)$ by computing $P_{\Psi}$ using the algorithm of Lim & Wright (2016).

We use the CIFAR-10 and CIFAR-100 datasets, with $n=10$ and $n=100$ classes, respectively. Following Cuturi et al. (2019), we use a vanilla CNN (4 Conv2D with 2 max-pooling layers, ReLU activation, 2 fully connected layers with batch norm on each), the ADAM optimizer (Kingma & Ba, 2014) with a constant step size of $10^{-4}$, and set $k=1$. Similarly to Cuturi et al. (2019), we found that squashing the scores $\bm{\theta}$ to $^{n}$ with a logistic map was beneficial.

Results.

Our empirical results, averaged over $12$ runs, are shown in Figure 4 (left, center). On both CIFAR-10 and CIFAR-100, our soft rank formulations achieve comparable accuracy to the OT formulation, though significantly faster, as we elaborate below. Similarly to Cuturi et al. (2019), we found that the soft top-$k$ loss slightly outperforms the classical cross-entropy (logistic) loss for these two datasets. However, we did not find that the All-pairs formulation could outperform the cross-entropy loss.

The training times for 600 epochs on CIFAR-100 were 29 hours (OT), 21 hours ($r_{Q}$), 23 hours ($r_{E}$) and 16 hours (All-pairs). Training times on CIFAR-10 were similar. While our soft operators are several hours faster than OT, they are slower than All-pairs, despite its $O(n^{2})$ complexity. This is due the fact that, with $n=100$, All-pairs is very efficient on GPUs, while our PAV implementation runs on CPU.

2 Runtime comparison: effect of input dimension

To measure the impact of the dimensionality $n$ on the runtime of each method, we designed the following experiment.

Results.

Run times for one batch computation with backpropagation disabled are shown in Figure 4 (Right). While their runtime is reasonable in small dimension, OT and All-pairs scale quadratically with respect to the dimensionality $n$ (note the log scale on the $y$-axis). Although slower than a softmax, our formulations scale well, with the dimensionality $n$ having negligible impact on the runtime. OT and All-pairs go out-of-memory starting from $n=2000$ and $n=3000$, respectively. With backpropagation enabled, they go out-of-memory at $n=1000$ and $n=2500$, due to the need for recording the computational graph. This shows that the lack of memory available on GPUs is problematic for these methods. In contrast, our approaches only require $O(n)$ memory and comes with the theoretical Jacobian (they do not rely on differentiating through iterates). They therefore suffer from no such issues.

3 Label ranking via soft Spearman’s rank correlation coefficient

We now consider the label ranking setting where supervision is given as full rankings (e.g., $2\succ 1\succ 3\succ 4$) rather than as label relevance scores. The goal is therefore to learn to predict permutations, i.e., a function $f_{\bm{w}}\colon\mathcal{X}\to\Sigma$. A classical metric between ranks is Spearman’s rank correlation coefficient, defined as the Pearson correlation coefficient between the ranks. Maximizing this coefficient is equivalent to minimizing the squared loss between ranks. A naive idea would be therefore to use as loss $\frac{1}{2}\|\bm{r}-r(\bm{\theta})\|^{2}$, where $\bm{\theta}=g_{\bm{w}}(\bm{x})$. This is unfortunately a discontinuous function of $\bm{\theta}$. We therefore propose to rather use $\frac{1}{2}\|\bm{r}-r_{\Psi}(\bm{\theta})\|^{2}$, hence the name differentiable Spearman’s rank correlation coefficient. At test time, we replace $r_{\Psi}$ with $r$, which is justified by the order-preservation property (Proposition 4).

We consider the $21$ datasets from (Hüllermeier et al., 2008; Cheng et al., 2009), which has both semi-synthetic data obtained from classification problems, and real biological measurements. Following (Korba et al., 2018), we average over two 10-fold validation runs, in each of which we train on 90% and evaluate on 10% of the data. Within each repetition, we run an internal 5-fold cross-validation to grid-search for the best parameters. We consider linear models of the form $g_{W,\bm{b}}(\bm{x})=W\bm{x}+\bm{b}$, and for ablation study we drop the soft ranking layer $r_{\Psi}$.

Results.

Due to the large number of datasets, we choose to present a summary of the results in Figure 5. We postpone detailed results to the appendix (Table 1). Out of 21 datasets, introducing a soft rank layer with $\Psi=Q$ works better on 15 datasets, similarly on 4 and worse on 2 datasets. We can thus conclude that even for such simple model, introducing our layer is beneficial, and even achieving state of the art results on some of the datasets (full details in the appendix).

4 Robust regression via soft least trimmed squares

The classical ridge regression can be cast as

To empirically validate our proposal, we compare cross-validated results for increasing percentage of outliers of the following methods:

Least trimmed squares, with truncation parameter $k$,

Soft least trimmed squares (26), with truncation parameter $k$ and regularization parameter $\varepsilon$,

Ridge regression (25), with regularization parameter $\varepsilon$,

Huber loss (Huber, 1964) with regularization parameter $\varepsilon$ and threshold parameter $\tau$, as implemented in scikit-learn (Pedregosa et al., 2011).

We consider datasets from the LIBSVM archive (Fan & Lin, 2011). We hold out 20% of the data as test set and use the rest as training set. We artifically create outliers, by adding noise to a certain percentage of the training labels, using $y_{i}\leftarrow y_{i}+e$, where $e\sim\mathcal{N}(0,5\times\text{std}({\bm{y}}))$. We do not add noise to the test set. For all methods, we use L-BFGS (Liu & Nocedal, 1989), with a maximum of $300$ iterations. For hyper-parameter optimization, we use $5$-fold cross-validation. We choose $k$ from $\{\lceil 0.1n\rceil,\lceil 0.2n\rceil,\dots,\lceil 0.5n\rceil\}$, $\varepsilon$ from $10$ log-spaced values between $10^{-3}$ and $10^{4}$, and $\tau$ from $5$ linearly spaced values between $1.3$ and $2$. We repeat this procedure $10$ times with a different train-test split, and report the averaged $R_{2}$ scores (a.k.a. coefficient of determination).

Results.

The averaged $R_{2}$ scores (higher is better) are shown in Figure 6.3. On all datasets, the accuracy of ridge regression deteriorated significantly with increasing number of outliers. Least trimmed squares (hard or soft) performed slightly worse than the Huber loss on housing, comparably on bodyfat and much better on cadata. We found that hard least trimmed squares (i.e., $\varepsilon=0$) worked well on all datasets, showing that regularization is less important for sorting operators (which are piecewise linear) than for ranking operators (which are piecewise constant). Nevertheless, regularization appeared useful in some cases. For instance, on cadata, the cross-validation procedure picked $\varepsilon>1000$ when the percentage of outliers is less than $20\%$, and $\varepsilon<10^{-3}$ when the percentage of outliers is larger than $20\%$. This is confirmed visually on Figure 6.3, where the soft sort with $\Psi=Q$ works slightly better than the hard sort with few outliers, then performs comparably with more outliers. The interpolation effect enabled by $\varepsilon$ therefore allows some adaptivity to the (unknown) percentage of outliers.

Conclusion

Building upon projections onto permutahedra, we constructed differentiable sorting and ranking operators. We derived exact $O(n\log n)$ computation and $O(n)$ differentiation of these operators, a key technical contribution of this paper. We demonstrated that our operators can be used as a drop-in replacement for existing $O(n^{2})$ ones, with an order-of-magnitude speed-up. We also showcased two applications enabled by our soft operators: label ranking with differentiable Spearman’s rank correlation coefficient and robust regression via soft least trimmed squares.

Acknowledgements

We are grateful to Marco Cuturi and Jean-Philippe Vert for useful discussions, and to Carlos Riquelme for comments on a draft of this paper.

References

Appendix A Additional empirical results

We include in this section the detailed label ranking results on the same $21$ datasets as considered by Hüllermeier et al. (2008) as well as Cheng et al. (2009).

Spearman’s rank correlation coefficient for each method, averaged over $5$ runs, is shown in the table below.

Appendix B Proofs

and in particular for ${\bm{w}}=\bm{\rho}$. The second claim follows from

B.2 Proof of Proposition 6 (Linear programming formulations)

where in the second equality we used $\mathcal{P}(\bm{\theta})=\operatorname*{conv}(\Sigma(\bm{\theta}))$ and the fundamental theorem of linear programming (Dantzig et al., 1955, Theorem 6). For the second claim, we have similarly

Setting ${\bm{w}}=\bm{\rho}$ and using $\bm{\rho}_{r(\bm{\theta})}=\bm{\rho}_{{\sigma^{-1}}(\bm{\theta})}={\sigma^{-1}}(-\bm{\theta})=r(-\bm{\theta})$ proves the claim.

B.3 Proof of Proposition 4 (Properties of soft sorting and ranking operators)

Let $\mathcal{C}$ be a closed convex set and let $\bm{\mu}^{\star}(\bm{z})\coloneqq\operatorname*{argmax}_{\bm{\mu}\in\mathcal{C}}\langle\bm{\mu},\bm{z}\rangle-\Psi(\bm{z})$. If $\Psi$ is strongly convex over $\mathcal{C}$, then $\bm{\mu}^{\star}(\bm{z})$ is Lipschitz continuous. By Rademacher’s theorem, $\bm{\mu}^{\star}(\bm{z})$ is differentiable almost everywhere. Furthermore, since $P_{\Psi}(\bm{z},{\bm{w}})=\nabla\Psi(\bm{\mu}^{\star}(\bm{z}))$ with $\mathcal{C}=\mathcal{P}(\nabla\Psi^{-1}({\bm{w}}))$, $P_{\Psi}(\bm{z},{\bm{w}})$ is differentiable a.e. as long as $\Psi$ is twice differentiable, which is the case when $\Psi\in\{Q,E\}$.

Order preservation.

Proposition 1 of Blondel et al. (2019) shows that $\bm{\mu}^{\star}(\bm{z})$ and $\bm{z}$ are sorted the same way. Furthermore, since $P_{\Psi}(\bm{z},{\bm{w}})=\nabla\Psi(\bm{\mu}^{\star}(\bm{z}))$ with $\mathcal{C}=\mathcal{P}(\nabla\Psi^{-1}({\bm{w}}))$ and since $\nabla\Psi$ is monotone, $P_{\Psi}(\bm{z},{\bm{w}})$ is sorted the same way as $\bm{z}$, as well. Let $\bm{s}=s_{\varepsilon\Psi}(\bm{\theta})$ and $\bm{r}=r_{\varepsilon\Psi}(\bm{\theta})$. From the respective definitions, this means that $\bm{s}$ is sorted the same way as $\bm{\rho}$ (i.e., it is sorted in descending order) and $\bm{r}$ is sorted the same way as $-\bm{\theta}$, which concludes the proof.

Asymptotic behavior.

We will now characterize the behavior for sufficiently small and large regularization strength $\varepsilon$. Note that rather than multiplying the regularizer $\Psi$ by $\varepsilon>0$, we instead divide $\bm{s}$ by $\varepsilon$, which is equivalent. {lemma}Analytical solutions of isotonic optimization in the limit regimes

If $\varepsilon\leq\varepsilon_{\text{min}}(\bm{s},{\bm{w}})\coloneqq\min_{i\in[n-1]}\frac{s_{i}-s_{i+1}}{w_{i}-w_{i+1}}$, then

If $\varepsilon>\varepsilon_{\text{max}}(\bm{s},{\bm{w}})\coloneqq\max_{i<j}\frac{s_{i}-s_{j}}{w_{i}-w_{j}}$, then

where $\text{LSE}(\bm{x})\coloneqq\log\sum_{i}e^{x_{i}}$.

We start with the $\varepsilon\leq\varepsilon_{\text{min}}(\bm{s},{\bm{w}})$ case. Recall that $\bm{s}$ is sorted in descending order. Therefore, since we chose $\varepsilon$ sufficiently small, the vector $\bm{v}=\bm{s}/\varepsilon-{\bm{w}}$ is sorted in descending order as well. This means that $\bm{v}$ is feasible, i.e., it belongs to the constraint sets in Section 5. Further, note that $v_{i}=\gamma_{Q}(\{i\};\bm{s}/\varepsilon,{\bm{w}})=\gamma_{E}(\{i\};\bm{s}/\varepsilon,{\bm{w}})=s_{i}/\varepsilon-w_{i}$ so that $\bm{v}$ is the optimal solution if we drop the constraints, which completes the argument.

Next, we tackle the $\varepsilon>\varepsilon_{\text{max}}(\bm{s},{\bm{w}})$ case. Note that the claimed solutions are exactly $\gamma_{Q}([n];\bm{s},{\bm{w}})$ and $\gamma_{E}([n];\bm{s},{\bm{w}})$, so the claim will immediately follow if we show that $[n]$ is an optimal partition. The PAV algorithm (cf. §B.6) merges at each iteration any two neighboring blocks $B_{1},B_{2}$ that violate $\gamma_{\Psi}(B_{1};\bm{s}/\varepsilon,{\bm{w}})\geq\gamma_{\Psi}(B_{2};\bm{s}/\varepsilon,{\bm{w}})$, starting from the partitions consisting of singleton sets. Let $k\in\{1,\dots,n-1\}$ be the iteration number. We claim that the two blocks, $B_{1}=\{1,2,\ldots,k\}$ and $B_{2}=\{k+1\}$, will always be violating the constraint, so that they can be merged. Note that in the quadratic case, they can be merged only if

which is indeed satisfied when $\varepsilon>\varepsilon_{\text{max}}(\bm{s},{\bm{w}})$. In the KL case, they can be merged only if

This will be true if the $i^{\text{th}}$ term on the left-hand side is smaller than the $i^{\text{th}}$ term on the right-hand side, i.e., when $(s_{i}-s_{k+1})/\varepsilon<w_{i}-w_{k+1}$, which again is implied by the assumption. ∎

We can now directly characterize the behavior of the projection operator $P_{\Psi}$ in the two regimes $\varepsilon\leq\varepsilon_{\text{min}}(s(\bm{z}),{\bm{w}})$ and $\varepsilon>\varepsilon_{\text{max}}(s(\bm{z}),{\bm{w}})$. This in turn implies the results for both the soft ranking and sorting operations using (12) and (13). {proposition}Analytical solutions of the projections in the limit regimes

If $\varepsilon\leq\varepsilon_{\text{min}}(s(\bm{z}),{\bm{w}})$, then

If $\varepsilon>\varepsilon_{\text{max}}(s(\bm{z}),{\bm{w}})$, then

Therefore, in these two regimes, we do not even need PAV to compute the optimal projection.

B.4 Proof of Proposition 5 (Reduction to isotonic optimization)

Before proving Proposition 5, we need the following three lemmas.

Note that $s_{2}-v_{2}\geq s_{2}-v_{1}\geq s_{1}-v_{1}$ and $s_{2}-v_{2}\geq s_{1}-v_{2}\geq s_{1}-v_{1}$. This means that we can express $s_{2}-v_{1}$ and $s_{1}-v_{2}$ as a convex combination of the endpoints of the line segment $[s_{1}-v_{1},s_{2}-v_{2}]$, namely

Solving for $\alpha$ and $\beta$ gives $\alpha=1-\beta$. From the convexity of $f$, we therefore have

Dual formulation of a regularized linear program

Let $\Omega\coloneqq\Psi+\Phi$, where $\Psi$ is strongly convex and $\Phi$ is convex. We have

which is the infimal convolution of $\Phi^{*}$ with $\Psi^{*}$. Moreover, $\nabla\Omega^{*}(\bm{z})=\nabla\Psi^{*}(\bm{z}-\bm{u}^{\star})$. The results follows from choosing $\Phi(\bm{\mu})=I_{\mathcal{C}}(\bm{\mu})$ and noting that $I_{\mathcal{C}}^{*}=s_{\mathcal{C}}$. ∎

For instance, with $\Psi=Q$, we have $\Psi^{*}=Q$, and with $\Psi=E$, we have $\Psi^{*}=\exp$.

The next lemma shows how to go further by choosing $\mathcal{C}$ as the base polytope $\mathcal{B}(F)$ associated with a cardinality-based submodular function $F$, of which the permutahedron is a special case. The polytope is defined as (see, e.g., Bach (2013))

Reducing dual formulation to isotonic regression

The support function $s_{\mathcal{B}(F)}(\bm{u})$ is known as the Lovász extension of $F$. For conciseness, we use the standard notation $f(\bm{u})\coloneqq s_{\mathcal{B}(F)}(\bm{u})$. Applying Lemma B.4, we obtain

Moreover, $\langle\bm{f},\bm{u}\rangle=f(\bm{u})$.

Let us fix $\sigma$ to the permutation that sorts $\bm{u}^{\star}$. Following the same idea as from (Djolonga & Krause, 2017), since the Lovász extension is linear on the set of all vectors that are sorted by $\sigma$, we can write

This is an instance of isotonic optimization, as we can rewrite the problem as

with $\bm{u}_{\sigma}^{\star}=\bm{v}^{\star}\Leftrightarrow\bm{u}^{\star}=\bm{v}_{\sigma^{-1}}^{\star}$.

Let $\bm{s}\coloneqq\bm{z}_{\sigma}$. It remains to show that $s_{1}\geq\dots\geq s_{n}$, i.e., that $\bm{s}$ and the optimal dual variables $\bm{v}^{\star}$ are both in descending order. Suppose $s_{j}>s_{i}$ for some $i<j$. Let $\bm{s}^{\prime}$ be a copy of $\bm{s}$ with $s_{i}$ and $s_{j}$ swapped. Since $\psi^{*}$ is convex, by Lemma B.4,

which contradicts the assumption that $\bm{v}^{\star}$ and the corresponding $\sigma$ are optimal. A similar result is proven by Suehiro et al. (2012, Lemma 1) but for the optimal primal variable $\bm{\mu}^{\star}$. ∎

We now prove Proposition 5. The permutahedron $\mathcal{P}({\bm{w}})$ is a special case of $\mathcal{B}(F)$ with $F(\mathcal{S})=\sum_{i=1}^{|\mathcal{S}|}w_{i}$ and $w_{1}\geq w_{2}\geq\dots\geq w_{n}$. In that case, $\bm{f}_{\sigma}=(f_{\sigma_{1}},\dots,f_{\sigma_{n}})=(w_{1},\dots,w_{n})={\bm{w}}$.

For $\mathcal{P}(\nabla\Psi^{*}({\bm{w}}))$, we thus have $\bm{f}_{\sigma}=\nabla\Psi^{*}({\bm{w}})$. Finally, note that if $\Psi$ is Legendre-type, which is the case of both $Q$ and $E$, then $\nabla\Psi^{*}=(\nabla\Psi)^{-1}$. Therefore, $\nabla\Psi(\bm{\mu}^{\star})=\bm{z}-\bm{u}^{\star}$, which concludes the proof.

B.5 Relaxed dual linear program interpretation

We show in this section that the dual problem in Lemma B.4 can be interpreted as the original dual linear program (LP) with relaxed equality constraints. Consider the primal LP

As shown by Bach (2013, Proposition 3.2), the dual LP is

Moreover, let $\sigma$ be a permutation sorting $\bm{z}$ in descending order. Then, an optimal $\bm{\lambda}$ is given by (Bach, 2013, Proposition 3.2)

Now let us restrict to the support of $\bm{\lambda}$ and do the change of variable

The non-negativity constraints in $\mathcal{C}$ become $v_{1}\geq v_{2}\geq\dots\geq v_{n}$ and the equality constraints in $\mathcal{C}$ become $\bm{z}_{\sigma}=\bm{v}$. Adding quadratic regularization $\frac{1}{2}\|{\bm{y}}\|^{2}$ in the primal problem (55) is equivalent to relaxing the dual equality constraints in (56) by smooth constraints $\frac{1}{2}\|\bm{z}_{\sigma}-\bm{v}\|^{2}$ (this can be seen by adding quadratic regularization to the primal variables of Bach (2013, Eq. (3.6))). For the dual objective (56), we have

where in the second line we used (Bach, 2013, Eq. (3.2)). Altogether, we obtain $\min_{v_{1}\geq\dots\geq v_{n}}\frac{1}{2}\|\bm{z}_{\sigma}-\bm{v}\|^{2}+\langle\bm{f}_{\sigma},\bm{v}\rangle$, which is exactly the expression we derived in Lemma B.4. The entropic case is similar.

B.6 Pool adjacent violators (PAV) algorithm

Let $g_{1},\dots,g_{n}$ be convex functions. As shown in (Best et al., 2000; Lim & Wright, 2016),

can be solved using a generalization of the PAV algorithm (note that unlike these works, we use decreasing constraints for convenience). All we need is a routine for solving, given some set $\mathcal{B}$ of indices, the “pooling” sub-problem

Thus, we can use PAV to solve (53), as long as $\Psi^{*}$ is separable. We now give the closed-form solution for two special cases. To simplify, we denote $\bm{s}\coloneqq\bm{z}_{\sigma}$ and ${\bm{w}}\coloneqq\bm{f}_{\sigma}$.

We have $g_{i}(v_{i})=\frac{1}{2}(s_{i}-v_{i})^{2}+v_{i}w_{i}$. We therefore minimize

Entropic regularization.

We have $g_{i}(v_{i})=e^{s_{i}-v_{i}}+v_{i}e^{w_{i}}$. We therefore minimize

where $\text{LSE}(\bm{x})\coloneqq\log\sum_{i}e^{x_{i}}$.

Although not explored in this work, other regularizations are potentially possible, see, e.g., (Blondel et al., 2019).

B.7 Proof of Proposition 5 (Jacobian of isotonic optimization)

Let $\mathcal{B}_{1},\dots,\mathcal{B}_{m}$ be the partition of $[n]$ induced by $\bm{v}\coloneqq\bm{v}_{\Psi}(\bm{s},{\bm{w}})$. From the PAV algorithm, for all $i\in[n]$, there is a unique block $\mathcal{B}_{l}\in\{\mathcal{B}_{1},\dots,\mathcal{B}_{m}\}$ such that $i\in\mathcal{B}_{l}$ and $v_{i}=\gamma_{\Psi}(\mathcal{B}_{l};\bm{s},{\bm{w}})$. Therefore, for all $i\in[n]$, we obtain

Therefore, the Jacobian matrix is block diagonal, i.e.,

The multiplication with the Jacobian uses the fact that each block is constant column-wise.

The expression above is for points $\bm{s}$ where $\bm{v}$ is differentiable. For points where $\bm{v}$ is not differentiable, we can take an arbitrary matrix in the set of Clarke’s generalized Jacobians, the convex hull of Jacobians of the form $\lim_{\bm{s}_{t}\to\bm{s}}\partial\bm{v}/\partial\bm{s}_{t}$. The points of non-differentiability occur when a block of the optimal solution can be split up into two blocks with equal values. In that case, the two directional derivatives do not agree, but are derived for quadratic regularization by Djolonga & Krause (2017).