Preference Completion: Large-scale Collaborative Ranking from Pairwise Comparisons

Introduction

This paper considers the following recommendation system problem: given a set of items, a set of users, and non-numerical pairwise comparison data, find the underlying preference ordering of the users. In particular, we are interested in the setting where data is of the form “user $i$ preferes item $j$ over item $k$”, for different ordered user-item-item triples $i,j,k$. Pairwise preference data is wide-spread; indeed, almost any setting where a user is presented with a menu of options – and chooses one of them – can be considered to be providing a pairwise preference between the chosen item and every other item that is presented.

Crucially, we are interested in the collaborative filtering setting, where (a) on the one hand the number of such pairwise preferences we have for any one user is woefully insufficient to infer anything for that user in isolation; and (b) on the other hand, we aim for personalization, i.e. for every user to possibly have different inferred preferences from every other. To reconcile these two requirements, our method relates the preferences of users to each other via a low-rank matrix, which we (implicitly) assume governs the observed preferences. Essentially, we fit a low-rank users $\times$ items score matrix $X$ to pairwise comparison data by trying to ensure that $X_{ij}-X_{ik}$ is positive when user $i$ prefers item $j$ to item $k$.

We present two algorithms to infer the score matrix $X$ from training data; once inferred, this can be used for predicting future preferences. While there has been some recent work on fitting low-rank score matrices to pairwise preference data (which we review and compare to below), in this paper we present the following two contributions: (a) A statistical analysis for the convex relaxation: we bound the generalization error of the solution to our convex program. Essentially, we show that the minimizer of the empirical loss also almost minimizes the true expected loss. We also give a lower bound showing that our error rate is sharp up to logarithmic factors. (b) A large-scale non-convex implementation: We provide a non-convex algorithm that we call Alternating Support Vector Machine (AltSVM). This non-convex algorithm is more practical than the convex program in a large-scale setting; it explicitly parameterizes the low-rank matrix in factored form and minimizes the hinge loss. Crucially, each step in this algorithm can be formulated as a standard SVM that updates one of the two factors; the algorithm proceeds by alternating updates to both factors. We apply a stochastic version of dual coordinate descent with lock-free parallelization. This exploits the problem structure and ensures it parallelizes well. We show that our algorithm outperforms several existing collaborative ranking algorithms in both speed and prediction accuracy, and it achieves significant speedups as the number of cores increases.

1 Related Work

Ranking/learning preferences is a classical problem that has been considered in a large amount of work. There are many different settings for this problem, which we discuss below.

The main problem in this community has been to estimate a ranking function from given feature vectors and relevance scores. Depending on its application, a feature vector may correpond to a user-item pair or a single item. While there have been algorithms that use pairwise comparisons of the training samples, our setting is different in that our data consists only of pairwise comparisons. We refer the reader to the survey .

In a single-user model, we are asked to learn a single ranking given pairwise comparisons. Jamieson & Nowak and Ailon consider an active query model with noiseless responses; Jamieson & Nowak give an algorithm for exactly recovering the true ranking under a low-rank assumption similar to ours, while Ailon approximately recovers the true ranking without such an assumption. Wauthier et al. and Negahban et al. learn a ranking from noisy pairwise comparisions; Negahban et al. consider a Bradley-Terry-Luce model similar to ours and attempt to learn an underlying score vector, while Wauthier et al. get by without structure assumptions, but only attempt to learn the ranking itself. Hajek et al. considered a problem to learn a single ranking given a more generalized partial rankings from the Plackett-Luce model and provided a minimax-optimal algorithm.

Given multiple users with different rankings, one could of course attempt to learn their rankings by simply applying an algorithm from the previous section to each user individually. However, it is more efficient – both statistically and computationally – to postulate some global structure and use it to relate the many users’ rankings. This is the same idea that has been applied so successfully in collaborative filtering. Rendle et al. and Liu et al. were the first to take this approach. They modeled the observations as coming from a BTL model with low-rank structure (i.e., very similar to our model) and gave algorithms for learning the model parameters. Yi et al. took a purely optimization-based approach. Rather than assuming a probabilistic model, they minimized a convex objective using the hinge loss on a low-rank matrix. In a slightly different model, Hu et al. and Shi et al. consider the problem of learning from latent feedback. Recently, Lu & Negahban analyzed an algorithm which is very similar to ours for the Bradley-Terry-Luce model independently from our work.

