0.0.0.1 Funding information

This work is supported by grants of the Dutch Research Council (NWO), research program Vidi, grant numbers VI.Vidi.195.141 (Andreas Alfons) and VI.Vidi.211.032 (Ines Wilms).

1 Introduction↩︎

Recommender systems play a critical role in modern digital platforms, particularly in online retail and advertising, where they help users discover relevant products efficiently and drive substantial commercial value. Improvements in recommendation accuracy, even marginally, can yield significant financial gains in industries such as e-commerce and computational advertising [1]. However, malicious agents may launch so-called attacks with fake user profiles to deliberately manipulate recommendations to their advantage (e.g., [2]; see also [3], and [4], for overviews). For instance, already in the early 2000s, platforms like Amazon and eBay have reported repeated attempts to subvert their recommender systems through deceptive practices such as fake reviews and purchased ratings [5]. Today, these concerns remain more relevant than ever, as highlighted by [6], and they are drawing increasing attention from regulatory authorities such as the [7], the [8], or the [9]. For instance, the [10] recently introduced new measures and regulations aimed at protecting consumers, which include urging online platforms to proactively detect and address fake profiles and suspicious behavioral patterns—such as those generated by artificial intelligence. This paper introduces a novel matrix completion method that can withstand adversarial manipulation in recommender systems; our experiments offer a statistically-sound yet realistic blueprint for future studies to evaluate performance.

From a statistical perspective, recommender systems are often framed as matrix completion problems for a sparse user-item rating matrix (e.g., ranging from one star to five stars), with predictions being used to recommend new items to users. A prominent example is the famous Netflix Prize competition [11]. While early work on matrix completion dates back to, among others, [12], [13], and [14], it has been actively studied since [15] and [16]. For overviews, we refer to, e.g., [17], [18], and [1].

A multitude of methods for matrix completion exist, including weighted low-rank matrix approximation through singular value decomposition [19]; matrix factorization (e.g., [13], [14]; or [20], for recent computational advances); rank minimization, among which nuclear norm minimization [15], [16], [21], singular value thresholding [22]–[24], iteratively reweighted least squares minimization [25]; combinations of the former [26], [27]; or dynamic approaches over time [28].

We contribute to the literature, first methodologically, by presenting a novel matrix completion procedure called robust discrete matrix completion (RDMC) that is tailored towards the discrete nature of incomplete rating-scale data while being robust to adversarial manipulation. We assess RDMC against the widely-used procedure Soft-Impute [21], [26], as well as a variant that we adjusted with a post-hoc discretization step to fit the discrete nature of rating-scale data. We study these methods in empirical case studies and extensive simulations. These experiments are carefully designed to be more realistic than previous studies, thereby providing an important second conceptual contribution to the literature, as our experiments may serve as realistic yet statistically-sound blueprints for future studies in the context of matrix completion and recommender systems to further boost transparency, reproducibility and generalizability. With our contributions, we in particular address the following issues.

First, the vast majority of existing methods for matrix completion are formulated over the real number domain and produce continuous predictions, despite recommender systems applications featuring (discrete and bounded) rating-scale data. Exceptions that include discreteness or box constraints can be found in [29], [30], [31], [32], [33], or [34]. Although the observed ratings may be interpreted as discrete measurements of a latent continuous sentiment, the underlying continuous distributions are not identified without additional assumptions [35], invalidating comparisons of predictions across different items. The proposed procedure RDMC therefore restricts the predictions to the given rating scale via a discreteness constraint.

Second, only a relatively small subset of the literature on matrix completion considers the presence of outliers or corrupted/manipulated observations, such as fake profiles in recommender systems. From a perspective of robust statistics, there are primarily two approaches for handling outliers to avoid biased analyses. One approach involves explicitly identifying and removing outliers, while the generally favored approach focuses on robust methods that can accommodate the presence of outliers without compromising performance [36]. In the context of matrix completion and recommender systems, an example of the first strategy is given in [37], while existing procedures for the latter are commonly based on robust nuclear norm minimization [29], [38]–[42] or robust matrix factorization [43]–[47]. Among these, only [29] address discrete data. Moreover, [6] cautions against the prevalent business practice of merely identifying and removing fraudulent reviews, as their study demonstrates that such an approach can result in persistent, long-term negative effects. This highlights the need for a robust statistical approach to protect against adversarial manipulation in matrix completion for discrete rating-scale data, which RDMC achieves by incorporating a robust loss function.

Third, much of the literature does not investigate the effect of missing data mechanisms other than missing completely at random (MCAR) on matrix completion procedures, with some notable exceptions being [48] for missing at random (MAR) as well as [49] and [50] for missing not at random (MNAR). Yet missing values in recommender systems applications are not MCAR in practice [1]: users predominantly rate items that they have consumed or purchased and thus expected to like, intrinsically linking the probability of missingness to their preferences. Hence, we study the performance of the proposed method under realistic MNAR settings.

For transparency and reproducibility—key concerns in the current scientific landscape as discussed in [1]—we provide implementations of the methodology in package RMCLab [51] for the statistical computing environment R [52], which is available from https://CRAN.R-project.org/package=RMCLab. The main part of the code is thereby written in C++ to improve computational efficiency. Furthermore, replication files of all analyses will be made publicly available upon acceptance of this manuscript.

The remainder of the paper is structured as follows. Section 2 introduces the proposed procedure RDMC for robust matrix completion with discrete rating-scale data. In Section 3, we introduce a simulation design that mimics recommender systems applications, and we investigate the performance of RDMC. Section 4 then presents two empirical case studies. The final Section 5 provides a concluding discussion.

