Spectral Ranking Regression

Rank Aggregation. The problem of rank aggregation [Dwork et al. 2001], in which a total ordering of samples is regressed from ranking observations, is classic; literature on the subject is vast—see, e.g., the surveys by Fligner and Verducci [1993], Cattelan [2012], and Marden [2014]. Probabilistic inference in this setting typically assumes (a) that a “true” total ordering of samples exists, and (b) that ranking observations exhibit a stochastic transitivity property [Agarwal 2016]: A sample is more likely to be ranked higher than another when this event is consistent with the underlying total ordering. The noisy permutation model is a non-parametric model for this setting: Pairwise comparisons consistent with the underlying total ordering are observed under i.i.d. Bernoulli noise. MLE is NP-hard in this setting. A polytime algorithm by Braverman and Mossel [2008] recovers the underlying ordering in \(\Theta (n\log n)\) comparisons, w.h.p.; this is tightened by several recent works [Wauthier et al. 2013; Mao et al. 2018a;, 2018b]. The Mallows model [Mallows 1957] assumes that the probability of a ranking observation is a decreasing function of its distance from the underlying total ordering, under appropriate notions of distance (e.g., Kendall-Tau (KT)); MLE can be approached, e.g., via EM [Lu and Boutilier 2011]. Shah et al. [2016b] learn the full matrix of pairwise comparison probabilities via a minimax optimal estimator. Rajkumar and Agarwal [2016] learn the matrix via matrix completion, requiring \(\Theta (nr \log n)\) comparisons, where \(r \ll n\) is the rank. Ammar and Shah [2011] assume that comparisons are sampled from an unknown distribution over total orderings and propose an entropy maximization algorithm requiring \(\Theta (n^2)\) comparisons.

Parametric Rank Aggregation. We focus on parametric models, as they are more natural in the context of regressing rankings from sample features. In both PL [Plackett 1975] and Thurstone [Thurstone 1927], each sample is parametrized by a score. In the Thurstone model, observations result from comparing scores after the addition of Gaussian noise. Vojnovic and Yun [2016] and Shah et al. [2016a] estimate Thurstone scores via MLE and provide sample complexity bounds that are inversely proportional to the smallest non-zero eigenvalue of the Laplacian of a graph modeling comparisons. In PL, the probability that a sample is chosen over a set of alternatives is proportional to its score. Hunter [2004] proposes a Minorization-Maximization (MM) approach to estimate PL scores via MLE, earlier used by Dykstra [1960] on pairwise comparisons (i.e., on the Bradley–Terry (BT) setting). Hajek et al. [2014] provide an upper bound on the error in estimating the PL scores via MLE and show that the latter is minimax-optimal. Negahban et al. [2018] propose a latent factor model, estimating parameters via a convex relaxation of the corresponding rank penalty and providing sample complexity guarantees. Vojnovic et al. [2020] show that gradient-descent and MM algorithms that solve the MLE problem converge with linear rate.

PL Inference via Spectral Methods. Our focus on PL is due to the recent emergence of spectral algorithms for inference in this setting. Negahban et al. [2012] propose the Rank Centrality (RC) algorithm for the BT setting and derive a minimax error bound. Chen and Suh [2015] propose a spectral MLE algorithm extending RC with an additional stage that cyclically performs MLE for each score. Soufiani et al. [2013] and Jang et al. [2017] extend RC to rankings of two or more samples by breaking rankings into independent comparisons. Improved bounds, applying also to broader noise settings, are provided by Rajkumar and Agarwal [2014]. Khetan and Oh [2016] generalize the work by Soufiani et al. [2013] by breaking rankings into independent shorter rankings, and building a hierarchy of tractable and consistent estimators. Blanchet et al. [2016] model sequential choices by state transitions in an MC, where transitions are functions of choice probabilities.

Bridging the above approaches with MLE, Maystre and Grossglauser [2015] show that the MLE of PL scores can be expressed as the stationary distribution of an MC. Their proposed ILSR algorithm estimates the PL scores faster than traditional optimization methods, such as, e.g., Hunter’s [Hunter 2004] and Newton’s method [Nocedal and Wright 2006], and more accurately than prior spectral rank aggregation methods. More recently, Ragain and Ugander [2016] showed that a spectral approach applies even after relaxing the assumption that the relative order of any two samples is independent of the alternatives. Agarwal et al. [2018] proposed another spectral method called accelerated spectral ranking by departing from the exact equivalence between MLE and MC approximation and demonstrate faster convergence than ILSR.

Ranking Regression via PL. We depart from all aforementioned methods by regressing ranked choices from sample features. Ranking regression via the PL [Plackett 1975] model has a long history, mostly focusing on pairwise comparisons, i.e., restricted to \(K=2\); this case is also known as the BT model [Bradley and Terry 1952]. The resulting algorithms have been employed in a broad array of domains including medicine [Tian et al. 2019], information retrieval on search engines [Joachims 2002; Burges et al. 2005], image-based aesthetic recommendations [Chang et al. 2016; Sun et al. 2017], click prediction for online advertisements [Sculley 2010], and image retrieval [Chen et al. 2015], to name a few.

Among shallow regressors, RankSVM [Joachims 2002] learns a target ranking from features via a linear Support Vector Machine (SVM), with constraints imposed by all possible comparisons. Pahikkala et al. [2009] propose a regularized least-squares based algorithm for learning to rank from comparisons. Several works [Joachims 2002; Pahikkala et al. 2009; Guo et al. 2018; Tian et al. 2019] regress comparisons via MLE over BT [Bradley and Terry 1952] models.

Deeper models for regressing comparisons have also been extensively explored in the BT setting [Burges et al. 2005; Chang et al. 2016; Dubey et al. 2016; Doughty et al. 2018; Yıldız et al. 2019]. Burges et al. [2005] estimate comparisons via a fully connected neural network called RankNet, using the BT model to construct a cross-entropy loss. Several works [Chang et al. 2016; Dubey et al. 2016; Doughty et al. 2018; Yıldız et al. 2019] learn from comparisons via siamese networks [Bromley et al. 1994]. Chang et al. [2016] and Yıldız et al. [2019] learn from comparisons using MLE on the BT model. Dubey et al. [2016] predict image comparisons, via a loss function combining cross-entropy and hinge loss. Doughty et al. [2018] learn similarities and comparisons of videos via hinge loss. All of the above methods generalize to rankings under the PL model, but suffer from the complexity issues outlined in the introduction.

Cao et al. [2007], Xia et al. [2008], and Ma et al. [2020] focus on learning from star/relevance ratings rather than rankings; they train a neural network called ListNet that treats ratings as PL scores. Particularly, given a dataset of \(n\) ratings, ListNet constructs all \(\frac{n!}{(n-K)!}\) possible \(K\)-sized rankings; it is then trained by minimizing the cross entropy between the PL probabilities that a given \(K\)-ranking comprises the top samples among all \(n\), with scores being (a) the predicted scores and (b) the ground-truth ratings, respectively. Hence, although ListNet solves a very different problem than ranking regression, it is related to the siamese methods mentioned above through its (exponential in \(K\)) objective.

ADMM Applications in DNN Training. When training DNNs, ADMM is traditionally used to enforce sparsity. One application is compressed sensing, in which a signal is recovered from undersampled measurements via a DNN transformation [Sun et al. 2016; Li et al. 2018; Liu et al. 2018; Ma et al. 2019; Yang et al. 2020]. The training objective combines a standard regression loss with a norm penalty on measurements, including, e.g., \(\ell _1\), \(\ell _2\), Frobenius, or nuclear norm. Another application is model pruning [Ye et al. 2019;, 2019; Zhao and Liao 2019]: ADMM is used to enforce weight sparsity by replacing norm penalties with (possibly non-convex) low-cardinality set constraints. Zhao et al. [2018] train a DNN classifier to recover samples distorted by additive adversarial noise, using an objective that combines a classification loss with an \(\ell _0\), \(\ell _1\), \(\ell _2\), or \(\ell _\infty\) norm penalty.