Instead of moving to pairwise comparisons, some work has suggested avoiding the difficulties of numerical ratings by instead asking users to give 1-bit ratings to items; that is, each user only indicates whether they like or dislike an item. In this setting, the work of Davenport et al. is most closely related to ours, in that they assume an underlying low-rank structure and give an algorithm based on convex optimization. Also, our theoretical analysis owes a lot to their work. Xu et al. consider a slightly different goal: rather than attempting to recover the preferences of each user, they try to cluster similar users and similar items together. Yun et al. proposed an optimization problem motivated from robust binary classification and used stochastic gradient descent to solve the problem in a large-scale setting.

The goal in this setting is the same as ours, except that the data is in the form of numerical ratings instead of pairwise comparisons. Weimer et al. attempted to directly optimize Normalized Discounted Cumulative Gain (NDCG), a widely used performance measure for ranking problems. Balakrishnan & Chopra , and Volkovs & Zemel converted this problem into a learning-to-rank problem and solved it using the existing algorithms. While these works considered the low-rank matrix model, different models are proposed by Weston et al. and Lee et al. . Weston et al. proposed a tensor model to rank items for different queries and users, and proposed a weighted sum of low-rank matrix models.

Empirical Risk Minimization (ERM)

Let us first formulate the problem mathematically. The task is to estimate rankings of multiple users on multiple items. We denote the numbers of users by $d_{1}$, and the number of items by $d_{2}$. We are given a set of triples $\Omega\subset[d_{1}]\times[d_{2}]\times[d_{2}]$, where the preference of user $i$ between items $j$ and $k$ is observed if $(i,j,k)\in\Omega$. The observed comparison is then given by $\{Y_{ijk}\in\{1,-1\}:(i,j,k)\in\Omega\}$ where $Y_{ijk}=1$ if user $i$ prefers item $j$ over item $k$, and $Y_{ijk}=-1$ otherwise. Let $\Omega_{i}=\{(j,k):(i,j,k)\in\Omega\}$ denote the set of item pairs that user $i$ has compared.

We propose (as have others) that $X$ is low-rank or close to low-rank, the intuition being that each user bases their preferences on a small set of features that are common among all the items. Then the empirical risk minimization (ERM) framework can naturally be formulated as

where $\mathcal{L}(\cdot)$ is a monotonically non-increasing loss function which induces $X_{ij}>X_{ik}$ if $Y_{ijk}=1$, and $X_{ij}<X_{ik}$ otherwise. (e.g., hinge loss, logistic regression loss, etc.)

Solving (1) is NP-hard because of the rank constraint. As a first alternative, we propose a straightforward convex relaxation.

Convex Relaxation

Our first method is the convex relaxation of (1), which involves a nuclear norm constraint.

Here, for any matrix $X$, the nuclear/trace norm $\|X\|_{*}$ denotes the sum of its singular values; it is a well-recognized convex surrogate for low-rank structure (most famously in matrix completion).

The only parameter of this algorithm is $\lambda$, which governs the trade-off between better optimizing the likelihood of the observed data, and the strictness in imposing approximate low-rank structure. Since we motivated our algorithm with the assumption that $X$ has low rank, we should point out how our algorithm’s parameter $\lambda$ compares to the rank: note that if $X$ is a $d_{1}\times d_{2}$ rank-$r$ matrix whose largest absolute entry is bounded by $C$ then $\|X\|_{*}\leq\sqrt{r}\|X\|_{F}\leq C\sqrt{rd_{1}d_{2}}$. In other words, $\lambda$ is a parameter that takes into account both the rank of $X$ and the size of its elements, and it is roughly proportional to the rank.