2 A robust procedure for discrete matrix completion↩︎

We start by formulating the regularized optimization problem for matrix completion in Section 2.1 and presenting the corresponding algorithm in Section 2.2. Section 2.3 provides further algorithmic details, while Section 2.4 introduces different choices of robust loss functions. Finally, the selection of the regularization parameter is discussed in Section 2.5.

2.1 Problem formulation↩︎

Suppose that we observe an incomplete rating matrix \(\boldsymbol{R}\) of \(n\) rows (representing individuals providing the ratings) and \(p\) columns (representing the rated items) with elements \(R_{ij} \in \{1, 2, \dots, K\}\). For practical reasons (see Section 2.3 for further discussion), we consider a column-centered matrix \(\boldsymbol{X}\) with elements \(X_{ij} \in \mathcal{C}_{j} = \{c_{1}^{j}, \dots, c_{K}^{j}\} \subset \mathbb{R}\), where \(c_{1}^{j} < \dots < c_{K}^{j}\) for \(j = 1, \dots, p\). Note that the values of the rating categories may vary between columns of \(\boldsymbol{X}\). Hence, the assumption that the original rating matrix \(\boldsymbol{R}\) takes values in \(\{1, 2, \dots, K\}\) is purely for simplicity and without loss of generality, and the proposed procedure extends in a straightforward manner to settings where even the number of rating categories may differ per column.

Let \(\Omega\) denote the index set of observed entries in \(\boldsymbol{X}\) and define the projection \(P_{\Omega}\) to be the \((n \times p)\)-dimensional matrix whose elements are given by \[\left( P_\Omega(\boldsymbol{X}) \right)_{ij} = \begin{cases} X_{ij} & \text{ if } (i,j) \in \Omega, \\ 0 & \text{ otherwise.} \end{cases}\] Matrix completion is often based on minimizing the mean squared error for the observed elements subject to a nuclear norm constraint as a relaxation of a low-rank constraint. Specifically, the widely-used procedure Soft-Impute [21], [26] solves the optimization problem (in Lagrangian form) \[\label{eq:SI} \min_{\boldsymbol{L}} \frac{1}{2} \| P_{\Omega}(\boldsymbol{X}) - P_{\Omega}(\boldsymbol{L}) \|_{F}^{2} + \lambda \| \boldsymbol{L} \|_{*},\tag{1}\] where \(\| \cdot \|_{F}\) denotes the Frobenius norm and \(\| \cdot \|_{*}\) the nuclear norm, while \(\lambda \geq 0\) is a regularization parameter.

In order to preserve the discrete nature of the ratings, one could add a discreteness constraint to 1. However, requiring the elements of \(\boldsymbol{L}\) to be discrete may not yield a solution, since such a discreteness constraint and a rank constraint are unlikely to be fulfilled simultaneously. To resolve this issue, an ancillary continuous matrix \(\boldsymbol{Z}\) can be introduced into problem 1 such that the discreteness constraint operates on \(\boldsymbol{L}\) and the nuclear norm constraint on \(\boldsymbol{Z}\), while ensuring that \(\boldsymbol{L}\) and \(\boldsymbol{Z}\) remain as close as possible. In addition, the squared Frobenius norm in 1 does not protect against corrupted observations (such as fake profiles in recommender systems). Indeed, this norm penalizes large errors quadratically, implying that corrupted observations with large errors may dominate the loss function and significantly influence the recommender system. For increased protection, we replace it with a robust loss function that puts less emphasis on large errors, thereby preventing them from disproportionately affecting the recommender system. More specifically, we replace it with a (pseudo-)norm \(\| \boldsymbol{Y} \|_{\rho} = \sum_{i,j} \rho(Y_{ij})\) based on a robust loss function \(\rho\) [46]. Putting all of this together, we obtain the following optimization problem in augmented Lagrangian form: \[\label{eq:RDMC} \begin{array}{r l} {\displaystyle \min_{\boldsymbol{L}, \boldsymbol{Z}}} & {\displaystyle \| P_{\Omega}(\boldsymbol{X}) - P_{\Omega}(\boldsymbol{L}) \|_{\rho} + \lambda || \boldsymbol{Z}||_{*} + \langle \boldsymbol{\Theta}, \boldsymbol{L} - \boldsymbol{Z} \rangle_{F} + \frac{\mu}{2} || \boldsymbol{L} - \boldsymbol{Z}||_{F}^{2}} \\ \text{subject to} & L_{ij} \in \mathcal{C}_{j}, \quad i = 1, \dots, n, j = 1, \dots, p, \end{array}\tag{2}\] where \(\langle \cdot,\cdot \rangle_{F}\) denotes the Frobenius inner product, \(\boldsymbol{\Theta}\) is a multiplier adjusting for the discrepancy between \(\boldsymbol{L}\) and \(\boldsymbol{Z}\), and \(\mu\) is an additional regularization parameter. It should be noted that 2 generalizes the optimization problem formulated in [29] to a wider class of robust loss functions and by allowing different values of the rating categories between columns of the (column-centered) rating matrix.

2.2 Algorithm↩︎

The formulation of the optimization problem 2 lends itself to an alternating direction method of multipliers (ADMM) algorithm [53], which follows along similar lines as that of [29]. For a given value of the regularization parameter \(\lambda\), the following steps are iterated until convergence (after initialization as described in Section 2.3).