Relationship to Previous Works. The present article is a combination and extension of two of our earlier works [Yıldız et al. 2020;, 2021], which dealt with shallow and deep model settings separately. We unify the key theoretical results of these works, and include all proofs omitted from these conference versions. We extend our experimental evaluations, and provide a more detailed technical preliminary.

As introduced in Section 2.2.1, MLE of the PL scores \(\mathbf {\pi }\in \mathbb {R}^n_+\) amounts to minimizing the negative log-likelihood:

To regress scores \(\mathbf {\pi }\) from sample features \(\mathbf {X}= {[\mathbf {x}_1,..,\mathbf {x}_n]}^T \in \mathbb {R}^{n\times d}\), we assume that there exists a function \(\tilde{{\pi }}:\mathbb {R}^d\rightarrow [0,1]\), parametrized by \(\mathbf {w}\in \mathbb {R}^{d^{\prime }}\), such thatThen, MLE amounts to minimizing the following objective w.r.t. \(\mathbf {w}\in \mathbb {R}^{d^{\prime }}\):We tackle the regression problem via the MLE objective (15) w.r.t. three regression functions. Our notation is summarized in Table 1.

Affine Regression. We assume that there exist \(\mathbf {a}\in \mathbb {R}^d\) and \({b}\in \mathbb {R}\), such that \(\mathbf {\pi }= \tilde{{\pi }}(\mathbf {X}; \mathbf {w})\equiv \mathbf {X}\mathbf {a}+ \mathbf {1} {b}\,\) for \(\mathbf {w}=(\mathbf {a},{b})\in \mathbb {R}^{d+1}\). Then, MLE of parameters \(\mathbf {w}\) amounts to solving:where \(\mathcal {L}\) is given by (15). We note that Problem (16) is not convex, as the objective is not convex in \((\mathbf {a},{b})\in \mathbb {R}^{d+1}\).

Logistic Regression. In the logistic case, we assume that there exists \(\mathbf {w}\in \mathbb {R}^d\) s.t. \(\mathbf {\pi }=\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\equiv [e^{\mathbf {w}^\top \mathbf {x}_i}]_{i\in \left[n\right]}.\) As \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\ge 0\) by definition, MLE corresponds toThe objective of (17) is convex in \(\mathbf {w}\), so an optimal solution can be found via, e.g., Newton’s method [Nocedal and Wright 2006]. The convexity of Problem (17) w.r.t. \(\mathbf {w}\) follows by the facts that (a) the reparametrization \({\pi }_i = e^{\theta _i}\), \(i \in \left[n\right]\) makes \(\mathcal {L}\) convex, and (b) the composition of convex and affine is convex [Boyd and Vandenberghe 2004].

DNN Regression. Finally, we regress scores \(\mathbf {\pi }\) via a generic DNN function \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) by solving the following MLE problem:where \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\in \mathbb {R}^n_+\) is the vector map whose coefficients \(\tilde{{\pi }}_i\in \mathbb {R}_+\) are given by \(\tilde{{\pi }}_i=\tilde{{\pi }}(\mathbf {x}_{i}; \mathbf {w})\), i.e., \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\equiv [\tilde{{\pi }}(\mathbf {x}_{i}; \mathbf {w})]_{i \in \left[n\right]}\). For instance, \(\tilde{{\pi }}(\cdot , \mathbf {w})\) can be a fully connected feed-forward neural network if \(\mathbf {X}\) contains numerical features, or a convolutional neural network if \(\mathbf {X}\) is an image.

As also discussed in Section 1, virtually all state-of-the-art ranking regression methods minimize Equation (15) via classic gradient-based optimization [Joachims 2002; Burges et al. 2005; Pahikkala et al. 2009; Chang et al. 2016; Dubey et al. 2016; Doughty et al. 2018; Guo et al. 2018; Tian et al. 2019; Yıldız et al. 2019]. Note that the optimization of objective (15) implies the following process: For each observation \((c_\ell ,A_\ell)\) in \(\mathcal {D}\), the regression function \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) needs to be evaluated and updated over all samples in \(A_\ell\). Overall, the ranking regression objective (15) implies that \(K\) samples need to be loaded in RAM to process a single ranking or maximal choice query. To make matters worse, the number \(m\) of possible rankings/maximal queries in \(\mathcal {D}\) grows as \(\left({{n}\atop {K}}\right) = O(n^K)\). This means that the length of a training epoch of stochastic gradient descent (and hence, the number of expensive gradient computations) is \(m=O(n^K)\), i.e., grows exponentially in \(K\). Combined, these two effects can make the cost of a single optimization step over Equation (15) prohibitive (c.f. Figures 1 and 5(c)). This motivates us to explore efficient spectral algorithms to solve regression problems (16)–(18), as we discuss next.

Given ILSR’s significant computational benefits in Section 2, we wish to develop analogues in the regression setting. In contrast to the feature-less setting, it is not a priori evident how to solve Problems (16)–(18) via a spectral approach. Taking the affine case as an example, and momentarily ignoring issues of non-convexity, the stationary points of the optimization problem (16) cannot be expressed via the balance equations of an MC. Our main contribution is to circumvent this problem by using the ADMM [Boyd et al. 2011]. Intuitively, ADMM allows us to decouple the optimization of scores \(\mathbf {\pi }\) from model parameters \(\mathbf {w}\), encapsulating them in a proximal penalty: The latter becomes amenable to a spectral approach after a series of manipulations that we outline below (see Theorem 3.2).

We wish to devise an efficient algorithm to minimize the regression objective (15). To do so, we extend the standard ADMM [Boyd et al. 2011] introduced in Section 2. To this end, we rewrite the minimization of (15) as

To solve (19) via ADMM, consider the following augmented Lagrangian (see Section 2):where \(\mathbf {y}\in \mathbb {R}^n\) is the Lagrangian dual variable corresponding to the equality constraint in (19b), \(\rho \ge 0\) is a penalty parameter, and \(D_p(\cdot ||\cdot)\) is a proximal penalty term. This term satisfies: (i) \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }}) \ge 0\) for all \(\mathbf {\pi },\tilde{\mathbf {\pi }}\in \mathbb {R}^n_+\), and (ii) \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }}) = 0\) if and only if \(\tilde{\mathbf {\pi }}=\mathbf {\pi }\). Classic ADMM involves using the usual \(\ell _2\) proximal penalty, i.e., \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }}) = \frac{1}{2} \parallel \! \mathbf {\pi }- \tilde{\mathbf {\pi }}\!\parallel _2^2\). We depart from this by considering more generic \(D_p\); as we discuss in Section 3.5, using KL divergence instead, i.e., \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }}) = \sum _{i=1}^{n} {\pi }_i \log \frac{{\pi }_i}{\tilde{{\pi }}_i}\) has significant advantages in our setting.

In its general form, ADMM alternates between optimizing \(\mathbf {\pi }\) and \(\mathbf {w}\) until convergence via the following primal-dual algorithm on the augmented Lagrangian (c.f. Section 2):

This has the following immediate computational advantages. First, step (21b) is equivalent toNote that this operation does not depend on dataset \(\mathcal {D}\): The model \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) is regressed directly from the present score estimates \(\mathbf {\pi }^{k+1}\) via penalty \(D_p\), with the additional linear dual term. In particular, as discussed below, \(D_p\) leads to a least-squares regression in the case of the \(\ell _2\)-proximal penalty, and a max-entropy penalty in the case of KL-divergence; both can be solved efficiently. Crucially, an optimization method that solves (22) iterates over the \(O(n)\) samples, instead of the \(O(n^K)\) ranking observations, which would be the case when Equation (15) is minimized directly without ADMM decoupling.