We analyze (2) by assuming a standard model for pairwise comparisons. Then we provide a statistical guarantee of the method under the model.

Assume that each user-item-item triple $(i,j,k)$ independently belongs to $\Omega$ with probability $p_{i,j,k}$, and let $m=\sum_{i,j,k}p_{i,j,k}$ be the expected size of $\Omega$. We will assume that the $p_{i,j,k}$ are approximately balanced in the sense that no user-item pair is observed too frequently:

There is a constant $\kappa>0$ such that for every $i,j$,

Note that if $\kappa=1$ in Assumption 3.1 then the $p_{i,j,k}$ are all equal, meaning that each user-item-item triple has an equal chance to be observed.

Our main upper bound shows that if $m$ is sufficiently large then our algorithm finds a solution with almost minimal risk. Given a loss function $\mathcal{L}$, define the expected risk of $X$ by

where the expectation is with respect to the distribution parametrized by the true parameters $X^{*}$.

We recall that the parameter $\lambda$ is related to rank in that if $X$ is a $d_{1}\times d_{2}$ rank-$r$ matrix whose largest absolute entry is bounded by $C$ then $\|X\|_{*}\leq\sqrt{r}\|X\|_{F}\leq C\sqrt{rd_{1}d_{2}}$. In other words, $\lambda$ is a parameter that takes into account both the rank of $X^{*}$ and the size of its elements, and it is roughly proportional to the rank. In particular, Theorem 3.1 shows that once we observe $m\sim r(d_{1}+d_{2})\log^{2}(d_{1}+d_{2})$ pairwise comparisons, then we can accurately estimate the probability of any user preferring any item over any other. In other words, we need to observe about $r(1+d_{2}/d_{1})\log^{2}(d_{1}+d_{2})$ comparisons per user, which is substantially less than the $rd_{2}\log(d_{2})$ comparisons that we would have required if each user were modelled in isolation. Moreover, our lower bound (below) shows that at least $r(1+d_{2}/d_{1})$ comparisons per user are required, which is only a logarithmic factor from the upper bound.

where $c>0$ is a constant depending only on $\mathcal{L}$.

Together, Theorems 3.1 and 3.2 show that (up to logarithmic factors) if $X^{*}$ has rank $r$ then about $r(1+d_{2}/d_{1})$ comparisons per user are necessary and sufficient for learning the users’ preferences.

By specializing the loss function $\mathcal{L}$, Theorem 3.1 has a simple corollary for maximum-likelihood estimation of $X^{*}$. Recall that if $\mu$ and $\nu$ are two probability distributions on a finite set $S$ the the Kullback-Leibler divergence between them is

under the convention that $0\log 0=0$. We recall that although $D(\cdot\|\cdot)$ is not a metric it is always non-negative, and that $D(\mu\|\nu)=0$ implies $\mu=\nu$.

Note that the loss function in Corollary 3.3 is exactly the negative logarithm of the logistic function, and so $\hat{X}$ in Corollary 3.3 is the maximum-likelihood estimate for $X^{*}$. Thus, Corollary 3.3 shows that the distribution induced by the maximum-likelihood estimator is close to the true distribution in Kullback-Leibler divergence.

Large-scale Non-convex Implementation

While the convex relaxation is statistically near optimal, it is not ideal for large-scale datasets because it requires the solution of a convex program with $d_{1}\times d_{2}$ variables. In this section we develop a non-convex variant which both scales and parallelizes very well, and has better empirical performance as compared to several existing empirical baseline methods.

Our approach is based on the following steps:

We solve the non-convex problem by alternating between updating $U$ while keeping $V$ fixed, and vice versa. With the hinge loss (which we found works best in experiments), each of these becomes an SVM problem - hence we call our algorithm AltSVM.

The problem is of course not symmetric in $U$ and $V$ because users rank items but not vice versa. For the $U$ update, each user vector naturally decouples and can be done in parallel (and in fact just reduces to the case of rankSVM ).

For the $V$ update, we show that this can also be made into an SVM problem; however it involves coupling of all item vectors, and all user ratings. We employ several tricks (detailed below) to speed up and effectively parallelize this step.