First, consider \(\boldsymbol{L}\) fixed and solve 2 for \(\boldsymbol{Z}\). Combining the terms for the Frobenius inner product and Frobenius norm, dropping constant terms, and division by \(\mu\) yields the equivalent minimization problem \[\min_{\boldsymbol{Z}} {\textstyle \frac{1}{2} \| (\boldsymbol{L} + \frac{1}{\mu} \boldsymbol{\Theta}) - \boldsymbol{Z} \|_{F}^{2} + \frac{\lambda}{\mu} \| \boldsymbol{Z} \|_{*} }.\] The solution is given by the soft-thresholded singular value decomposition (SVD) of \(\boldsymbol{L} + \frac{1}{\mu} \boldsymbol{\Theta}\) [22]. That is, with \(\boldsymbol{L} + \frac{1}{\mu} \boldsymbol{\Theta} = \boldsymbol{U} \boldsymbol{D} \boldsymbol{V}^{\top}\), where \(\boldsymbol{U}\) is \((n \times q)\)-dimensional, \(\boldsymbol{V}\) is \((p \times q)\)-dimensional, \(\boldsymbol{D}\) is a \((q \times q)\)-dimensional diagonal matrix whose diagonal elements are denoted by \(d_{1}, \dots, d_{q}\), and \(q \leq \min(n, p)\) denoting the rank, we obtain \[\boldsymbol{Z} = \boldsymbol{U} S(\boldsymbol{D}) \boldsymbol{V}^{\top}, \label{eq:Zupdate}\tag{3}\] where \(S(\boldsymbol{D}) = \text{diag} \left( (d_{1} - \frac{\lambda}{\mu})_{+}, \dots, (d_{q} - \frac{\lambda}{\mu})_{+} \right)\) with \((y)_{+} = \max(y, 0)\).

Second, consider \(\boldsymbol{Z}\) fixed and solve 2 for \(\boldsymbol{L}\). Using similar operations as described above but without rescaling, we arrive at the equivalent minimization problem \[\begin{array}{r l} {\displaystyle \min_{\boldsymbol{L}}} & \| P_{\Omega}(\boldsymbol{X}) - P_{\Omega}(\boldsymbol{L}) \|_{\rho} + \frac{\mu}{2} \| \boldsymbol{L} - \boldsymbol{Z} + \frac{1}{\mu} \boldsymbol{\Theta} \|_{F}^{2} \\ \text{subject to} & L_{ij} \in \mathcal{C}_{j}, \quad i = 1, \dots, n, j = 1, \dots, p. \end{array}\] This can be solved element-wise, with the solution being given by \[L_{ij} = \begin{cases} {\displaystyle \mathop{\mathrm{arg\,min}}_{c_{k} \in \mathcal{C}_{j}}} \; \rho(c_{k}- X_{ij}) + \frac{\mu}{2} (c_{k} - Z_{ij} + \frac{1}{\mu} \Theta_{ij})^2 & \text{ for } (i,j) \in \Omega, \\ {\displaystyle \mathop{\mathrm{arg\,min}}_{c_{k} \in \mathcal{C}_{j}}} \; (c_{k} - Z_{ij} + \frac{1}{\mu} \Theta_{ij})^2 & \text{ for } (i,j) \notin \Omega. \end{cases} \label{eq:updateL}\tag{4}\]

Third, update the discrepancy parameter \(\boldsymbol{\Theta} \leftarrow \boldsymbol{\Theta} + \mu (\boldsymbol{L} - \boldsymbol{Z})\) and the regularization parameter \(\mu \leftarrow \delta \mu\) with \(\delta > 1\). Hence, \(\mu\) grows exponentially in the number of iterations to speed up convergence.

Pseudo-code for the algorithm is presented in Algorithm 1. Following [29] and [46], we initialize \(\mu = 0.1\) and set \(\delta = 1.05\). As convergence criterion, we take the relative change in the objective function from 2 falling below a given threshold, which we set to \(\varepsilon_\text{tol} = 10^{-4}\). Furthermore, we set the maximum number of iterations to \(t_\text{max}=100\), as this sufficed for convergence in all our numerical experiments.

2.3 Column centering and initialization↩︎

Since the algorithm contains a soft-thresholded SVD step, it is important that the observed (incomplete) rating matrix \(\boldsymbol{R}\) is centered. One possibility is to center each column by the midpoint of the rating scale (i.e., the mean of the minimum and maximum rating category). While this can be expected to work well if the rating distributions of the columns are symmetric around the midpoint, such a setting is unrealistic in practical applications. In recommender systems, the distributions of popular items, for instance, are typically skewed towards higher ratings, with the maximum rating often being the most frequent. Hence, we center the \(j\)th column in the given data matrix \(\boldsymbol{R}\) by the median \(M_j\) of its observed cells, i.e., we obtain \(\boldsymbol{X}\) by setting \(X_{ij} = R_{ij} - M_{j}\), \(i = 1, \dots, n\), \(j = 1, \dots, p\). Accordingly, we transform the set of rating categories to \(C_{j} = \{c_{1}^{j}, \dots, c_{K}^{j}\}\) with \(c_{k}^{j} = k - M_{j}\) for \(k = 1, \dots, K\). This highlights the need for a procedure that allows for different rating categories in different columns, a feature that our proposed procedure RDMC accommodates.