Most importantly, step (21b)—which does depend on \(\mathcal {D}\)—admits a highly efficient spectral implementation for a proximal penalty \(D_p\); ADMM (21) therefore delegates solving the “expensive” portion of the problem to a highly efficient algorithm. To establish this, we begin with the following lemma, which we prove in Appendix C.

Contrasting Equation (26) to Equation (4), it is not immediately evident that the former corresponds to the balance equations of an MC as, in general, \(\mathbf {\sigma }(\mathbf {\pi }) = [\sigma _i(\mathbf {\pi })]_{i \in \left[n\right]} \ne \mathbf {0}\). Nevertheless, for an appropriate definition of transition rates, we prove that this is indeed the case:

The proof is in Appendix D. By Lemma 3.1 and Theorem 3.2, we conclude that a stationary \(\mathbf {\pi }\in \mathbb {R}^n_+\) satisfying (23) is also the stationary distribution of the continuous-time MC with transition rates:where \(\mu _{ji}(\mathbf {\pi })\) are given by Equation (27). Motivated by these observations, and mirroring ILSR (Equation (7)), we compute a solution to (21b) viawhere \(\mathbf {M}(\cdot)\) is given by Equation (28) and \(\mathtt {ssd}(\cdot)\) is an algorithm computing the steady state distribution. We refer to this procedure as ILSRX (“ILSR with features”).

Our spectral algorithm solving Equation (19) is summarized in Algorithm 1 and Figure 2. We initialize all samples with equal scores, i.e., \(\mathbf {\pi }= \tilde{\mathbf {\pi }}= \frac{1}{n}\mathbf {1}\), and set the initial dual variable as \(\mathbf {y}^0 = \mathbf {0}\). We iteratively update \(\mathbf {\pi }\), \(\mathbf {w}\), and \(\mathbf {y}\) via Equation (21) until convergence. At each ADMM iteration \(k+1\), we initialize \(\mathbf {\pi }\) with \(\mathbf {\pi }^{k}\), and update \(\mathbf {\pi }\) via ILSRX given by Equation (29). We elaborate on the resulting implementation specifics for each of the regression problems (16)–(18) below.

In this section, we apply Algorithm 1 on affine regression of PL scores viaWe introduce \(\mathbf {w}=(\mathbf {a}, {b}) \in \mathbb {R}^{d+1}\) and \(\mathbf {\tilde{X}}= [\mathbf {X}| \mathbf {1}] \in \mathbb {R}^{n\times (d+1)}\) for notational simplicity. Recall from Section 3.1 that the resulting MLE problem is given by

In this setting, we employ standard ADMM with quadratic proximal penalty:This leads to additional theoretical guarantees explained in Section 3.7. Employing Algorithm 1 with quadratic \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }})\) and affine \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\), the regression step (21b) of ADMM becomes:Equation (33) follows from completing the square by combining the linear and quadratic terms into a single quadratic objective. This has the following immediate computational advantage: The regression step (33) is a quadratic minimization and admits a closed form solution.

Crucially, the score update step (21b) is amenable to the spectral solution established in Theorem 3.2, whereHaving adjusted the transition matrix \(\mathbf {M}(\mathbf {\pi })\) thusly, \(\mathbf {\pi }\) can be obtained by repeated iterations of ILSRX (29), with \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) and \(D_p\) given by Equations (30) and (32), respectively. We refer to this instance of Algorithm 1 as Spectral Ranking-l2 (SR-l2).

We note that, as Problem (19) is non-convex, selecting a good initialization point is important in practice. We discuss initialization below, leaving the additional theoretical guarantees to Section 3.7.

Initialization. We initialize \(\mathbf {w}\) so that \(\mathbf {\tilde{X}}\mathbf {w}\) is a good approximation of PL scores. We use a technique akin to Saha and Rajkumar [2018], applied to our affine setting. Given a distribution over queries \(A\subseteq \left[n\right]\), let \(P_{ij}=\mathbb {E}_A[c=i|\lbrace i,j\rbrace \subseteq A]\) be the probability that \(i\) is chosen given a query \(A\) that contains both \(i\) and \(j\). By (1), for \(i, j \in \left[n\right]\), \(\frac{{P}_{ij}}{{P}_{ji}} = \frac{{\pi }_i}{{\pi }_j} = \frac{\mathbf {x}_i^\top \mathbf {a}+ {b}}{\mathbf {x}_j^\top \mathbf {a}+ {b}}\), orMotivated by (35), we estimate \(P_{ij}\) empirically from \(\mathcal {D}\), and obtain our initialization \(\mathbf {w}^0=(\mathbf {a}^0,{b}^0)\in \mathbb {R}^{d+1}\) by solving (35) in the least-square sense; that is,Note that this is a convex quadratic program. Given \(\mathbf {w}^0\), we generate the initial PL scores via the affine parametrization \(\mathbf {\pi }^0=\mathbf {\tilde{X}}\mathbf {w}^0\).

We describe here how to apply Algorithm 1 on logistic regression of PL scores viawhere \(\log \mathbf {\pi }=[\log {\pi }_i]_{i\in \left[n\right]}\) is the \(\mathbb {R}^n\) vector generated by applying \(\log\) to \(\mathbf {\pi }\) element-wise. Recall from Section 3.1 that the resulting MLE problem is given bywhich is convex, and can thus be solved by Newton’s method. Nevertheless, we would like to accelerate its computation via a spectral method akin to ILSR.

Mirroring the affine case, we employ Algorithm 1 with quadratic proximal penalty:As a result, the regression step (21b) corresponds towhich is again a quadratic minimization and admits a closed form solution.

More importantly, the score update step (21b) is amenable to the spectral solution established in Theorem 3.2. Mutatis mutandis, following the same manipulations in Lemma 3.1, a stationary point of the objective in each step (21b) can be cast as the stationary distribution of the continuous-time MC with transition rates \(\mu _{ji}(\mathbf {\pi })\), \(i,j\in \left[n\right]\) (Equation (27)), where vector \(\mathbf {\sigma }= [\sigma _i]_{i\in \left[n\right]}\) is now given byHaving adjusted the transition matrix \(\mathbf {M}(\mathbf {\pi })\) thusly, \(\mathbf {\pi }\) can again be obtained by repeated iterations of ILSRX (29), with \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) and \(D_p\) given by Equations (37) and (39), respectively. We refer to this instance of Algorithm 1 as Spectral Ranking-l2-log (SR-l2-log).

Initialization. Similar to the initialization of SR-l2 (c.f. Equation (36)), we initialize \(\mathbf {w}\) so that the initial scores obey the PL model, mirroring the approach by Saha and Rajkumar [2018]. Defining \(P_{ij}\), \(i,j\in \left[n\right]\) the same way, and using the logistic parametrization in Equation (17), we have thatAccordingly, we initialize \(\mathbf {w}\) aswhere \(\hat{P}_{ij}\), \(i,j\in \left[n\right]\), are again empirical estimates obtained from dataset \(\mathcal {D}\). Given \(\mathbf {w}^0\), we generate the initial PL scores via the logistic parametrization \(\mathbf {\pi }^0=[e^{\mathbf {x}_i^T \mathbf {w}^0}]_{i\in \left[n\right]}\).

In this section, we employ Algorithm 1 for DNN regression of PL scores via \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) (c.f. Problem (18)).