where we replace the nuclear norm regularizer using the property $\|X\|_{*}=\min_{X=UV^{\top}}\frac{1}{2}(\|U\|_{F}^{2}+\|V\|_{F}^{2})$ . $u_{i}^{\top}$ and $v_{i}^{\top}$ denote the $i$th rows of $U$ and $V$, respectively. While this is a non-convex algorithm for which it is hard to find the global optimum, it is computationally more efficient since only $(d_{1}+d_{2})r$ variables are involved. We propose to use L2 hinge loss, i.e., $\mathcal{L}(x)=\max(0,1-x)^{2}$.

In the alternating minimization of (3), the subproblem for $U$ is to solve

while $V$ is fixed. This can be decomposed into $n$ independent problems for $u_{i}$’s where each solves for

This part is in general a small-scale problem as the dimension is $r$, and the sample size is $|\Omega_{i}|$ for each user $i$.

On the other hand, solving for $V$ with fixed $U$ can be written as

We note that the feature matrices $\{A^{(i,j,k)}:(i,j,k)\in\Omega\}$ are highly sparse since in each feature matrix only $2r$ out of the $d_{2}r$ elements are nonzero. This motivates us to apply the stochastic dual coordinate descent algorithm , which not only converges fast but also takes advantages of feature sparsity in linear SVMs. Each coordinate descent step takes $O(r)$ computation, and iterations over $|\Omega|$ coordinates provide linear convergence .

where $\mathcal{L}^{*}(z)$ is the convex conjugate of $\mathcal{L}$. At each coordinate descent step for $\alpha_{ijk}$, we find the value of $\alpha_{ijk}$ minimizing (7) while all the other variables are fixed. If we maintain $u_{i}=\sum_{(j,k)\in\Omega_{i}}\alpha_{ijk}Y_{ijk}(v_{j}-v_{k})$, then the coordinate descent step is simply to find $\delta^{*}$ minimizing

and update $\alpha_{ijk}\leftarrow\alpha_{ijk}+\delta^{*}$.

where $\beta$ is the dual vector for the subproblem (6). Similarly to $\alpha_{ijk}$, the coordinate descent step for $\beta_{ijk}$ is to replace $\beta_{ijk}$ by $\beta_{ijk}+\delta^{*}$ where $\delta^{*}$ minimizes

and maintain $V=\sum_{(i,j,k)\in\Omega}\beta_{ijk}Y_{ijk}A^{(i,j,k)}$.

The detailed description of AltSVM is presented in Algorithm 1. In each subproblem, we run the stochastic dual coordinate descent, in which a pairwise comparison $(i,j,k)\in\Omega$ is chosen uniformly at random, and the dual coordinate descent for $\alpha_{ijk}$ or $\beta_{ijk}$ is computed. We note that each coordinate descent step takes the same $O(r)$ computational cost in both subproblems, while the subproblem sizes are much different.

For each subproblem, we parallelize the stochastic dual coordinate descent algorithm asynchronously without locking. Given $T$ processors, each processor randomly sample a triple $(i,j,k)\in\Omega$ and update the corresponding dual variable and the user or item vectors. We note that this update is for a sparse subset of the parameters. In the user part, a coordinate descent step for one sample updates only $r$ out of the $rd_{1}$ variables. In the item part, one coordinate descent step for a sample update only $2r$ out of the $rd_{2}$ variables. This motivates us not to lock the variables when updated, so that we ignore the conflicts. This lock-free parallelism is shown to be effective in for stochastic gradient descent (SGD) on the sum of sparse functions. Moreover, in , it is also shown that the stochastic dual coordinate descent scales well without locking. We implemented the algorithm using the OpenMP framework. In our implementations, we also parallelized steps 3 and 13 of Algorithm 1. We show in the next section that our proposed algorithm scales up favorably.

2 Remark on the implementation