Using the median-centered matrix \(\boldsymbol{X}\), we initialize the matrix \(\boldsymbol{L} = P_{\Omega}(\boldsymbol{X})\) (corresponding to median imputation) and the discrepancy parameter \(\boldsymbol{\Theta}\) as an \((n \times p)\)-dimensional matrix of zeros. However, the algorithm is typically applied for a grid of values for \(\lambda\) (see Section 2.5). Then the obtained solutions for \(\boldsymbol{L}\) and \(\boldsymbol{\Theta}\) for a given value of \(\lambda\) are used as starting values for the next value of \(\lambda\). As \(\boldsymbol{\Theta}\) adjusts for the discrepancy between \(\boldsymbol{L}\) and \(\boldsymbol{Z}\) (which should increase with increasing \(\lambda\)), the values of \(\lambda\) should thereby be sorted in ascending order.

2.4 Loss functions↩︎

In order to reduce the influence of large errors due to corrupted observations such as fake profiles, we consider the following robust loss functions for the (pseudo-)norm \(\|\cdot\|_{\rho}\) in 2 , which are common choices in the literature on robust methods:

• The pseudo-Huber loss \(\rho(y) = \tau^2 (\sqrt{1 + (y / \tau)^2} - 1)\) with parameter \(\tau\). The inclusion of this loss function is motivated by its successful application to rating-scale data in a different context, namely in autoencoder neural networks for the detection of careless responding in surveys from the behavioral sciences [54], [55]. Here, we suggest to link the parameter \(\tau\) to the step size between rating categories by setting \(\tau = 1\).

• The absolute loss \(\rho(y) = |y|\).

• A truncated variant of the absolute loss given by \(\rho(y) = \min (|y|, \tau)\). We suggest to set \(\tau = (K - 1) / 2\), i.e., the loss is truncated at half the range of the rating categories.

Since the update step for \(\boldsymbol{L}\) yields separable problems for its elements (see Section 2.2), the use of a nonconvex function such as the truncated absolute loss comes at no cost to computational complexity. A nonconvex loss may yield greater robustness in selecting the regularization parameter \(\lambda\) when minimizing the loss on a test set, which is discussed next.

2.5 Selection of the regularization parameter↩︎

To select the regularization parameter \(\lambda\) from a grid of candidate values based on out-of-sample prediction performance, some of the observed elements of \(\boldsymbol{R}\) can be set to missing values for subsequent use as a validation set. For measuring the prediction error on the validation set, a natural choice is to apply the same loss function \(\rho\) that is used for fitting the algorithm on the training set. It is possible to apply a cross-validation scheme whereby the observed elements are randomly divided into blocks, with each block being used as validation set once. However, even if the rows are independent, elements within the same row are not. It is therefore unclear if cross-validation would reduce the correlations among prediction errors on the different validation sets, compared to repeated holdout validation whereby a proportion of the observed elements are randomly selected in each replication to form the validation set. Hence, we prefer repeated holdout validation, as the proportion of observations in the validation set and the number of replications can be chosen independently.

3 Simulations↩︎

Before applying the proposed procedure in empirical case studies, we thoroughly evaluate it via simulations. Crucially, our simulation design closely mimics recommender systems applications so that it may serve as blueprint for future studies. We perform 100 replications.

3.1 Data generation↩︎

We simulate data of \(n = 300\) user ratings on \(p = 200\) items. We start by simulating latent continuous data from the low-rank matrix factorization model \(\boldsymbol{Z} = \boldsymbol{A} \boldsymbol{B}^{\top} + \mathcal{E}\), where \(\boldsymbol{A}\) is of dimension \(n \times q\), \(\boldsymbol{B}\) of dimension \(p \times q\), and \(\boldsymbol{\mathcal{E}}\) of dimension \(n \times p\), with rank \(q = 20\) and the elements of \(\boldsymbol{A}\), \(\boldsymbol{B}\), and \(\boldsymbol{\mathcal{E}}\) being independent and standard normally distributed. Subsequently, we rescale \(\boldsymbol{Z}^{*} = \boldsymbol{Z} / \sqrt{q+1}\) so that the elements of \(\boldsymbol{Z}^{*}\) have variance 1.

To create popular items (in general higher ratings) and unpopular items (in general lower ratings), we add random mean shifts \(s_{1}, \dots, s_{p}\) to the respective columns of \(\boldsymbol{Z}^{*}\) before discretization into \(K \in \{3, 5, 10\}\) ordinal rating categories encoded as values \(\{1, \dots, K\}\). These mean shifts are randomly drawn from the interval \([-s_{\text{max}}, s_{\text{max}}]\), where \(s_{\text{max}}\) depends on the number of categories and breakpoints in the discretization. Specifically, \(s_{\text{max}}\) is chosen so that the corresponding mean shift results in 40% of the ratings being expected in the maximum rating category.

For discretizing the resulting continuous data matrix in order to obtain a rating matrix \(\boldsymbol{R}\), we use the following breakpoints depending on the number of rating categories:

• \(K = 3\): inspired by Netflix’ asymmetric rating scale (“Not for me”, “I like this”, “Love this”), we set the breakpoints between categories at 0 and 1.5.

• \(K = 5\): to simulate common 1 to 5 star ratings, we set the breakpoints at \(-1.5\), \(-0.5\), \(0.5\), and \(1.5\).

• \(K = 10\): to simulate common 1 to 10 star ratings (or, equivalently, ratings up to 5 stars but allowing for half-stars), we set the breakpoints at \(-2\), \(-1.5\), \(-1\), \(-0.5\), \(0\), \(0.5\), \(1\), \(1.5\), and \(2\).

3.2 Missing values↩︎

We generate missing values in the rating matrix \(\boldsymbol{R}\) using two different settings:

• Missing not at random (MNAR): the negated mean shifts \(-s_{1}, \dots, -s_{p}\) are mapped to the interval \([0.4, 0.99]\) in order to obtain the proportion of observations in each item to be replaced by missing values. For instance, the most popular item (large positive mean shift, high ratings) contains 40% missing values, while the most unpopular item (large negative mean shift, low ratings) contains 99% missing values.² Hence, low ratings are much more likely to be missing than high ratings.

• Missing completely at random (MCAR): a proportion \(\gamma = 0.7\) of cells in \(\boldsymbol{R}\) is randomly selected and replaced by missing values. Here, \(\gamma\) is chosen so that the overall proportion of missing values is similar to the MNAR setting. Note that this setting is included for reference purposes since MCAR is unrealistic in recommender systems applications.

3.3 Attacks↩︎

To investigate the robustness of the methods against adversarial manipulation with profile injection attacks (also known as shilling attacks), we focus on so-called nuke attacks aimed at demoting a certain target item, i.e., decreasing the probability of it being recommended. We apply three efficient attack schemes, denoted by average, reverse bandwagon, and love/hate (see [57] for a detailed overview of common attack schemes).

To select the target item, we first make a pre-selection of items with a mean shift larger than \(0.9 \cdot s_{\text{max}}\). Although it requires knowledge of unobserved information, this pre-selection reflects that the attacker has some prior notion of popular items. Among those, the item with the highest average observed rating is targeted, with fake profiles assigning the minimum rating to this item. To keep the effect of the attacks comparable across missing data settings, the number of fake profiles is determined as a proportion \(\varepsilon = 0.2\) of the number of observed ratings in the target item.

Table 1 summarizes the strategy in the three attack schemes regarding so-called selected items (which are chosen by the attacker based on specific characteristics) and filler items (which are randomly chosen). In the average attack, filler items are assigned ratings based on the item mode. The idea is that the fake profiles have typical ratings in other items but dislike the target item. We thereby use the item mode instead of the item average, following the recommendation of [58] for discrete ratings. The reverse bandwagon attack is intended to associate the target item with disliked items. Hence, it uses unpopular items as selected items and assigns the minimum rating to both the selected and filler items. As unpopular items, we take the items with the lowest average observed ratings, provided that at least 20 ratings are observed. In the love/hate attack, the filler items receive the maximum rating so that the fake profiles love any other item while hating the target item.

3.4 Methods↩︎

For the proposed method RDMC, we consider the three loss functions described in Section 2.4 (denoted by pHuber for the pseudo-Huber loss, absolute for the absolute loss, and truncated for the truncated absolute loss) and set other parameters as described in Section 2.2.

We compare RDMC with the popular procedure Soft-Impute (SI) [21], [26]. Since this method yields continuous predictions, we introduce a variant that discretizes the obtained predictions to the given rating scale via the mapping function \(m_{j}(y) = \min(\max([y], c_{\min}), c_{\max})\), where \(c_{\min}\) and \(c_{\max}\) are the minimum and maximum rating, respectively, and \([\cdot]\) denotes rounding to the nearest integer (SI-discretized). We use the SVD-based algorithm due to its theoretical convergence guarantees and set the same convergence threshold as for RDMC.

For selecting the regularization parameter \(\lambda\), both RDMC and SI use repeated holdout validation with 5 replications and 10% of the observed cells to be randomly selected as validation set. The candidate values are thereby given by a logarithmic grid of ten values between 0.01 and 1, which are scaled by the largest singular value of \(P_{\Omega}(\boldsymbol{X}_{\text{train}})\), with \(\boldsymbol{X}_{\text{train}}\) denoting the median-centered (for RDMC) or mean-centered (for SI) training data.

Additional benchmark methods are median imputation (denoted by median), a discretized variant of median imputation (median-discretized; in case the median of a column falls in between rating categories, the predictions are randomly sampled from the corresponding two categories), and mode imputation (mode; in case of a column with multiple modes, one of them is selected at random for each missing cell).

3.5 Evaluation criteria↩︎

In the absence of an attack, the methods are evaluated using the mean absolute error (MAE) over the predictions of all missing cells, i.e., \[\label{eq:MAE} MAE = \frac{1}{|\Gamma|} \sum_{(i,j) \in \Gamma} |R_{ij} - \widehat{R}_{ij}|,\tag{5}\] where \(\Gamma\) denotes the index set of missing cells in \(\boldsymbol{R}\), and \(\widehat{R}_{ij}\) is the prediction of the cell \(R_{ij}\). For the different attack settings, we evaluate the methods via the mean prediction shift (MPS) in the target variable [57], i.e., \[MPS = \frac{1}{|\Gamma_{\text{target}}|} \sum_{i \in \Gamma_{\text{target}}} \left (\widehat{R}_{ij_{\text{target}}}^{\text{after}} - \widehat{R}_{ij_{\text{target}}}^{\text{before}} \right),\] where \(\Gamma_{\text{target}}\) is the index set of missing cells in the target variable \(j_{\text{target}}\) of \(\boldsymbol{R}\), while \(\widehat{R}_{ij}^{\text{before}}\) and \(\widehat{R}_{ij}^{\text{after}}\) denote the prediction of the cell \(R_{ij}\) before and after the attack, respectively.

3.6 Results↩︎