To regress scores via a generic DNN \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\), we generalize the traditional quadratic proximal penalty \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }}) = \frac{1}{2} \parallel \! \mathbf {\pi }- \tilde{\mathbf {\pi }}\!\parallel _2^2\) to KL-divergence, which is naturally suited to Problem (19): This is because the scores computed by ILSRX (29) form, by construction, a distribution over \(\left[n\right]\). This generalization to KL penalty is still amenable to a spectral solution via ILSRX (29), as Theorem 3.2 holds for any proximal penalty \(D_p\); in practice, this also leads to a significantly improved performance over an \(\ell _2\)-proximal penalty (by up to \(56\%\) Top-1 accuracy and \(25\%\) KT correlation, as discussed in Section 4.9). As shown in Equation (48) below, the KL-divergence penalty leads to training \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) with a max-entropy loss and an additional linear term at each ADMM iteration. Compared to the \(\ell _2\)-proximal penalty, max-entropy has been observed to converge faster and lead to better fitting [Caruana and Niculescu-Mizil 2004; Golik et al. 2013; Bosman et al. 2020]. Finally, ADMM with Bregman divergence proximal penalties, including KL divergence, has been shown to offer local convergence guarantees for non-convex problems (c.f. Section 2.2.3); this further motivates us to employ KL over other distance measures between probability distributions.

We implement DSR with the two proximal penalty functions \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }})\) mentioned above: KL divergence [Kullback and Leibler 1951] and \(\ell _2\) norm. We describe implementation specifics for each of these cases below.

\(\ell _2\) proximal penalty. For ADMM with quadratic proximal penalty:the regression step (21b) of DSR becomes a least-squares minimization of the form:As a result, we train the DNN regressor \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) via SGD optimization over the loss function (45).

Moreover, the score update step (21b) is amenable to the spectral approach established in Theorem 3.2 and solved via repeated iterations of ILSRX (29), wherefor all \(i \in \left[n\right]\) in the transition rates (27). We refer to this instance of Algorithm 1 as Deep Spectral Ranking-l2 (DSR-l2).

KL proximal penalty. For ADMM with KL proximal penalty:the optimization step (21c) to train \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) corresponds to minimizing a max-entropy loss with an additional linear term:

Moreover, the score update step (21b) is again solved via repeated iterations of ILSRX (29), wherefor all \(i \in \left[n\right]\) in the transition rates (27). We refer to this instance of Algorithm 1 as Deep Spectral Ranking-KL (DSR-KL).

Initialization. For both cases of penalties, we initialize all samples with equal scores, i.e., \(\mathbf {\pi }= \tilde{\mathbf {\pi }}= \frac{1}{n}\mathbf {1}\). At each ADMM iteration \(k\), we initialize \(\mathbf {\pi }\) with \(\mathbf {\pi }^{k-1}\), and update \(\mathbf {\pi }\) via ILSRX. Then, we initialize the DNN parameters \(\mathbf {w}\) with \(\mathbf {w}^{k-1}\), and fine-tune the DNN \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) via SGD over Equations (45) or (48).

Each iteration of our spectral algorithm (c.f. Algorithm 1) involves the three steps in Equation (21). One iteration of ILSRX at step (21b) is \(O(Km+ n^2)\) for constructing the transition rates via Equation (27) and for finding the stationary distribution \(\mathbf {\pi }\) via, e.g., a power method [Lei et al. 2016], respectively. It is important to note that the number of ranking observations \(m\) appears only in the construction of the transition rates through the terms \(\mu _{j, i}\) in Equation (21b), which, in turn, are computed via sums w.r.t. these \(m\) quantities (see Equation (27)). This is important, as the (potentially exponential in \(K\)) ranking observations affect the complexity only through scalar summations, which themselves can be easily parallelized (via, e.g., reduce-by-key operations). This is in sharp contrast to the classic gradient-based optimization (c.f. Section 3.1 and Figure 1), in which \(m\) directly determines the size of each epoch.

We discuss the complexity of the regression step (21c) separately for shallow and DNN model regression cases. For shallow model regression via affine or logistic parametrizations, i.e., SR-l2 and SR-l2-log, the regression step (21c) takes the specific forms (33) and (40), respectively. In both cases, \(\mathbf {w}\) update is \(O\big (n(d+1)\big)\) as a matrix-vector multiplication, since the matrix \({(\mathbf {X}^T\mathbf {X})}^{-1}\mathbf {X}^T\) can be precomputed. For DNN regression via DSR, the regression step (21c) is equivalent to training \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) via Equation (22). Most importantly, the length of each training epoch, is linear in number of samples \(n\): Each optimization step goes over \(O(n)\) samples to update the \({d^{\prime }}\) parameters. Again, this is in contrast to the \(m=O(n^K)\)-length epochs if one was minimizing the likelihood function directly via standard gradient descent methods. Finally, given \(\tilde{\mathbf {\pi }}\) and \(\mathbf {\pi }\), the update of \(\mathbf {y}\) at step (21a) is \(O(n)\).

Overall, updating \(\mathbf {\pi }\) via ILSRX at step (21b) and \(\mathbf {w}\) via gradient-based training at step (21c) are both iterative processes, for which the computational complexity until convergence have not yet been theoretically shown. That said, the literature on the complexity and convergence of iterative optimization over such highly non-convex objectives illustrates the challenges behind such a theoretical analysis [Jain and Kar 2017; Taheri et al. 2021]. Nevertheless, our Theorem 3.4 below, establishing convergence guarantees in the affine regression case, is a step toward this direction. This is in sharp contrast to, e.g., siamese networks, in which \(m\) characterizes the length of an entire training epoch.

As the ranking regression problem (19) is, in general, non-convex, we aim at obtaining convergence guarantees for our spectral algorithm (c.f. Algorithm 1). To do so, we focus on the affine regression case explained in Section 3.3. Given this case, we wish to establish conditions such that (i) the scores computed at each ADMM iteration form a global optimum w.r.t. step (21a), and (ii) the scores and model parameters generated by Algorithm 1 converge to a stationary point of Problem (16).

Recall from Section 3.3 that the affine regression problem is given byIn this setting, we employ ADMM with standard \(\ell _2\) proximal penalty \(D_p(\mathbf {\pi }||\tilde{\mathbf {\pi }}) = \frac{1}{2} \parallel \! \mathbf {\pi }- \tilde{\mathbf {\pi }}\!\parallel _2^2\). As Problem (19) is non-convex, condition (23) is, in general, not sufficient for optimality w.r.t. step (21a). To show this, we require the following technical assumption:

In other words, we assume that the scores computed at each ADMM iteration are strictly positive. Under this assumption, we show that stationarity implies optimality w.r.t. (21a) for large enough \(\rho\):

The proof is in Appendix E. Theorem 3.3 implies that given large enough \(\rho\), the stationary scores satisfying condition (23) also form a global optimum of the ADMM step (21a). Using this result, we further establish the following convergence guarantee for Algorithm 1 in the affine regression case:

The proof is in Appendix F. By Theorem 3.4, we conclude that the scores \(\mathbf {\pi }\) and the affine regression parameters \(\mathbf {w}\) generated by Algorithm 1 converge to a stationary point of Problem (16) for large enough \(\rho\).

We evaluate our algorithms on synthetic and real-life datasets, summarized in Table 2.

Synthetic Datasets. We generate the feature vectors \(\mathbf {x}_i \in \mathbb {R}^{d}, \;\; i \in \left[n\right]\) from \(\mathcal {N} (\mathbf {0},\sigma ^2_{x}\mathbf {I}_{d\times d})\) and a common parameter vector \({\mathbf {w}} \in \mathbb {R}^{d}\) from \(\mathcal {N} (\mathbf {0},\sigma ^2_{{\beta }}\mathbf {I}_{d\times d})\). Then, we generate the PL scores via the logistic parametrization \(\mathbf {\pi }=[e^{\mathbf {x}_i^T \mathbf {w}}]_{i\in \left[n\right]}\). We normalize the resulting scores, so that \(\mathbf {1}^\top \mathbf {\pi }=1\). We set \(\sigma ^2_{x}=\sigma ^2_{{\beta }}=0.8\) in all experiments. Given \(\mathbf {\pi }\), we generate each observation in \(\mathcal {D}\) as follows: We first select \(|A_\ell |=2\) samples out of \(n\) samples uniformly at random. Then, we generate the choice \(c_l\), \(l \in \left[m\right]\) from the PL model given by Equation (1).