In Algorithm 1, the subproblem for $V$ comes first, and then it solves for the user vectors $U$. We empirically observed that this order gives better convergence on practical datasets. We also note that each subproblem reuses the dual variables in the previous outer iteration. When almost converged, the features ($V$ for solving $U$, and $U$ for solving $V$) do not change too much. By reusing the dual variables in the previous iteration we can start with a feasible solution close to the optimum.

Experimental results

We used the MovieLens 100k dataset, which contains 100,000 ratings given by 943 users on 1682 movies. The ratings are given as integers from one to five, but we converted them into preference data by declaring that a user preferred one movie to another if they gave it a higher rating (if two movies received the same rating, we treated it as though the user did not provide a preference). Then we held out $20\%$ of the data as a test set.

We compared our algorithm to the following two:

Bayesian Personalized Ranking (BPR) : This algorithm is based on a similar model to ours, but a different optimization procedure (essentially, a variant of stochastic gradient descent).

Matrix completion from pairwise differences : A standard matrix completion algorithm that observes – for various triples $(i,j,k)\in\Omega$ – the difference between user $i$’s ratings for item $j$ and item $k$. Note that this algorithm has an advantage over (2) because it sees the magnitude of this difference instead of only its sign. Nevertheless, the matrix completion algorithm does not perform any better than (2). A similar phenomenon was also observed in .

We evaluate our performance by computing the proportion of pairwise comparisons in the test set $\mathcal{T}$ for which we correctly infer the user’s preference.

This is similar to the AUC statistic measured by Rendle et al. , and if the data were fully observed then it would measure Kendall’s distance between each user’s true preferences and the learned ones. However, our main reason for choosing this measure of performance is that, as an average accuracy over all pairwise comparisions, it resembles the quantity that we study in our theoretical bounds.

Unsurprisingly, we were more accurate at correctly inferring strong preferences; therefore, we have also shown the accuracy obtained by only measuring performance on pairs whose rankings differ by two or more. Both the methods we considered do measurably better at predicting these orderings.

2 Large-scale experiments on rating data

Now we demonstrate that our algorithm performs well as a collaborative ranking method on rating data. We used the datasets specified in Table 1. Given a training set of ratings for each user, our algorithm will only use non-tying pairwise comparisons from the set, while other competing algorithms use the ratings themselves. Hence, they have more information than ours. The competing algorithms are those with publicly available codes provided by the authors.

CofiRank http://www.cofirank.org, The dimension and the regularization parameter are set as suggested in the paper. For the rest of the parameters, we left them as provided. This algorithm uses alternating minimization to directly optimize NDCG.

Local Collaborative Ranking (LCR) http://prea.gatech.edu, We run the code with each of the 48 sets of loss function and parameters given in the main code, and the best result is reported. We could not run this algorithm on the Netflix dataset due to time constraint. : The main idea is to predict preferences from the weighted sum of multiple low-rank matrices model.

RobiRank https://bitbucket.org/d_ijk_stra/robirank, We used the part for collaborative ranking from binary relevence score. We left the parameter settings as provide with the implementation. : This algorithm uses stochastic gradient descent to optimize the loss function motivated from robust binary classification.

Global Ranking : To see the effect of personalized ranking, we compare the results with a global ranking of the items. We fixed $U$ to all ones and solved for $V$.

The algorithms are compared in terms of two standard performance measures of ranking, which are NDCG and Precision@$K$. NDCG@$K$ is the ranking measure for numerical ratings. NDCG@$K$ for user $i$ is defined as

and $\pi_{u}(k)$ is the index of the $k$th ranked item of $\mathcal{T}_{i}$ in our prediction. $M_{ij}$ is the true rating of item $j$ by user $i$ in the given dataset, and $\pi_{u}^{*}$ is the permutation that maximizes DCG@$K$. This measure counts only the top $K$ items in our predicted ranking and put more weights on the prediction of highly ranked items. We measured NDCG@$10$ in our experiments. Precision@$K$ is the ranking measure for binary ratings. Precision@$K$ for user $i$ is defined as

where $M_{ij}$ is the binary rating on item $j$ by user $i$ given in the dataset. This counts the number of relevant items in the predicted top $K$ recommendation. These two measures are averaged over all of the users.