For RDMC, we focus on the pseudo-Huber loss, as preliminary simulations indicate that results for the absolute loss and the truncated absolute loss are similar (see Appendix 6). Figure 2 displays box plots of the MAE for the setting without an attack. In case of three rating categories, SI without the discretization step does not improve upon median and mode imputation. The discretized variant, on the other hand, yields a clear improvement in MAE, with RDMC falling in between. For five rating categories, SI and RDMC considerably outperform the benchmark methods, with RDMC being in between the two SI approaches for the MNAR setting, but performing equally well as SI-discretized in the MCAR setting. For ten rating categories, RDMC yields the best performance in terms of MAE, while both SI approaches exhibit large variability in the MCAR setting.

We now turn to stability of the predictions under adversarial manipulation, with Figure 3 containing box plots of the MPS for the MNAR setting. Results for the MCAR setting are qualitatively similar and can be found in Appendix 6. Both SI approaches are strongly influenced by the attacks, with the average attack having the biggest impact. The MPS is thereby large enough that the target item would be recommended to far fewer users—for instance, the predictions on average drop by at least one full rating category in the setting with ten categories. RDMC is far more stable in the presence of attacks, with an MPS much closer to zero in all settings. Furthermore, while median regression is stable for three and five rating categories, it yields an MPS of \(-1\) in most instances with ten rating categories.

On a final note, we also included a variant of RDMC with the pseudo-Huber loss that we do not iterate until convergence. Instead, we apply a liberal stopping criterion, limiting the procedure to a maximum of 10 iterations (denoted by RDMC-pHuber-liberal in the figures). Interestingly, the MAE in case of no attack is often smaller compared to iterating until convergence, whereas the MPS in the presence of an attack is slightly larger in absolute value in case of three rating categories.

4 Empirical case studies↩︎

We consider two publicly available datasets in the context of recommender systems. The first in Section 4.1 is used to investigate the impact of attacks on empirical data (rather than stylized simulated data), while the second in Section 4.2 allows to compare two missing data mechanisms (MNAR and MCAR). In both case studies, the observed ratings are split randomly into subsets of 80% for training the algorithms and 20% for testing prediction performance.

Since computational burden forms an important practical consideration across many recommender system applications, we also analyze the stability of RDMC and SI with respect to their stopping criteria in order to assess whether reducing the computational burden comes at a cost of accuracy. We use two stopping criteria for each method, denoted by liberal and strict. For RDMC, the liberal criterion corresponds to stopping after a maximum of 10 iterations, whereas the strict criterion iterates until convergence (cf. our simulations from Section 3). For SI, on the other hand, we follow the strategy of [26] by setting a larger convergence threshold of \(10^{-3}\) for the liberal criterion, while the strict criterion uses the threshold \(10^{-4}\) (as in the simulations). The other parameters are kept the same as in Section 3, except that the number of replications in repeated holdout sampling for selecting the regularization parameter is increased from 5 to 10.

4.1 MovieLens 100K data: The impact of attacks↩︎

The famous MovieLens 100K dataset [56] is extensively analyzed in previous studies on recommender systems [46], [48], [59]. It consists of 100,000 movie ratings on a scale from 1 to 5, of \(n = 943\) users on 1,682 movies. We restrict our analysis to the \(p = 939\) movies that are rated at least 20 times. The observed ratings follow a left-skewed distribution with a majority of medium to high ratings, as shown in Figure 4.

First, Figure 5 (left) displays the prediction performance of the different methods in terms of the MAE on the test set. RDMC (with any of the three loss functions) and SI clearly outperform median and mode imputation with an MAE of around 0.725, while SI-discretized presents an even smaller MAE around 0.69. Applying the liberal stopping criterion thereby hardly affects prediction performance of RDMC and SI compared to the strict criterion. However, Figure 5 (right) shows that it yields a substantial reduction in computation time for RDMC—by more than a factor of three—resulting in similar computation time to SI.

Second, we investigate the sensitivity of the methods to adversarial manipulation. To this end, we use the three nuke attack schemes described in Section 3.3. We additionally vary the number of fake profiles, based on a proportion \(\varepsilon \in \{0.10,0.15,0.20\}\) of the number of observed ratings in the target variable. Figure 6 displays the results for the MPS. All RDMC approaches outperform both SI approaches by a considerable margin, for all attack schemes. Moreover, RDMC provides effective protection against the attacks with an MPS close to 0 in almost all cases, irrespective of the loss function. In the presence of the love/hate attack, however, the strict stopping criterion is necessary for RDMC to achieve the best result. As expected, the MPS of the two SI approaches deteriorates as the size of the attack increases. Interestingly, SI-discretized falls behind SI when the attack size exceeds 15%.

4.2 Yahoo! Music data: The impact of missingness mechanisms↩︎

The Yahoo! Music ratings data contain two datasets of users rating songs on a scale of 1 to 5. In the first dataset (User Selected), the users choose the songs they want to rate, whereas in the second dataset (Randomly Selected), 10 songs from a pool of \(p_{2} = 1,000\) are randomly assigned to them. Hence, the former constitutes an MNAR mechanism while the latter represents the MCAR mechanism. To ensure a fair comparison across datasets, we restrict our analysis to the \(n = 2,361\) users who are present in both datasets and rated at least 20 songs in the first dataset. Furthermore, we limit the first dataset to the \(p_{1} = 998\) songs with at least 20 ratings. In both datasets, we have more than 95% of missing values (96% for MNAR and 99% for MCAR). Figure 7 displays the distributions of the observed ratings in the two datasets. The missing data mechanism clearly leads to different rating distributions. When the users can choose the songs they want to evaluate (MNAR, left plot), the distribution of the ratings is almost bimodal with 1-star ratings being the most frequent, followed by 5-star ratings. When the songs to be rated are randomly assigned (MCAR, right plot), the distribution is highly right-skewed with decreasing frequency for increasing ratings, and very few high ratings.