Retinopathy of Prematurity (ROP). The ROP dataset contains \(n=100\) vessel-segmented retina images. Experts are provided with two images and are asked to choose the image with higher severity of the ROP disease. Five experts independently label 5,941 image pairs; the resulting dataset contains \(m=29,705\) pairwise comparisons. Note that some pairs are labeled more than once by different experts. We employ two ROP datasets, ROP-num and ROP-img for shallow and DNN regression, respectively: Each sample in ROP-num has \(d=143\) extracted numerical features [Ataer-Cansızoğlu 2015], while each sample in ROP-img is an image with dimensions \(d=224 \times 224\).

Filter Aesthetic Comparison (FAC). The FAC dataset [Sun et al. 2017] contains 1,280 unfiltered images pertaining to eight different categories. A total of 22 different image filters are applied to each image. Labelers are provided with two filtered images and are asked to identify which image has better quality. We select \(n=1,000\) images within one category, as only the filtered image pairs that are within the same category are compared. The resulting dataset contains \(m=728\) pairwise comparisons. Moreover, for each image, we extract features via a state-of-the-art convolutional neural network architecture, namely GoogLeNet [Szegedy et al. 2015], with weights pre-trained on the ImageNet dataset [Deng et al. 2009]. We select \(d=50\) of these features by Principal Component Analysis [Jolliffe 1986].

SUSHI. The SUSHI Preference dataset [Kamishima et al. 2009] contains \(n=100\) sushi ingredients with \(d=18\) features. Each of the 5,000 customers independently ranks 10 ingredients according to her preferences. We select the rankings provided by 10 customers, where an ingredient is ranked higher if it precedes the other ingredients in a customer’s ranked list. We generate two datasets: triplet Sushi containing \(m=1,200\) rankings of \(|A_\ell |=3\) ingredients, and pairwise Sushi containing \(m=450\) pairwise comparisons.

International Conference on Learning Representations (ICLR). The ICLR dataset contains abstracts and reviewer ratings of 2,561 papers that are submitted to ICLR 2020 conference and are available on OpenReview website [Sun 2020]. We choose the top \(n=100\) papers, and extract \(d=768\) numerical features from each abstract using the Deep Bidirectional Transformers (BERT) [Devlin et al. 2019] architecture, pre-trained on the Books Corpus dataset [Zhu et al. 2015] and English Wikipedia. We normalize \(\mathbf {X}\) to have 0 mean and unit variance over samples \(\left[n\right]\). We generate all possible \(m=120,324 (2,248,524)\) \(K=3(4)\)-way rankings w.r.t. the relative order of the average reviewer ratings. We add noise to the resulting rankings following the same process as Movehub-Cost.

Movehub-Cost. The Movehub-Cost dataset contains the total ranking of 216 cities w.r.t. cost of living [Blitzer 2017]. Each city is associated with \(d=6\) numerical features, which are average costs for cappuccino, cinema, wine, gasoline, rent, and disposable income. We normalize \(\mathbf {X}\) to have 0 mean and unit variance over samples \(\left[n\right]\). We select \(n=50\) cities and generate all \(m=230,298 (2,118,756)\) \(K=4(5)\)-way rankings w.r.t. the relative order of the queried cities in the total ranking. To mimic the real-life noise introduced by human labeling, we apply the following post-processing to the resulting rankings: For each ranking, we sample a value uniformly at random in \([0,1]\). If the value is less than 0.1, we add noise to the ranking by a cyclic permutation of the ranked samples.

Movehub-Quality. The Movehub-Quality dataset contains total ranking of the same 216 cities as Movehub-Cost, this time w.r.t. quality of life. Each city is associated with \(d=5\) numerical features, including overall scores for purchase power, healthcare, pollution, quality of life, and crime. We normalize \(\mathbf {X}\) to have 0 mean and unit variance over samples \(\left[n\right]\). We select \(n=50\) cities and generate all \(m=230,298 (2,118,756)\) \(K=4(5)\)-way rankings w.r.t. the relative order of the queried cities in the total ranking. We add noise to the rankings following the same process as Movehub-Cost.

IMDB. The IMDB Movies Dataset contains IMDB ratings of 14,762 movies, each of which is associated with \(d=36\) numerical features [Leka 2016]. We normalize \(\mathbf {X}\) to have 0 mean and unit variance over samples \(\left[n\right]\). We select \(n=50\) movies and generate all possible \(m=85,583\) \(K=4\)-way rankings w.r.t. the relative order of the ratings of queried movies. We add noise to the resulting rankings following the same process as Movehub-Cost.

We partition each dataset into training and test sets in two ways. In rank partitioning, we partition the dataset w.r.t. \(\left[m\right]\), using \(90\%\) of the \(m\) observations for training, and the remaining \(10\%\) for testing. In sample partitioning, we partition \(\left[n\right]\), using \(90\%\) of the \(n\) samples for training, and the remaining \(10\%\) for testing. In this setting, observations containing samples from both training and validation/test sets are discarded. We perform cross validation on the resulting training sets. We use 10 folds for synthetic datasets, ROP-num, FAC, and Triplet Sushi, and 3 folds for the other datasets. For synthetic datasets, we also repeat experiments over five random generations.

To evaluate the convergence speed of each algorithm, we measure the elapsed time, including time spent in initialization, in seconds (Time). Moreover, we measure the prediction performance using Top-1 accuracy (Top-1 Acc.) and KT correlation [Kendall 1938] on validation and test sets (c.f. 4.7). For synthetic datasets, we also measure the quality of convergence by the norm of the difference between estimated and true PL scores (\(\triangle \mathbf {\pi }\)); lower values indicate better estimation. We report averages and standard deviations over folds for all algorithms except for the siamese network method: As training takes several hours, we execute only one fold.

We implement¹ seven competing algorithms. Four are feature methods, i.e., algorithms that regress PL scores from features: SR-l2 described in Section 3.3, SR-l2-log described in Section 3.4, sequential least-squares quadratic programming (SLSQP), that solves (16), and Newton on \({{\mathbf {w}}}\), that solves the convex problem (17) via Newton’s method. The remaining three are featureless methods, i.e., algorithms that learn the PL scores from the choice observations alone: ILSR described by Equation (7), the MM algorithm [Hunter 2004], and Newton on \(\mathbf {\theta }\) that solves Equation (3) via the reparametrization \({\pi }_i = e^{\theta _i}\), \(i \in \left[n\right]\) using Newton’s method on \(\mathbf {\theta }=[\theta _i]_{i\in \left[n\right]}\). We expand upon the implementation of all algorithms below.

SR-l2. We compute the stationary distribution at each iteration of ILSRX (c.f. Equation (29)) using the power method [Lei et al. 2016]. As the stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\) and \(\parallel \!\! \mathbf {\tilde{X}}\mathbf {w}^k - \mathbf {\tilde{X}}\mathbf {w}^{k-1} \!\!\parallel _2 \,\lt \, {{r_\mathtt {tol}}} \parallel \!\! \mathbf {\tilde{X}}\mathbf {w}^k \!\!\parallel _2\). We set the relative tolerance \(r_\mathtt {tol}=10^{-4}\) for all experiments. We use the same relative tolerance for the stopping criterion of the power method. We set \(\rho =1\) in our experiments, which is a standard choice in the ADMM literature [Boyd et al. 2011]. In our experiments, we consistently observe that Equation (26) is satisfied. That is why, we use \(c_j = \frac{- {\pi }_j \sigma _j}{\sum _{i\in \left[n\right]_{-}} {\pi }_i \sigma _i}\), \(j \in \left[n\right]_{+}\) to calculate the transition rates (27).