We first compare our algorithm with numerical rating based algorithms, CofiRank and LCR. We follow the standard setting that are used in the collaborative ranking literature . For each user, we subsampled $N$ ratings, used them for training, and took the rest of the ratings for test. The users with less than $N+10$ ratings were dropped out. Table 2 compares AltSVM with numerical rating based algorithms. While $N=20$ is too small so that a global ranking provides the best NDCG, our algorithm performs the best with larger $N$. We also ran our algorithm with subsampled pairwise comparions with the largest numerical gap (AltSVM-sub), which are as many as $N$ for each user (the number of numerical ratings used in the other algorithms). Even with this, we could achieve better NDCG. We can also observe that the statistical performance is better with the hinge loss than with the logistic loss.

We have also experimented with collaborative ranking on binary ratings. We compare our algorithm against RobiRank , which is a recently proposed algorithm for collaborative ranking with binary ratings. We ran an experiment on a binarized version of the Movielens1m dataset. In this case, the movies rated by a user is assumed to be relevant to the user, and the other items are not. Since it is inefficient to take all possible comparisons which are in average a half million per user, we subsampled $C$ comparisons for each user. Both algorithms are set to estimate rank-100 matrices. Table 3 shows that our algorithm provides better performance than RobiRank.

3 Computational speed and Scalability

We now show the computational speed and scalability of our practical algorithm, AltSVM. The experiments were run on a single 16-core machine in the Stampede Cluster at University of Texas.

Figures 2a and 2b show NDCG@10 over time of our algorithms with 1, 4, and 16 threads, compared to CofiRank. Figure 2c shows Precision@10 over time of our algorithm with $C=5000$. We note that our algorithm converges faster, while the sample size $|\Omega|$ for our algorithm is larger than the number of training ratings that are used in the competing algorithms. Table 4 shows the scalability of AltSVM. We measured the time to achieve $10^{-5}$ tolerance on the binarized MovieLens1m dataset. As can be seen in the table, we could achieve significant speedup.

Conclusion

We considered the collaborative ranking problem where one fits a low-rank matrix to the pairwise comparisons by multiple users. We showed that the convex relaxation of the empirical risk minimization provides good generalization guarantees. For the large-scale practical settings, we also proposed a non-convex algorithm, which alternately solves two SVM problems. Our algorithm was shown to outperform the existing ones and parallelizes well.

References

Appendix A Proof of Theorem 3.1

We write $L(X)$ for the function being optimized; i.e.,

By some algebraic of manipulations $L$, we reduce the problem to showing a uniform law of large numbers for the family of functions $\{L(X):X\in\sqrt{\lambda d_{1}d_{2}}K\}$.

Using symmetrization and duality properties of $K$, we reduce the problem to bounding the norm of a matrix $M$ whose entries are sums of random signs.

We bound the norm of $M$ using various concentration inequalities and a theorem of Seginer .

In other words, it suffices to show a uniform law of large numbers for $\{L(X):X\in\sqrt{\lambda d_{1}d_{2}}K\}$.

Let $\epsilon_{i,j,k}$ be i.i.d. $\pm 1$-valued variables and let $\xi_{i,j,k}$ be the indicator that $(i,j,k)\in\Omega$. By Giné-Zinn’s symmetrization (as in ),

Since $\mathcal{L}$ is 1-Lipschitz, we obtain

where in the last line, we recognized that $\epsilon_{i,j,k}Y_{i,j,k}$ has the same distribution as $\epsilon_{i,j,k}$. Now, let $M$ denote the matrix where $M_{ij}=\sum_{k}(\xi_{i,j,k}\epsilon_{i,j,k}-\xi_{i,k,j}\epsilon_{i,k,j})$. Then

Together with the following lemma (which we prove in Appendix B), this completes the proof of Theorem 3.1

Appendix B Proof of Lemma A.1

We will decompose $M$ into two parts, $M=M^{(1)}-M^{(2)}$, with