Figure 8 shows the MAE on the test set for both missing data mechanisms. As expected, the MAE of all methods is lower in case of MCAR compared to MNAR. Overall, RDMC outperforms the SI approaches in both scenarios. For RDMC, the loss function and the number of iterations seem to have an effect only for the MNAR scenario, with a small advantage for the absolute loss function and perhaps a minor improvement with the liberal stopping criterion. For SI, the discretized variant indicates only a minimal improvement over the standard algorithm, but the stopping criterion has a major impact for the MNAR scenario: the strict stopping criterion is required to ensure performance similar to RDMC. Nevertheless, as illustrated in Figure 9, RDMC with the liberal stopping criterion remains competitive with SI in terms of computation time.

5 Conclusions and discussion↩︎

In a recent review paper, [1] conclude that too few statisticians are engaged in research on recommender systems despite their potential to make consequential contributions. But to fulfill this potential, it is vital that recommender systems are studied under more realistic scenarios. Our simulations and case studies may thereby provide a blueprint—or at least useful starting values—on how to design such studies.

Moreover, we present a novel method for robust matrix completion with discrete rating-scale data (RDMC), as well as an extension of the popular method Soft-Impute (SI) with a post-hoc discretization step (SI-discretized). RDMC thereby combines a robust loss function on the errors for the observed ratings with a discreteness constraint on the predictions and a low-rank constraint on an ancillary continuous matrix. Extensive simulations and two empirical case studies focused on attacks with fake user profiles and missing not at random (MNAR) mechanisms demonstrate that RDMC protects well against adversarial manipulation, while paying only a small price in terms of performance in the absence of corrupted observations compared to SI-discretized. The choice among the three considered loss functions does not yield considerable differences in performance, with the pseudo-Huber loss being recommended overall due to in general requiring fewer iterations for convergence. Moreover, performing only a small number of iterations (rather than iterating until convergence) can reduce computational effort, in many cases without a detrimental effect on performance. Regarding Soft-Impute, the discretization step results in an improvement over its standard variant in the absence of corrupted observations, in particular in settings with few rating categories (five or less), but it remains equally vulnerable to fake profiles.

Our method shows promise for broader applicability across diverse domains. For instance, empirical research in social and behavioral sciences relies heavily on rating-scale surveys, but listwise deletion remains the predominant method of handling missing data in some fields [60]. Preliminary simulations suggest that RDMC represents a notable advancement in robust imputation of rating-scale data (see Appendix 7), even in the presence of inattentive respondents or bot respondents—a critical issue in online surveys. Furthermore, our approach holds potential for contexts involving so-called strategic missing data situations [61]—such as financial reporting, college admissions, job applications, and marketing—where individuals may intentionally withhold information to achieve favorable outcomes. Finally, recent research has shown a notable rise in matrix completion methods for panel data [49], [62], [63], providing a promising direction for future adaptation of our robust methodology.

Data availability↩︎

The MovieLens 100K data are publicly available from https://grouplens.org/datasets/movielens/100k/. The Yahoo! Music data are available through the Yahoo! Webscope data sharing program at https://webscope.sandbox.yahoo.com/catalog.php?datatype=r. Various datasets are available through this link; we use the dataset titled “R3 - Yahoo! Music ratings for User Selected and Randomly Selected songs, version 1.0”, downloaded on September 9, 2024.

6 Additional simulation results↩︎

Figure 10 contains preliminary results for RDMC with different loss functions from 50 replications of the simulated recommender system with five rating categories, where the left plot displays the MAE for the setting without an attack and the right plot presents the MPS for different attack schemes. All three loss functions yield similar results, while the pseudo-Huber loss typically requires fewer iterations to converge (see Figure 11).

Figure 12 shows the MPS under different attacks for the simulated recommender system in the MCAR setting. The results are qualitatively similar to those of the MNAR setting, but a notable difference is that there is a small number of instances where RDMC yields a relatively large negative prediction shift. A possible explanation why this occurs in the MCAR setting but not in the MNAR setting is that in the former, all variables including the target variable contain equally sparse information. In practice, however, this is a rather unrealistic setting.

7 Rating-scale surveys↩︎

7.1 Introduction↩︎

Another relevant application of rating-scale data, although rarely considered in the literature on matrix completion, are surveys in the social and behavioral sciences. In such surveys, researchers measure each latent construct of interest (e.g., personality traits) by multiple rating-scale items (e.g., respondents may be asked how accurately certain adjectives describe their personality, with response categories ranging from 1 = “very accurate” to 5 = “very accurate”; [64]). Missing values are commonly occurring due to nonresponse in parts of the survey, with nonresponse often being missing not at random (MNAR). For instance, participants may abandon the survey yielding higher probability of missingness for later survey items, or they may not respond to items they find, e.g., confusing or inappropriate. In addition, careless responding (i.e., responses not based on item content due to inattention or misunderstanding; e.g., [65]) or bot responding [66] are common occurrences in online surveys, and a prevalence as low as 5% can invalidate research findings [64], [67]–[69].

In some fields, removing participants with missing responses is still the most common method of handling missing responses in survey data [60]. Yet this constitutes a loss of valuable information (the observed responses of the removed participants) and generally results in biased analyses [70]. Furthermore, the common presence of careless respondents or bot respondents requires imputation methods designed for discrete rating-scale data. We therefore investigate the use of the proposed matrix completion method RDMC in preliminary simulations motivated by surveys in the social and behavioral sciences.