SR-l2-log. As the stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\) and \(\parallel \!\! e^{\mathbf {X}{{\mathbf {w}}^k}} - e^{\mathbf {X}{{\mathbf {w}}^{k-1}}} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! e^{\mathbf {X}{{\mathbf {w}}^k}} \!\!\parallel _2\), where exponentiation is applied elementwise.

SLSQP. We initialize SLSQP the same as SR-l2 (c.f. Section 3.3). As stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\), where \(\mathbf {\pi }^k = \mathbf {X}\mathbf {a}^k + \mathbf {1} {b}^k\), \(k \in \mathbb {N}\). Each iteration of SLSQP is \(O(\sum _{\ell \in \mathcal {D}}(|A_\ell | (d+1)) + (d+1)^2)\) for constructing the gradient of Equation (15) w.r.t. \(\mathbf {w}\) and updating \(\mathbf {w}\), respectively.

Newton on \({{\mathbf {w}}}\). We initialize Newton on \({{\mathbf {w}}}\) the same as SR-l2-log (c.f. Section 3.4). As stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\), where \(\mathbf {\pi }^k = [e^{ \mathbf {x}_i^T \mathbf {w}^k}]_{i\in \left[n\right]}\), \(k \in \mathbb {N}\). Each iteration of Newton on \({{\mathbf {w}}}\) is \(O(\sum _{\ell \in \mathcal {D}}(|A_\ell |\, d^2) + d^2)\) for constructing the Hessian of Equation (15) w.r.t. \(\mathbf {w}\) and updating \(\mathbf {w}\), respectively.

ILSR. We initialize ILSR with \({\mathbf {\pi }}^0 = \frac{1}{n}\mathbf {1}\). We compute the stationary distribution at each iteration of ILSR using the power method. As the stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\). Each iteration of ILSR is \(O(\sum _{\ell \in \mathcal {D}}(|A_\ell |) + n^2)\) for constructing the transition matrix \(\mathbf {\Lambda }(\mathbf {\pi })\) (c.f. Equation (6)) and finding the stationary distribution \(\mathbf {\pi }\), respectively.

MM. We initialize MM with \({\mathbf {\pi }}^0 = \frac{1}{n}\mathbf {1}\). As the stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\). Each iteration of MM is \(O(\sum _{\ell \in \mathcal {D}}(|A_\ell |))\).

Newton on \(\mathbf {\theta }\). We initialize Newton on \(\mathbf {\theta }\) with \({\mathbf {\theta }}^0 = [\theta _i^0]_{i\in \left[n\right]} = \mathbf {0}\). As stopping criterion, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\), where \({\mathbf {\pi }}^k_i=e^{\theta _i^k}\), \(i \in \left[n\right]\), \(k \in \mathbb {N}\). Each iteration of Newton on \(\mathbf {\theta }\) is \(O(\sum _{\ell \in \mathcal {D}}(|A_\ell |^2) + n^2)\) for constructing the Hessian of Equation (15) w.r.t. \(\mathbf {\theta }\) and updating \(\mathbf {\theta }\), respectively.

We implement² five competing algorithms. DSR with KL penalty (DSR-KL), DSR with \(\ell _2\) norm penalty (DSR-l2), and siamese network competitor regress scores via a DNN \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) with parameters \(\mathbf {w}\). SR-l2 and SR-KL use an affine regressor \(\tilde{\mathbf {\pi }}=\mathbf {X}\mathbf {a}+ \mathbf {1} {b}\) with parameters \(\mathbf {a}\) and \({b}\), and solve (19) via ADMM with \(\ell _2\) and KL penalties, respectively.

DSR-KL and DSR-l2 are explained in Section 3.5. We compute the stationary distribution at each iteration of ILSRX (c.f. Equation (29)) using the power method [Lei et al. 2016]. At each ADMM iteration, we fine-tune the DNN regressor \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) via Adam optimization [Kingma and Ba 2015] over Equation (22). We check the convergence of \(\mathcal {L}\) on the training set and KT correlation evaluations on the validation set (c.f. 4.7) as the stopping criterion for: (i) fine-tuning \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) at each ADMM iteration, and (ii) overall DSR-KL algorithm. We declare convergence when KT on validation set does not change for 5 iterations, for maximum total of 50 iterations, with relative tolerance \(r_\mathtt {tol}=10^{-4}\). We use the same relative tolerance for the stopping criterion of the power method.

To aid convergence in practice [Boyd et al. 2011], we update the dual variable at each ADMM iteration with a multiplicative smoothing parameter \(\gamma ^k=1/k\). We also adapt the penalty parameter aswhere \(\tau =2\) and \(\beta =10\).

The siamese network competitor minimizes Equation (15) w.r.t. \(\mathbf {w}\) via SGD. We train a siamese network architecture with base network \(\tilde{{\pi }}\) on observations \(\mathcal {D}\) via Adam optimization [Kingma and Ba 2015] over Equation (15). We employ the same convergence criterion and experiment setup as DSR-KL and DSR-l2 to train and optimize the siamese network.

Finally, we implement two spectral algorithms that regress Plackett-scores via an affine model, i.e., \(\tilde{\mathbf {\pi }}=\mathbf {X}\mathbf {a}+ \mathbf {1}{b}\): SR-l2, and its variant SR-KL that uses KL proximal penalty instead of \(\ell _2\) norm penalty for ADMM. As the stopping criterion for both algorithms, we use \(\parallel \!\! {\mathbf {\pi }}^k- {\mathbf {\pi }}^{k-1} \!\!\parallel _2 \,\lt \, {r_\mathtt {tol}} \parallel \!\! {\mathbf {\pi }}^k \!\!\parallel _2\) and \(\parallel \!\! (\mathbf {X}\mathbf {a}^k + \mathbf {1}{b}^k) - (\mathbf {X}\mathbf {a}^{k-1} +\mathbf {1} {b}^{k-1} \!\!\parallel _2 \,\lt \, {{r_\mathtt {tol}}} \parallel \!\! \mathbf {X}\mathbf {a}^k + \mathbf {1}{b}^k \!\!\parallel _2\). We set the ADMM penalty parameter as \(\rho =1\).

To evaluate each method on ROP-img, we choose \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) as a state-of-the-art convolutional neural network architecture, namely GoogLeNet [Szegedy et al. 2015], followed by a fully connected output layer comprising a single neuron with sigmoid activation. We initialize the convolutional layers with weights pre-trained on the ImageNet dataset [Deng et al. 2009]. For other datasets, we design \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\) as a fully connected architecture with relu activation for hidden layers and an output layer comprising a single neuron with sigmoid activation. We add \(\ell _2\) regularizers to all layers. For each configuration of (i) \(\ell _2\) regularization parameter varying in \([2\times 10^{-5},2\times 10^{-2}]\), (ii) learning rate varying in \([10^{-4},10^{-2}]\), and (iii) number of layers varying in \([1,10]\) for fully connected \(\tilde{{\pi }}(\mathbf {X}; \mathbf {w})\), we run each method until convergence (c.f. 4.4 for criteria). We determine the best set of hyperparameters w.r.t. the prediction performance via cross validation (c.f. 4.2).