Then $\|M\|\leq\|M^{(1)}\|+\|M^{(2)}\|$. Since $M^{(1)}$ and $M^{(2)}$ have the same distribution,

and so we are reduced to studying $M^{(1)}$, which has i.i.d. entries. Now, we apply Seginer’s theorem :

where $M^{(1)}_{i*}$ denotes the $i$th row of $M^{(1)}$ and $M^{(1)}_{*j}$ denotes the $j$th column, and $\|\cdot\|_{2}$ denotes the Euclidean norm.

Next, we will consider the size of the elements in $M^{(1)}$. First of all, $M^{(1)}_{ij}\leq Z_{ij}$ (this fairly crude bound will lose us a factor of $\sqrt{\log(d_{1}d_{2})}$). Now, Bernstein’s inequality applied to $Z_{ij}$ gives

Taking a union bound over $i$ and $j$, if $t\geq C\kappa\log(d_{1}d_{2})$ then

The same argument applies to $M^{(1)}_{*j}$ (but with $\sqrt{pd_{1}}$ instead of $\sqrt{pd_{2}}$), and so we conclude from (11) that

Appendix C Proof of Theorem 3.2

The proof of Theorem 3.2 uses Fano’s inequality.

On the other hand, $R_{1}(X^{1})$ and $R_{1}(X^{2})$ differ by $\Theta(\gamma^{2})$, because for a constant fraction of triples $i,j,k$, the chance that $Y_{i,j,k}$ is 1 differs by $O(\gamma)$ in $X^{1}$ and $X^{2}$, and on the event that $Y_{i,j,k}$ differs in these two models the loss differs by another $O(\gamma)$ factor.

C.2 Some concentration lemmas

We begin by quoting some standard concentration results (see, e.g. ).

One can easily show that the product of two subgaussian variables is subexponential:

If $X$ is $\sigma^{2}$-subgaussian and $Y$ is $\tau^{2}$-subgaussian then $XY$ is $C\sigma\tau$-subexponential for a universal constant $C$.

Moreover, one has a Bernstein-type inequality for sums of independent subexponential variables.

If $X_{1},\dots,X_{k}$ are i.i.d. $L$-subexponential then

C.3 Construction of a packing set

Let $0<\gamma<1$ be some parameter to be determined such that $B:=\lambda\gamma^{-2}$ is an integer.

Suppose that $\mathcal{L}^{\prime}(0)<0$. For every sufficiently small $\gamma$ (depending on $\mathcal{L}$), there exists a set $\mathcal{X}\subset\sqrt{\lambda d_{1}d_{2}}K$ of $\exp(cBd_{2})$ $d_{1}\times d_{2}$ matrices such that for any two $X^{1},X^{2}\in\mathcal{X}$,

Following Davenport et al., we construct this set $\mathcal{X}$ randomly: let $X$ be a random $B\times d_{2}$ matrix, where each element is chosen independently to be either $\gamma$ or $-\gamma$.

Let $X^{1}$ and $X^{2}$ be independent copies of $X$. Then with probability at least $1-\exp(-cBd_{2})$,

where $f(x)=e^{x}/(1+e^{x})$ is the logistic function, and the last line follows from a Taylor expansion of $D(f(x)\|f(y))$ around $x=y$, because all the $X^{1}_{ij}$ and $X^{2}_{ij}$ are bounded by $\gamma<1$. Together with (16), this proves the first inequality in Proposition C.4; the second inequality follows because each term of the form $D(f(X_{ij}-X_{ik})\|f(Y_{ij}-Y_{ik}))$ is bounded by a constant times $\gamma^{2}$. This proves the second inequality of Proposition C.4.

By Taylor expansion again, if $\gamma$ is sufficiently small (depending on $\mathcal{L}$) then

The same holds when $i,j,k$ is a triple for which $-2\gamma=X^{1}_{i,j}-X^{1}_{i,k}<X^{2}_{i,j}-X^{2}_{i,k}$. Finally, if $i,j,k$ is a triple such that $X^{1}_{i,j}-X^{1}_{i,k}=X^{2}_{i,j}-X^{2}_{i,k}$ then the expectation is zero. Summing over all triples, we see that on the event that Lemma C.5 holds,