7.2 Simulations↩︎

7.2.1 Data generation↩︎

We simulate rating-scale responses of \(n = 300\) participants to a survey on \(q = 10\) latent constructs, with each construct being measured by \(r \in \{4, 8\}\) items such that the total number of items is \(p \in \{40, 80\}\). For this purpose, we first simulate the underlying latent sentiments \(\boldsymbol{Z}\) from a multivariate normal distribution \(N(\boldsymbol{0}, \boldsymbol{\Sigma})\), where \(\boldsymbol{\Sigma}\) follows a block-Toeplitz structure with 1 on the diagonal and the remaining elements being given by \(\rho_{kl} = 0.6^{|k-l|+1}\) with \(k = 1, \dots, q\) and \(l = 1, \dots, q\) denoting the row and column indices of the blocks in \(\boldsymbol{\Sigma}\).

In survey-based research, skewness towards one side of the rating scale is common, for instance, in organizational and management research [71], [72]. Hence, we generate items with skewed distributions by adding random mean shifts to the columns of \(\boldsymbol{Z}\)—with all items within the same construct receiving the same mean shift—before discretization into \(K \in \{5, 7, 9\}\) answer categories interpreted as values \(\{1, \dots, K\}\). Similar to the previous simulations, the mean shifts are drawn from the interval \([0, s_{\text{max}}]\), where \(s_{\text{max}}\) is chosen to result in 40% of the responses being expected in the highest answer category.

The resulting continuous data matrix is then discretized using equispaced breakpoints \(-K/2 + 1, -K/2 + 2, \dots, K/2 -2, K/2 -1\). To simulate reverse-keyed items, which are commonly used in practice to help detect response inattention, half the items of each construct are reversed by reassigning answer category 1 to \(K\), 2 to \(K-1\), and so on. Note that this implies that the data contains both left-skewed and right-skewed items. As is common in the behavioral sciences, we randomize the order of the items in the survey (with the same order being used for all participants), resulting in the final rating-scale data matrix \(\boldsymbol{R}\).

7.2.2 Missing values and careless responding↩︎

While online data collection tools can be set up to require participants to respond to individual questions in order to progress in the survey, participants may abandon the survey altogether. We therefore inject missing values into the data matrix \(\boldsymbol{R}\) by simulating respondents who stop responding at some point in the survey. We set the proportion of respondents who abandon the survey to \(\gamma \in \{0.2, 0.4, 0.6\}\). For each of those respondents, we randomly select the item at which they stop responding, with all responses from this item onwards being replaced by missing values. We emphasize that here \(\gamma\) does not indicate the proportion of missing cells in \(\boldsymbol{R}\), but the number of rows that contain missing values.

To simulate careless respondents, a proportion \(\varepsilon \in \{0, 0.1, 0.2\}\) of the rows in \(\boldsymbol{R}\) are replaced with careless respondents who randomly select from the two extreme answer categories \(\{1, l\}\). Although this is a common careless response style [54], [73], it may not be realistic that all careless respondents behave in this way. Yet we chose this careless response style because the high leverage of the responses should be particularly influential for matrix completion methods. Note that our simulation allows for the situation that a careless respondent abandons the survey such that a certain row may contain both careless responses and missing values.

7.2.3 Evaluation criteria↩︎

We compare the methods by computing the mean absolute error (MAE) over the predictions of missing cells as in 5 , but excluding missing cells in rows corresponding to careless respondents. The reason for excluding the latter is that in a practical setting, careless respondents should be handled by detection and removal, or better by applying robust methods that downweight such outlying observations [54].

7.2.4 Results↩︎

We focus on the results from Figures 13 and  14 for the settings with 20% of respondents who abandon the survey, as the results for 40% and 60% are similar (see Figures 15–18).

Box plots for the MAE in the simulated survey with \(p=40\) items are shown in Figure 13. In general, differences between the methods are smaller than in the simulations for recommender systems. On the other hand, we now see more differences among RDMC with different loss functions. Most notably, RDMC with the pseudo-Huber function remains close to SI-discretized in the absence of careless responding, while RDMC with the truncated absolute loss yields the best results for 20% careless respondents. In addition, with a higher number of answer categories, robustness benefits of RDMC over SI become more apparent. For instance, with 10% of careless respondents, RDMC and SI-discretized perform similarly for five answer categories, whereas all three variants of RDMC outperform SI-discretized in case of seven or nine answer categories.

Increasing the number of items to \(p=80\), Figure 14 contains the corresponding box plots of the MAE. We now observe bigger differences between the methods. For all settings regarding the number of answer categories, the MAE of RDMC remains in between that of the two SI approaches in case of no careless responding, while RDMC yields better performance in the presence of careless respondents. As previously, the pseudo-Huber loss is preferable in the absence of careless responding (with an MAE close to that of SI-discretized), while the truncated absolute loss is the most robust for larger prevalence of careless respondents.

Finally, we included variants of RDMC with the three loss functions for which we perform at most 10 iterations. Here, the MAE is at best similar but in most cases slightly higher compared to iterating until convergence. For simplicity, we thus omit these variants from the figures.

In summary, the choice of loss function seems to play a more important role in a survey context than in our simulations of recommender systems (see Section 3). The pseudo-Huber loss is thereby recommended if one expects relatively few careless respondents, whereas the truncated absolute loss is recommended if one expects a sizable number of careless respondents.

References↩︎