In Tables 3–5, we measure timing on an Intel Xeon CPU E5-2680v2 2.8 GHz with 128 GB RAM. Particularly for experiments on larger synthetic datasets (c.f. Figures 3–4), we use an Intel Xeon CPU E5-2680v4 2.4 GHz with 500 GB RAM. For all DNN regression experiments, we employ NVIDIA GPUs with Intel Gold [email protected] GHz CPUs and 128 GB RAM.

We measure the prediction performance by Top-1 accuracy (Top-1 Acc.) and KT correlation on the test set. Let the test set be \(\mathcal {D}_{\mathtt {rank}} = \lbrace (\alpha ^\ell , A_\ell) \,|\, \ell \in \lbrace 1,\ldots ,m_{\mathtt {test}}\rbrace \rbrace\), where \(\alpha ^\ell = \alpha _{1}^\ell \succ \alpha _{2}^\ell \succ \cdots \succ \alpha ^\ell _{K}\) is an ordered sequence of the samples in \(A_\ell\). Given \(A_\ell\), we predict the \(\ell\)th choice as \(\hat{c}_\ell = \text{arg,max}_{i \in A_\ell } {\pi }_i\). We calculate the Top-1 accuracy (Top-1 Acc.) asWe also predict the ranking as \(\hat{\alpha }^\ell = \text{arg,sort}[{\pi }_i]_{i \in A_\ell }\), i.e., sequence of the samples in \(A_\ell\) ordered w.r.t. their estimated scores. We calculate KT correlation [Kendall 1938] as a measure of the correlation between each true ranking \(\alpha ^\ell\) and predicted ranking \(\hat{\alpha }^\ell\), \(\ell \in \lbrace 1,\ldots ,m_{\mathtt {test}}\rbrace\). For observation \(\ell\), let \(T_\ell = \sum _{t=1}^{K} \sum _{s=1}^{K} \mathbb {1}(\hat{\alpha }_t^\ell \succ \hat{\alpha }_s^\ell \wedge \alpha _t^\ell \succ \alpha _s^\ell)\) be the number correctly predicted ranking positions, and \(F_\ell = \sum _{t=1}^{K} \sum _{s=1}^{K} \mathbb {1}(\hat{\alpha }_t^\ell \succ \hat{\alpha }_s^\ell \wedge \alpha _s^\ell \succ \alpha _t^\ell)\) be the number incorrectly predicted ranking positions. Then, KT is computed bywhere \(\left({{K}\atop {2}}\right)\) is the number of sample pairs.

Sample partitioning. We begin with the experiments on sample partitioning, in which we partition samples \(\left[n\right]\) into training and test sets. Table 3 shows the evaluations of all algorithms trained on a synthetic dataset with \(n=1,000\) samples, \(d=100\) features, \(m=1,000\) observations, and query size \(|A_\ell |=2\). SR-l2 and SR-l2-log converge 4–27 times faster than other feature methods, i.e., Newton on \({{\mathbf {w}}}\) and SLSQP. Recall that in sample partitioning, training observations contain only training samples: Test samples do not participate in any of the training observations. Thus, featureless methods ILSR, MM, and Newton on \(\mathbf {\theta }\) are no better than random predictors, with 0.5 Top-1 Acc. and 0.0 KT. By regressing the PL scores from features, SR-l2 and SR-l2-log significantly outperform the predictions of ILSR, MM, and Newton on \(\mathbf {\theta }\), by 16%–33% Top-1 Acc. and 16%– 30% KT.

Real datasets. We observe an equally significant speed gain on real datasets; Table 4 shows the evaluations on four real datasets partitioned w.r.t. sample partitioning. SR-l2 and SR-l2-log are 3–18 times faster than Newton on \({{\mathbf {w}}}\) and SLSQP. This speed gain is fundamentally due to the smaller per iteration complexity of SR-l2 and SR-l2-log. For instance, compared to Newton on \({{\mathbf {w}}}\), SR-l2-log requires about 2 times more iterations, but still converges 3 times faster than Newton on \({{\mathbf {w}}}\) on FAC. Moreover, while significantly decreasing the convergence time, SR-l2 or SR-l2-log consistently attain similar prediction performance to Newton on \({{\mathbf {w}}}\) and SLSQP, except for Sushi, for which they perform slightly worse (by \(13\%\)–\(20\%\) Top-1 Acc.), though the convergence time dividends are striking in comparison (\(78\%\)–\(95\%\)).

Aligned with the prediction performance on synthetic datasets, featureless methods ILSR, MM, and Newton on \(\mathbf {\theta }\) can only attain 0.5 Top-1 Acc. and 0.0 KT. By regressing the PL scores from features, SR-l2 and SR-l2-log significantly outperform the predictions of ILSR, MM, and Newton on \(\mathbf {\theta }\), by \(3\%\)–\(31\%\) Top-1 Acc. and \(5\%\)–\(78\%\) KT.

Rank partitioning. A sample can appear in both training and test observations in rank partitioning. Hence, featureless methods ILSR, MM, and Newton on \(\mathbf {\theta }\) should fare better than in sample partitioning. Nonetheless, in Table 3, as \(n=1,000\) is larger than \(d=100\), there are more scores to learn than parameters. As a result, feature methods are advantageous for good predictions compared to featureless methods. Particularly, SR-l2 and SR-l2-log outperform the predictions of ILSR, MM, and Newton on \(\mathbf {\theta }\) in rank partitioning on Table 3, by \(13\%\) Top-1 Acc. and \(9\%\) KT. The relative performance of feature vs. featureless methods is governed by the relationship among \(n\), \(d\), and \(m\). We do not include Newton on \(\mathbf {\theta }\) and MM in this analysis, as they are too slow.

Impact of \(d\). To assess the impact of number of parameters, we fix \(n=1,000\), \(m=250\), \(|A_\ell |=2\) and generate synthetic datasets with \(d\in \lbrace 10,100,1000,10000\rbrace\). Figure 3(a) shows the Time and Top-1 Acc. of SR-l2, SR-l2-log, ILSR, and Newton on \({{\mathbf {w}}}\). As \(m=250\) observations are not enough to learn \(n=1,000\) scores, SR-l2 leads to significantly better Top-1 Acc. compared to ILSR. When \(d=10\), SR-l2 and SR-l2-log outperform ILSR by \(18\%\)–\(28\%\) Top-1 Acc. Moreover, SR-l2 and SR-l2-log are consistently faster than Newton on \({{\mathbf {w}}}\), for all \(d\gt 100\). Particularly, for \(d=10,000\), SR-l2 and SR-l2-log converge 42–579 times faster than Newton on \({{\mathbf {w}}}\). Interestingly, the convergence time of SR-l2-log can even decrease with increasing \(d\). This is because the number of iterations until convergence decreases. While significantly decreasing the convergence time, SR-l2 consistently attains better Top-1 Acc. than Newton on \({{\mathbf {w}}}\), up to \(8\%\) for \(d=100\). Impact of \(n\). To assess the impact of number of samples, we fix \(d=100\), \(m=250\), \(|A_\ell |=2\) and generate synthetic datasets with \(n\in \lbrace 50,100,1000,10000,100000\rbrace\); Figure 3(b) shows evaluations on the resulting datasets. For \(n\gt d=100\), i.e., when there are more scores to learn than parameters, SR-l2 leads to significantly better Top-1 Acc. compared to ILSR. Particularly, for \(n=100,000\), SR-l2 outperforms ILSR by \(25\%\) Top-1 Acc. This confirms that, especially when the number of observations \(m\) is not sufficient to learn \(n\) scores, exploiting the features associated with the samples is crucial in attaining good prediction performance. As expected, convergence time of Newton on \({{\mathbf {w}}}\) is not significantly affected by \(n\). Despite this, SR-l2 and SR-l2-log are faster than Newton on \({{\mathbf {w}}}\) for all \(n\lt 1,\!000\). Particularly, for \(n=50\), SR-l2 and SR-l2-log converge 2–75 times faster than Newton on \({{\mathbf {w}}}\). While decreasing the convergence time, SR-l2 consistently attains better Top-1 Acc. than Newton on \({{\mathbf {w}}}\), up to \(19\%\) for \(n=100,\!000\).