After summing over all $\lceil d_{1}/B\rceil$ blocks, this proves the first inequality of Proposition C.4.

We may study each of the cross-terms separately: for the $X_{ij}Y_{ik}$ term, note that $\sum_{j}X_{ij}$ and $\sum_{k}Y_{ik}$ are both $\gamma^{2}d_{2}$-subgaussian (by Hoeffding’s inequality). Hence, $\sum_{jk}X_{ij}Y_{ik}$ is $C\gamma^{2}d_{2}$-subexponential (by Lemma C.2) and so by Lemma C.3,

The similar argument applies to the $X_{ij}X_{ik}$ term: $\sum_{j}X_{ij}$ is $\gamma^{2}d_{2}$-subgaussian and so $\sum_{ijk}X_{ij}X_{ik}=\sum_{i}(\sum_{j}X_{ij})^{2}$ is $C\gamma^{2}d_{2}$-subexponential; hence

Of course, the $Y_{ij}Y_{ik}$ term is identical. Finally, note that $\sum_{ijk}X_{ij}Y_{ij}=d_{2}\sum_{ij}X_{ij}Y_{ij}$. Since the terms in this sum are i.i.d., we may apply Hoeffding’s inequality to obtain

Putting everything together, we see that with high probability, the total of all the cross-terms in (17) is at most half of the first term. ∎

C.4 Completing the proof

Let $C$ denote the constant from Proposition C.4. Assume that $d_{1}\leq d_{2}$ and that $m$ is large enough so

Note that under the assumptions $\lambda\geq 1$ and $m\geq d_{1}+d_{2}$ from Theorem 3.2, the lower bound of (18) is satisfied. Moreover, if the upper bound of (18) is not satisfied then we may decrease $\lambda$ until it is; the conclusion of Theorem 3.2 will not be affected because as long as (18) fails, the minimum in Theorem 3.2 will be 1.

By the lower bound in (18), there is an integer $B$ such that

By the upper bound in (18), $\gamma\leq 1$.

Finally, note that by the first inequality in Proposition C.4, the error incurred by choosing the wrong $X\in\mathcal{X}$ is at least $c\gamma^{2}\asymp\sqrt{\frac{\lambda d_{2}}{m}}$.

Now, we have so far only discussed the case $d_{2}\geq d_{1}$. The case $d_{1}\leq d_{2}$ is not exactly equivalent because our model is not symmetric in its treatment of users and items. However, the proof of Theorem 3.2 does not change very much. We take horizontally stacked blocks of size $d_{1}\times B$ instead of $B\times d_{2}$. The main difference is in the calculation leading to (16): there are extra cross-terms appearing due to the fact that items in different blocks need to be compared with one another. However, all of these additional terms may be controlled with Lemmas C.2 and C.3 in much the same way as the existing terms are controlled.

Appendix D Comparison to Stochastic Gradient Descent

Another practical algorithm to optimize (3) is Stochastic Gradient Descent (SGD). We have experimented SGD on the same datasets in Table 1. We ran the algorithm with the same regularization parameters and different step sizes. The statistical results for SGD were observed to be no better than AltSVM, and hence we did not present them in the main paper.

Let us first describe the SGD procedure. At each step, ones chooses a triple $(i,j,k)\in\Omega$ uniformly at random and run a SGD step, which can be written as

where $\Omega^{(j)}$ denotes the number of comparisons in $\Omega$ which involve item $j$. $\eta$ is a step size and $g\in\partial\mathcal{L}(u_{i}^{\top}(v_{j}-v_{k}))$.

The following tables show the statistical result of SGD. The step size is chosen by $\eta=\frac{\alpha}{1+\beta t}$ as suggested in . $\alpha$ and $\beta$ were the powers of $10^{-1}$, and the best result is reported. The results are comparable to AltSVM, but it did not achieve better results. We note that this is the best result from several different step sizes, while AltSVM does not have any other parameter to choose except for the regularization parameter.