Impact of \(m\). Figure 4 shows the convergence time (Time) and top-1 test accuracy (Top-1 Acc.) of SR-l2, SR-l2-log, ILSR, and Newton on \({{\mathbf {w}}}\) when trained on synthetic datasets with number of observations \(m\in \lbrace 10,100,1000,10000,100000\rbrace\). Observations are partitioned w.r.t. rank partitioning, where number of samples is \(n=1,\!000\), number of parameters is \(d=100\), and size of each query is \(|A_\ell |=2\). As \(n\gt d\), SR-l2 benefits from being able to regress \(n\) scores from a smaller number of \(d\) parameters and leads to significantly better Top-1 Acc compared to ILSR in Figure 4. Especially when \(m\) is not enough to learn \(n=1,000\) scores, but to learn \(d=100\) parameters, SR-l2 gains the most performance advantage over ILSR, up to \(13\%\) Top-1 Acc. Moreover, SR-l2 and SR-l2-log are consistently faster than Newton on \({{\mathbf {w}}}\), for all number of observations \(m\gt 100\). Particularly, for \(m=100,000\), SR-l2 and SR-l2-log converge 4–60 times faster than Newton on \({{\mathbf {w}}}\).

Real datasets. Performance agrees with observations above regarding the dependence on \(n\) and \(d\). For datasets where \(n\gt m\gt d\), e.g., FAC, SR-l2, and SR-l2-log significantly outperform the predictions of ILSR, by \(10\%\) Top-1 Acc. and \(25\%\) KT. For datasets where \(m\) is much larger than \(n\) (c.f. Table 2), feature methods lead to similar prediction performance to each other and slightly lower performance than ILSR, MM, and Newton on \(\mathbf {\theta }\). Overall, SR-l2 and SR-l2-log consistently converge faster than Newton on \({{\mathbf {w}}}\) and SLSQP, by 3–27 times across all real datasets. Table 5 also confirms the observations from Figure 3(a) and (b) regarding the effect of \(n\) and \(d\). For FAC, \(n=1000 \gt m=728 \gt d=50\), i.e., there are enough observations to learn parameters, but not scores. As a result, SR-l2 and SR-l2-log significantly outperform the predictions of ILSR, by \(10\%\) Top-1 Acc. and \(25\%\) KT. For the other real datasets, \(m\) is much larger than \(n\) (c.f. Table 2). Thus, feature methods lead to similar prediction performance to each other and lower performance than ILSR, MM, and Newton on \(\mathbf {\theta }\).

Training Time vs. Prediction Performance. Figure 5(a) and (b) show the training time of DSR-KL and siamese network vs. Top-1 Acc. and KT on test sets of Movehub-Cost-4, Movehub-Quality-4, and IMDB-4 datasets, partitioned with rank partitioning. For all datasets, DSR-KL results lie on the top left region of both Top-1 Acc. and KT plots: DSR-KL consistently leads to much faster training and better predictions than the siamese counterpart.

Columns 3–5 in Table 6 show the training time and test set prediction performance of DSR-KL, DSR-l2, siamese network, SR-l2, and SR-KL trained on datasets partitioned with rank partition- ing. DSR-KL and DSR-l2 are 1.5–142 times faster than siamese network over all datasets. Moreover, DSR-KL consistently attains equivalent or better prediction performance than both siamese network and DSR-l2 w.r.t. both Top-1 Acc. and KT. DSR-KL leads to particularly better performance in ranking predictions, by up to \(6\%\) higher KT than siamese and \(25\%\) higher KT than DSR-l2 on IMDB-4.

Deeper regression methods consistently outperform the predictions of shallow regression methods. Particularly, our deep spectral algorithms DSR-KL and DSR-l2 lead to significantly better predictions than the shallow versions SR-l2 and SR-KL, up to \(38\%\) Top-1 Acc. and \(41\%\) KT on IMDB-4. The training times of DSR-KL and DSR-l2 are also not noticeably larger than the ones of shallow versions; SR-KL converges even slower than DSR-KL, by up to 12 times on Movehub-Quality-5. Unlike SR-l2 that solves a least-squares problem with closed form solution, parameter update step of SR-KL (c.f. 48) requires an iterative optimization at each ADMM iteration.

Columns 6–8 in Table 6 show the training time and test set prediction performance of DSR-KL, DSR-l2, siamese network, SR-l2, and SR-KL trained on all datasets, partitioned with sample partitioning. Note that this is the setting when regressing scores from features is essential, as training samples cannot participate in any observations in validation or test sets. Agreeing with the speed gain in rank partitioning, DSR-KL and DSR-l2 are 8–175 times faster than the siamese network over all datasets. Moreover, DSR-KL leads to particularly better performance in maximal-choice predictions, by up to \(26\%\) higher Top-1 Acc. than siamese on IMDB-4 and \(56\%\) higher Top-1 Acc. than DSR-l2 on Movehub-Cost-4. Finally, deeper regression methods DSR-KL, DSR-l2, and siamese network, outperform the predictions of shallow counterparts SR-l2 and SR-KL, by up to \(39\%\) Top-1 Acc. and \(11\%\) KT on Movehub-cost-5 and ICLR-3, respectively.

Impact of \(K\). Figures 1 and 5(c) show training time and prediction performances of DSR-KL, DSR-l2, and siamese network vs. query size (\(K\)) on Movehub-Cost and Movehub-Quality datasets partitioned w.r.t. rank partitioning. The speed gain of DSR-KL and DSR-l2 over siamese reach up to 43 times as \(K\) increases, as each epoch of siamese grows exponentially with \(K\). Moreover, agreeing with Table 6, DSR-KL consistently outperforms the predictions of DSR-l2.

Details on Convergence. Figure 6 shows the log-likelihood \(-\mathcal {L}\) on training and Top-1 Acc. on validation sets of Movehub-Cost-4, Movehub-Quality-4, and IMDB-4 datasets, partitioned w.r.t. rank partitioning. Each point for DSR-KL and DSR-l2 correspond to an overall iteration of ADMM (c.f. (21)), while each point for siamese corresponds to a training epoch. DSR-KL consistently attains higher log-likelihood and better validation performance than both siamese network and DSR-l2, while converging faster.

Comparison to Naïve Approach. A naïve spectral algorithm can be constructed by a single iteration of the primal steps in (21): (i) solving (21a) to learn \(\mathbf {\pi }\) via repeated iterations of (29), and (ii) given \(\mathbf {\pi }\), solving (21b) by training the DNN regressor \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\) via SGD over Equation (22). Intuitively, this ignores/does not exploit the fact that samples with similar features ought to have similar scores. We denote this naive approach as One-iter.-DSR-KL. Unlike One-iter.-DSR-KL, our algorithm DSR-KL has the capability of repeatedly adapting both \(\mathbf {\pi }\) and \(\mathbf {w}\) by solving (19) via ADMM. This advantage is illustrated in Figure 6(a) and (b), where not only \(\tilde{\mathbf {\pi }}(\mathbf {X}; \mathbf {w})\), but also \(\mathbf {\pi }\) are adjusted until convergence. Moreover, Figure 6(d) shows the corresponding test set predictions of DSR-KL vs. One-iter.-DSR-KL on Movehub-Cost-4 and Movehub-Quality-4; DSR-KL outperforms One-iter.-DSR-KL by 13–21% Top-1 Acc. and 6–15% KT.