extfieldsofstudy:
Computer Science
Mathematics

text:
Multi-output Gaussian Processes for Uncertainty-aware Recommender Systems

Recommender systems are often designed based on a collaborative filtering approach, where user preferences are predicted by modelling interactions between users and items. Many common approaches to solve the collaborative filtering task are based on learning representations of users and items, including simple matrix factorization, Gaussian process latent variable models, and neural-network based embeddings. While matrix factorization approaches fail to model nonlinear relations, neural networks can potentially capture such complex relations with unprecedented predictive power and are highly scalable. However, neither of them is able to model predictive uncertainties. In contrast, Gaussian Process based models can generate a predictive distribution, but cannot scale to large amounts of data. In this manuscript, we propose a novel approach combining the representation learning paradigm of collaborative filtering with multi-output Gaussian processes in a joint framework to generate uncertainty-aware recommendations. We introduce an efficient strategy for model training and inference, resulting in a model that scales to very large and sparse datasets and achieves competitive performance in terms of classical metrics quantifying the reconstruction error. In addition to accurately predicting user preferences, our model also provides meaningful uncertainty estimates about that prediction.

INTRODUCTION
Collaborative filtering (CF) provides a powerful solution to recommender systems. Recommending a new item to a * Work done for Siemens AG user is based on the assumption that users demonstrating similar rating or purchasing patterns are interested in similar items. A database describing such user-item interactions often takes the form of a matrix, where each entry describes the interaction between one user and one item. The overall rating or purchasing pattern of a user can therefore be described by the corresponding row in such a matrix. However, since there are typically large numbers of users and items in the database, and each user is usually only interested in a small subset of items, this user-item matrix is often large and sparse. It is therefore inefficient to define the similarity between users in the high dimensional feature space defined by all items. Instead, it is more advantageous to derive abstract feature vectors that represent users and items, which inspired a large variety of low-rank matrix decomposition models such as non-negative matrix decomposition [Zhang et al., 2006], biased matrix decomposition ] and non-parametric decomposition [Yu et al., 2009]. These methods aim at learning low dimensional representations for all users and items, allowing for the prediction of the unobserved interaction between a new pair of user and item.
Most of these state-of-the-art representation learning methods in form of matrix decomposition, however, focus on point estimates and do not quantify any predictive uncertainty. This limits the applicability of these methods in safety-critical situations. For instance, the recommendation of a treatment in the context of clinical decision support needs to be interpreted with much more caution than the recommendation of movies, and transparent uncertainty estimates are crucial to assess the quality of such recommendation. Similarly, when recommender systems are used to detect novel drug-target interactions to facilitate new drug discoveries, it is crucial for such systems to be uncertainty aware, since large costs are associated with studies required for target validation. Other applications could arise in the context of multi-armed bandits for exploration and exploitation during recommendation. Furthermore, missing ratings in collaborative filtering tasks are usually not completely random [Schnabel et al., 2016, Marlin et al., 2012, and good performance in terms of RMSE and MAE does not always promise good recommendations [Cremonesi et al., 2010]. In this context, the predictive uncertainty serves as a practical way to make this bias in the training data transparent in a principled manner. The model can communicate the uncertainty of a prediction to the user, such that the user knows which individual predictions are likely to translate into good recommendations.
Gaussian processes (GP) are a well known class of models that generates prediction uncertainty alongside point estimates via its kernel function, which measures the similarity between data samples. While Gaussian processes are primarily used to solve a wide range of supervised modelling tasks, they have also been applied to perform probabilistic matrix decomposition with a focus on dimensionality reduction, via the Gaussian Process Latent Variable Model (GPLVM). Here, a Gaussian process is used as a prior distribution for a function that maps a low-dimensional latent space to the high-dimensional data matrix. The individual components of this mapping function are modelled as independent draws from a Gaussian process, assuming all dimensions of the data matrix are independent. This independence assumption is problematic in the context of collaborative filtering, where modelling the interdependency of both users and items has proven beneficial [Dong et al., 2014, Nickel et al., 2011. While GPLVM has been used as matrix factorization method for collaborative filtering in the past [Lawrence and Urtasun, 2009], as for other matrix decomposition approaches, practical applicability can be limited also for computational reasons. Being designed primarily as a dimensionality reduction method for data with few missing values, it is not readily amenable to the triple-based learning paradigm of collaborative filtering. That is, one trains a learning model with data samples in form of triples of indices and the corresponding target values, e.g. (i, j, y i,j ). The learning model will have to query the corresponding rows i and j from the respective factor matrices and attempt to regress them to the target entry y i,j . This paradigm facilitates efficient modeling of large and highly sparse datasets such as large useritem matrices or even tensor representation of knowledge graphs [Nickel et al., 2015]. However, the GPLVM-based recommender systems typically require storing the dense data matrix. More specifically, in the triple-based learning paradigm we only need iterate through all observed, i.e., non-zero scalar entries in a sparse matrix.
In this work, we propose a novel approach to combine multioutput Gaussian processes with representation learning for collaborative filtering via matrix decomposition. We learn a latent representation vector for each entity type (i.e. users and items) jointly with a multi-output GP that predicts not only a point estimate of the user-item interaction, but also its predictive variance. We thereby explicitly model not only dependencies between users, but also between items. We motivate this design by connections to coregionalization with separable kernels in multi-output Gaussian processes as well as to the GPLVM. We exploit variational sparse approximations [Titsias, 2009, Hensman et al., 2015 and combine them with a sparse matrix representation so that our approach scales to datasets consisting of large numbers of users and items via minibatch learning. Furthermore, we propose a new approach for evaluating and visualizing the quality of predictive variance of such uncertainty-aware models by stratifying predictions by their associated uncertainty.

RELATED WORK
Matrix decomposition methods applied to collaborative filtering and recommender systems attempt to learn explanatory latent features for users and items. Non negative matrix factorization [Luo et al., 2014, Zhang et al., 2006 restrict the factor matrices to be non-negative, in order to enable better interpretability of the representations. Such methods are trained to minimize the reconstruction error in terms of the entire matrix, and are therefore no longer efficient or reliable in case of high sparsity in the matrix. In contrast, biased matrix factorization [Koren et al., 2009, Gomez-Uribe andHunt, 2015] has proven to be very efficient in decomposing sparse matrices by training only based on observed entries. However, these methods are prone to overfitting and require careful tuning and validation [Nickel et al., 2015].
The Gaussian Process Latent Variable Model (GPLVM) [Lawrence, 2004] applies GP regression in an unsupervised setting. Given a matrix Y , the GPLVM learns a latent representation vector a i for each row y i . A GP regression model is then used to map the unobserved latent variables a i to the multiple output dimensions y ·,j . Since the latent variables are of lower dimensionality than the outputs, GPLVM can be seen as a dimension reduction approach. While GPLVM models are typically applied to dimensionality reduction tasks, [Lawrence and Urtasun, 2009] apply GPLVM to largescale matrices to perform item recommendation by collaborative filtering. While the GPLVM is usually viewed as a matrix factorization approach and has also an interpretation as probabilistic PCA, it can also be interpreted as a multioutput GP regression model where the vectors of different outputs are drawn independently from the same Gaussian process prior.
Multi-output GP regression models can be used to predict vector-valued functions in a supervised manner by constructing a covariance matrix that describes the dependencies between all the input and output variables. Most implementations of multi-output GPs are formulated around the framework of the Linear Model of Coregionalization (LMC) [Alvarez et al., 2012], in the context of supervised learning. Within this framework, dependencies between individual regressors are modeled with separable kernels (or the sum thereof), which can be written as the product of two kernels. The first kernel measures similarity of samples in the input space, while the second kernel captures the similarity between each pair of output dimensions. [Dai et al., 2017] extend multi-output GPs by representing the output dimensions with unobserved latent variables. Thereby they avoid problems with overfitting caused by having to estimate all parameters of the full covariance matrix in standard LMC, allow for the prediction of new outputs at test time and facilitate a lower computational complexity. They further make sure their model scales to large datasets by utilizing the sparse Gaussian process approximation [Titsias, 2009, Titsias andLawrence, 2010] (SVGP). In the SVGP framework, the model is augmented by an auxiliary variable and the covariance matrix is computed on a set of inducing inputs that represent the whole dataset. The kernel distance between two input samples can then be factorized into their respective kernel distances to inducing points, thus avoiding the need to compute a full covariance matrix. For a systematic overview and comparison of sparse approximations, we refer to [Quiñonero-Candela and Rasmussen, 2005] and [Bui et al., 2017].

PRELIMINARY: MATRIX FACTORIZATION FOR COLLABORATIVE FILTERING
In collaborative filtering, the dataset often takes the form of a user-item matrix Y ∈ R I×J , assuming I users and J items. Each entry y i,j describes the interaction between user i and item j, for example in form of a rating or a purchase. Such a matrix is typically large and sparse, first due to the large number of users and items, and second due to the fact that each user is usually interested in a very small subset of items only. In order to derive new recommendations for a user, one has to predict their potential interest in an item with which no interaction has taken place yet, i.e. where there is no entry in the user-item matrix. A typical solution to this collaborative filtering task is via a matrix decomposition approach.
During training, one fits a decomposition model that can recover all entries in the user-item matrix to a certain extent. As a simple illustration, the user-item matrix can be mod- At inference time, any entry value that is not observed in the training data, can be predicted by reading the corresponding location on the reconstructed the user-item matrix: Y = AB T . The generalization power of the decomposition model lies in the low-rank matrices A and B, which can be interpreted as latent representations of the users and items, respectively. In order to achieve better modeling performance, it is possible to apply more complex functions than the dot product to join the latent representations. Writing the model definition as the function g can be realized in various ways. For instance, SVD defines g(a i , b j ) = (λ • a i ) T b j , where λ consists of the singular values and a i and b j are restricted to be orthogonal. In case of biased matrix decomposition [Koren, 2009] where α i and β j are user-and item-specific biases parameters, respectively. Multi-way neural networks [Dong et al., 2014], if applied to the matrix case, take the form of Here, the latent representations are concatenated and fed into a multi-layer perceptron.
The GPLVM [Lawrence, 2004[Lawrence, , 2007 solves the matrix decomposition task by fitting J independent GP regression models on unobserved latent variables x i with p(y ·,j ) = N (y ·,j | 0, K j + σ 2 I), where K j ∈ R I×I captures the covariances between each pair of latent representation (a i , a i ) defined by the covariance function k(a i , a i ). In other words, the GP-LVM can be interpreted as a multiple-output GP regression model where only the output data are given and the unobserved inputs are being optimized.
It is worth noting that the training target in case of standard GPLVM implementations is y ·,j . If the target matrix, and hence each column y ·,j , is sparse, such dense representations are computationally inefficient. In case of (biased) matrix decomposition and multi-way neural networks, however, one could formulate a training sample in form of a triple (i, j, y i,j ) and only iterate through the observed nonzero target values.

NOTATION
We denote a GP regression on target y as where σ is the hyper parameter that defines the random noise level. The vector f denotes the vector consisting of values produced by the regression model from all input data x i , which follows another multi-variate Gaussian distribution: where we denote the covariance matrix using K.

MODEL DEFINITION
The GPLVM treats all output dimensions in the matrix as independent, which is an assumption that may not always hold, in particular if the matrix describes user-item interactions. For instance, the items in which a user takes interest could very well be correlated or similar [Nickel et al., 2015, Agrawal et al., 1994. To capture such dependencies between output dimensions, we propose a new coregionalization kernel to perform multi-output Gaussian Process regression for unsupervised matrix decomposition in the spirit of the GPLVM.
Our proposed separable kernel can be written as the product of two individual kernels, where the first kernel measures similarity of samples in the input space and the second kernel captures the similarity between each pair of output dimensions. More formally, this kernel takes the form or, equivalently, where k A and k W are scalar kernels on A × A and {1, . . . , J} × {1, . . . , J} respectively and is a symmetric and positive semi-definite matrix, which models the dependency between each pair of outputs. While coregionalization is usually performed in the context of supervised regression, in our application the inputs a i are unobserved and, as for the GPLVM model, need to be optimized. W being the identity matrix implies independence between outputs and the model falls back to the standard GPLVM. There is a variety of approaches for choosing W , ranging from the design of a symmetric and positive semidefinite matrix based on suitable regularizers to the choice of covariance functions for the different output components [Alvarez et al., 2012].
Here, we choose to replace the coregionalization matrix W by a kernel on latent representations of items in vector space B, such that the covariance matrix can be written as a Kronecker product of the covariance matrix of the latent variables representing the items K B and the covariance matrix of the latent variables representing the users K A . We thereby combine representation learning with coregionalization in multi-output GPs for collaborative filtering via unsupervised matrix factorization. Note that we have briefly introduced Latent Variable Multiple Output Gaussian Processes (LVMOGP) [Dai et al., 2017] as a related method in Sec. 2. In the context of supervised learning, LVMOGP augments multi-output GPs by latent variables, capturing latent information on the outputs. While the authors also replace W by a kernel on latent representations, in their supervised scenario, inputs are observed, whereas in our model both kernels act on latent variables. In Tab. 1 we compare the most closely related methods to ours, including multi-output GP regression, GPLVM and LVMOGP, in three aspects of learning task, output and input treatment, respectively.
Taken together, the kernel of our proposed model can be written as For the sake of symmetry, we choose the same covariance function for k A and k B . As for matrix K A , each element in K B measures the similarity between a pair of output dimensions (or items) (j, j ). As the rightmost term of the equations above implies, we could treat all latent variables in matrices A and B as hyper parameters in a specific kernel [Snelson and Ghahramani, 2006] that takes as input the indices. From a functional perspective, our proposed kernel measures in fact the similarity between user-item pairs (i, j) and (i , j ). The measurement is carried out by evaluating the two kernel functions at (i, i ) and (j, j ), respectively, and calculating the product of both kernel values. The training samples can be thus formulated using only the indices and corresponding entries in the matrix: J]. This data format is especially advantageous in case of large number of items and high sparsity in matrix Y since it does not require storing the data in dense matrix format.
We finally write our model as: f ∼ N (0, K coreg + σ 2 I) ∈ R I·J , K coreg ∈ R (I·J)×(I·J) . The vector f has length I · J, indicating that it consists of the outcome of all possible combinations of I users and J items.

MODEL FITTING
The major challenge in fitting such a model lies in the size of the covariance matrix K coreg ∈ R (I·J)×(I·J) , since there are typically large numbers of users and items. Computing the inverse of the covariance matrix, which is required for computing the log marginal likelihood, has a complexity of O(n 3 ) where n is the number of training samples, in our case n = I · J. We take sparse GP approaches to address this challenge. A sparse GP model introduces m n inducing points, denoted as Z m that represent the entire dataset. These inducing points are optimized as additional parameters in a GP model. We apply SVGP [Titsias, 2009] and augment the model with inducing variables u with a Gaussian prior p(u), that contain the values of function f at inducing points Z m . With this augmentation, a variational lower bound of the marginal log likelihood can be derived Multi-output GPR GPLVM LVMOGP Ours Unsupervised task × × Dependency between outputs × Coregionalization by kernel × × Unobserved input × × Table 1: A comparison of our methods with mostly related approaches. Coregionalization by kernel refers to learning latent vectors that represent output dimensions. In contrast, standard coregionalization handles dependency between outputs using a full covariance matrix whose only constraint is being symmetric and positive semi-definite. [Hensman et al., 2015]: with q(u) = N (u|µ, Σ) being the variational posterior distribution on the inducing variables u, and Q nn = K nm K −1 mm K mn . K mm is the covariance function evaluated between all the inducing points and K nm is the covariance function between all inducing points and training points. This formulation enjoys the advantage of being able to factorize across data, i.e., to be updated using single samples or mini batches, and can thus scale to large training data.
The essential motivation of applying inducing points is to avoid calculating and inverting the full covariance matrix K ∈ R n×n . Instead, one only calculates the kernel between a training sample and all inducing points as in K mn , and the kernel between two inducing points as in K mm , which is defined as follows in our specific case of coregionalization: SVGP therefore reduces the computational complexity from O((I · J) 3 ) to O((I · J) · (m A · m B ) 2 ). Our triple-based formulation means that the set of training samples consists only of observed triples; by making use of the SVGP approximation, we only ever need to compute the kernel between such a training sample and all inducing points or between pairs of inducing points. This in turn means that for matrices with sparsity S (fraction of observed useritem interactions), the computational complexity reduces to O(S · (I · J) · (m A · m B ) 2 ). In practice, we can further reduce this computational complexity if we choose the same number of inducing points m in A and B: this allows us to learn coupled pairs of inducing points in both spaces (see next subsection for implementation details). By tying the parameters of both sets of latent representations, we can then reduce the effective size of K mm to m × m and further reduce computational complexity to O(S · (I · J) · m 2 ).

COUPLING INDUCING POINTS
The adjacency matrix characterizing interactions between users and items in collaborative filtering tasks is typically highly sparse. An efficient representation of such data is via a triple store where an observed training sample is represented as (i, j, y i,j ) ∀(i, j) ∈ [1, I] × [1, J] with users i, items j and corresponding entries in the user-item matrix y i,j . In this data formulation for sparse data, input i and output j of one training sample are inherently coupled. We would like for this coupling to be also reflected in the sparse approximation and in particular the choice of inducing points, in order to fully take advantage of the benefits of this data representation in terms of computational complexity. In standard formulations of SVGP for multi-output GPs, the covariance matrix of the inducing variables K mm , is computed as are computed on independent sets of inducing points Z A and Z B with kernel functions k A and k B that represent input space A and output space B respectively [Dai et al., 2017]. Similarly, the cross-covariance between training samples and inducing points K nm is constructed as K nm = K A nm ⊗ K B nm , with K A nm ∈ R I×m A and K B nm ∈ R J×m B , where I and J correspond to the number of inputs (users) and outputs (items) as before. K A nm is the cross-covariance between latent inputs A and Z A and K B nm the cross-covariance between latent outputs B and Z B . The covariance between training samples K nn is computed as K nn = K A nn ⊗ K B nn , with K A nn being the covariance matrix of A constructed with k A and K A nn being the covariance matrix of B constructed with k B . However, this formulation does not reflect the sparse triplebased data formulation and results in a computational complexity that is quartic in m (c.f. previous section). If we decide to choose the same number of inducing inputs for Z A and Z B , we can couple the inducing points for pairs of inputs and outputs and reformulate the construction of K nm , K nm and K nn , such that it reflects the coupling between input i and output j of a training sample (i, j). More specifically, when we compute the cross-covariance between the q-th training sample (i, j) * and the l-th pair of inducing points in order to build K nm , we do so by pairing the corresponding latent representations a i and b j and using a product kernel to compute similarities. That is, Using this paired formulation, it becomes immediately clear that becomes substantially smaller, with K paired nm ∈ R (I·J)×m compared to the standard case where the sets of inducing points are treated as independent and K nm ∈ R (I·J)×(m A ·m B ) . We construct the covariance matrix K mm using the same pairing approach and compute the covariance between the o-th and p-th inducing points as That means, using paired inducing points, K mm can be written as elementwise product between K A mm and K B mm : mm . As before, K mm becomes substantially smaller since it is constructed via an elementwise product rather than via the Kronecker product. Finally, we use the same pairing formalism when constructing K nn and compute the covariance between the q-th training sample (i, j) and the l-th training sample (i , j ) as Constructing K nn in this manner rather than using the Kronecker product between K A nn and K B nn is particularly suited for highly sparse data, where the number of training samples N is substantially lower than total number of possible samples I · J. In this case the paired approach directly leads to a natural construction of a small covariance matrix of observed training samples of size N × N † . Taken together, this means that the paired approach reduces the computational complexity of of our method from O(I · J · m 2 A · m 2 B ) to O(S · (I · J) · m 2 ) (with sparsity S, c.f. section 4.2).
We summarize the training algorithm of our approach as follows. Given the number of inducing pairs m and rank r as hyper parameters: We describe our algorithm as online * There is always a bijection / double-indexing relation between (i, j) and q. † Note: the full covariance matrix is not needed during inference. Model training using SVGP merely requires the computation of the trace of Knn Algorithm 1: Training algorithm input :Training N tuples (i, j, y i,j ) where i ∈ [1, I], j ∈ [1, J], and y i,j is the (i, j)-th entry in matrix Y ; output :A, B, Z A , Z B , µ, S and the RBF parameters; Initialize latent representations A ∈ R I×r , B ∈ R J×r , and inducing pairs Z A ∈ R m×r , Z B ∈ R m×r ; while not converged do for each training epoch do for a training sample (i, j, y i,j ) do Calculate the covariance matrix of inducing pairs K mm using Eq. (3); Derive the bijective index q of (i, j); Calculate its covariance values K nm [q, l] with all inducing pairs, using Eq. (2). Update the ELBO in Eq.
(1) w.r.t. a i , b j , z A l , z B l , µ, S and the RBF kernel parameters. end end end training with single samples for a more concise notation as well as to emphasize the scalability of the training. In the step of updating the ELBO, one only needs to iterate once through all inducing pairs. In practice, one would often like to train on batches. The complexity of this case merely becomes the product of batch size and number of inducing pairs. One should also note that in a standard SVGP setting, one updates the kernel parameters and inducing points. In our approach, however, we update additionally the input to the GP which are the latent representations.
We have implemented our method using tensorflow/GPflow [De G. Matthews et al., 2017] and provide an easy-to-use python package at https://github.com/Tuyki/ mogp-decomposition with experimental results. For all experiments, we used an RBF kernel, for both k A and k B . Experiments were run on hardware utilising 4 CPUs and one GeForce GTX TITAN X GPU.

EVALUATION
Due to the fact that the missing entries in CF tasks may be not be completely random, high performance in terms of MAE and RMSE might not always guarantee good recommendations [Marlin et al., 2012, Cremonesi et al., 2010. To this end, we propose an additional metrics to evaluate uncertainty-aware CF models, including ours and GPLVM, by sorting the point estimates based on the predictive variance. Specifically, we denote the ground truth values as (y i ) n i=1 , and the predictive distribution as (f i ∼ N (m i , s 2 i )) n i=1 . We calculate Q quantiles of all the standard deviations (s i ) n i=1 . For each quantile q ∈ { 1 Q , 2 Q , . . . , 1} where the predictive standard deviation s i is smaller than or equal to the q-quantile, we evaluate all point estimates m i against the corresponding ground truth values y i . Formally let, where R could be any evaluation metric defined on a set of pairs of ground truth and point estimate. Each d q is the evaluation of the predictive performance of the model on its top-q · n most confident predictions. Ideally, an uncertaintyaware prediction model is expected to have better performances on test samples about which it demonstrates higher confidence, i.e., lower predictive variance, and vice versa.
If we plot all d q values on y-axis against the corresponding q values on the x-axis, we would expect a monotonously increasing plot line. We name this evaluation approach QPplot, where Q stands for quantile and P for performance.

EXPERIMENTS
We conduct experiments on three well established datasets for evaluating recommender systems: movielens-1m (ML-1M), movielens-10m (ML-10M) and Jester joke recommendations. The task is to predict the rating value of an arbitrary pair of user and item, which can be later interpreted as a score quantifying the interest of the user in the item.
The ML-1M consists of 1000209 ratings of 6040 users on 3706 items, the ML-10M consists of 10000054 ratings of 69878 users on 10677 items, and the jester dataset consists of 1728847 ratings of 50692 users on 140 items. The ML-1M has as ratings only natural numbers between 1 and 5, while ML-10M has rating in the same range but the ratings also include one value between each pair of natural numbers, such as {1, 1.5, 2, 2.5, ..., 5}. The Jester dataset contains arbitrary real numbers between -10 and 10. Our model, based on a GP regression, models all three types of ratings as real values.
We performed a hyper-parameter search especially with regard to the size of the latent vectors and number of inducing points. On a grid search of batch sizes {2 8 , 2 10 , 2 12 , 2 14 , 2 16 , 2 18 }, a relatively large size of 2 16 produced the best performance. We also perform a grid search for the optimal number of inducing points. Among {32, 64, 128, 256, 512, 1024} inducing points, 128 outperforms other choices. This comparison is also illustrated as QP-plot in supplementary materials. We also find out that the model is robust with respect to the size of latent vector. For all the experiments, we use the best performing size of 8, which is a relatively small number compared to, e.g., SMA [Li et al., 2016], which requires a rank of up to 50.
For each dataset we report the performance in terms of MAE and RMSE on the test set in a 5-fold cross validation setting. Quantitative comparisons regarding the ML-1M, RMSE MAE CF-NADE [Zheng et al., 2016] 0.829 -GPLVM [Lawrence and Urtasun, 2009] 0.880 0.644 Sparse FC [Muller et al., 2018] 0.824 -GC [Berg et al., 2017] 0.832 -CWOCF [Lu et al., 2013] 0.958 0.761 LLORMA [Lee et al., 2013] 0   [Muller et al., 2018] 0.769 -GC [Berg et al., 2017] 0.777 -CWOCF [Lu et al., 2013] 0.903 0.701 LLORMA [Lee et al., 2013] 0.822 -SMA [Li et al., 2016] 0.  We also demonstrate performances of our method using the QP-plot in Fig. 1, 2  case of ML-10M, the best reported RMSE lies between our quantile 90% and 100%, the best reported MAE between 90% and 100%. In other words, for those 70% to 90% of test samples for which our model is most confident, it also produces equal or better point estimates than state-of-the-art recommender systems. It is also worth noting that our model outperforms the only uncertainty-aware model, GPLVM, in all quantiles of variances, presumably due to the explicit modeling of output dependency.
Our training time on the ml-1m dataset is 13.4 min (1 CPU, 1 GPU) for all 500 epochs. For the same setting, the biased MF model requires 16.1 min and SVD++ 290.8 min (1 CPU). The uncertainty-aware Bayesian GPLVM has a substantially longer training time of more than 10 hours (16 CPUs, wall time). We perform inference for all tuples in the test set simultaneously. In case of ml-1m, for a test set of size 200K, the average inference time is 0.506 seconds. This demonstrates the computational efficiency of our solution to couple the inducing points (section 4.4). We also realize that, since the kernel evaluation is independent for different pairs of vectors, our model can profit from the parallelization power of GPUs.

CONCLUSION AND DISCUSSION
In this manuscript, we propose a novel uncertainty-aware approach to perform collaborative filtering in the context of recommender systems. To this end, we introduce a coregionalization kernel to perform multi-output Gaussian Process regression for unsupervised matrix decomposition in the spirit of GPLVM. In contrast to the GPLVM we do not treat output dimensions as independent, but model a coregionalization matrix with a kernel on latent variables representing the output. Our approach not only accounts for dependencies between output dimensions, but also allows our model to handle extremely sparse data in form of triple-stores in an efficient manner.
Finally, we would like to stress that the aim of our method is not to improve prediction quality that can be measured by metrics such as RMSE or MAE. Instead, we are interested in exploring the possible solutions to add uncertainty awareness to collaborative filtering approaches and recommender Figure 2: QP plots regarding the ML-10M dataset, comparing our method with Bayesian GPLVM and best reported performances. The horizontal dashed lines correspond to reported state-of-the-art methods. RMSE: [Zheng et al., 2016].
In terms of MAE, the Biased MF outperforms our method only at quantile 100%. The GPLVM seems to suffer from the data volume and fails to converge in training. systems. More specifically, we enable these methods to express its uncertainty about every prediction. With the help of our evaluation method, the QP-plot, we can show that such uncertainty quantification are valid: for predictions that deviate from the ground truth, the corresponding predictive variances tend to be larger. With our experiments, we also show that Bayesian GPLVM fails to scale to large datasets in terms of both training and inference. Our approach in contrast merely requires training and inference time that is comparable to highly optimized linear models such as SVD++.


Intelligent Diagnosis of Cervical Cancer Based on Data Mining Algorithm

The intelligent diagnosis of cervical cancer by using a class of data mining algorithms has important practical significance. In particular, the useful information included in a significant quantity of medical data may not only discreetly boost the development of medical technology but also detect cervical cancer in the future. This paper improves the data mining algorithm and combines image recognition technology and data mining technology to extract and analyze image features. Moreover, this paper makes full use of the information contained in the image to realize the segmentation of the cervical cancer cell image, select the feature vector according to the characteristics of the cervical cancer cell, and use the statistical classification method to design the classifier. The test results show that the automatic recognition effect of this system is good, and it has a good auxiliary diagnosis effect. Therefore, it can be verified in clinical practice in the follow-up.

Introduction
Various kinds of malignant lesions may develop in different areas of the body, among the illnesses that pose a danger to people's lives. Cervical cancer is the fourth highest cause of death in women, and it is something that women have long been worried about and should not be disregarded. Cervical cancer therapy is generally determined by the patient's tumor stage. The early stages of cancer are treatable with surgery if detected early enough [1]. Radical surgery for cervical cancer coupled with radiation and chemotherapy is the sole option if the tumor has progressed to an advanced stage and most lymph nodes have been metastasized or invaded [2]. However, the majority of women who are diagnosed with breast cancer are already in the late stages of the disease. By determining the stage of a patient's cervical cancer, doctors can better anticipate the likelihood of recurrence and tailor treatment accordingly. As a result, predictive biomarkers are playing an increasingly significant role in determining prognosis [3].
With the rapid development of technology, our ultimate goal is to mine the information contained within the data.
The application fields of related data mining algorithms for processing big data are very wide, including different fields such as medical treatment, medicine, business activities, science, and social practice networks. This has a lot to do with what we do on a daily basis. To be sure, data mining algorithms that analyze large amounts of big data are intertwined with everything from online appointment registration in hospitals to intelligent illness diagnosis to traffic congestion management and catastrophe prediction. In addition to the enormous quantity, diversity of kinds, and fast updating, these data also represent the value that may be reaped from the data itself.
Data mining has slowly built up theoretical expertise over decades, such as classification and grouping. Among data mining's main areas of study, the classification issue entails creating a classification model using associated algorithms to evaluate and forecast categories of unknowns. Data sample attributes are used as inputs, and the category in which a data sample is classified is returned as an output. The categorization issue is a fictitious one created to distract attention from real-world issues. Cervical cancer mortality and morbidity rates are rising fast throughout the world, and the causes for this are many. While on the one hand reflecting the issue of an ageing population, it also represents some of the most important cancer risk factors. As an alternative, a growing number of nations are experiencing fast population growth and an ageing population, with cancer emerging as the leading cause of death. As the number of patients is increasing year by year, the complexity of cervical cancer lesions and the many uncertainties of suffering from cervical cancer and other factors have led to their clinical data showing a large number of clinical data, diversified data types, extremely fast development, and implicit features of high information value. Patients may be reluctant to provide some private information in relation to clinical data, resulting in the lack of specific clinical data; in other cases, the obtained clinical data contains noise and complex information. The abovementioned issues may be resolved by using conventional techniques. The restrictions are very severe, and the final result of solving the problem is not ideal. Therefore, the use of a type of data mining algorithm for the intelligent diagnosis of cervical cancer has important practical significance, especially the valuable information contained in a large amount of medical data cancer diagnosis promotes.

Related Work
Multiple logistic regression was used in the literature to determine the relationship between late diagnosis and Medicaid status before the disease was even identified [4]. The research examines whether Medicaid recipients are more likely than women who have not participated to be diagnosed with advanced cervical cancer. According to the findings, advanced illness was seen in 51% of Medicaid recipients and 42% of non-Medicaid recipients. Women from lower socioeconomic backgrounds and the elderly are more vulnerable. According to the data shown above, in order to guarantee that at-risk women get screening services, greater awareness is required [5]. There are about 370,000 new instances of cervical cancer each year, according to the research [6]. This accounts for approximately 10% of all new cases of cancer in the world. The first step in combating global cancer is to gather as much data as possible on its prevalence and death. Recent worldwide statistics [7] show that outcomes are comparable to those from 10 years ago, although there is some fluctuation in the figures. Regular and prompt examination inspection and treatment are very important. At the molecular level, the literature [8] studied the pathogenesis of cervical cancer through the implicit data of gene chips, explored malignant tumor markers, and provided a powerful solution for tumor prevention and treatment. The literature [9] investigated the practical usefulness of transvaginal real-time ultrasound elastography technology in the detection of cervical cancer. The combined use of tumor markers was addressed in the literature [10] for the diagnosis of cervical cancer. The effects of particular anti-HPV16E6 ribozymes on cervical cancer cell phenotypic and gene expression have been addressed in the literature [11].
As we all know, for women over 30 years of age, this is the best choice for average risk [12]. The literature [13] evaluates the practices of primary care providers, stimulus factors, barriers to the use of common testing methods, and prolonged screening intervals for low-income women. Few of the patients provided in the literature [14] performed the combined test of the two screening methods, and it is recommended to extend the time interval appropriately. Among them, the excessively frequent screening and the error caused by the misreporting harm are all within the budget allowable range, thereby balancing the obstacles caused by extending the screening interval. Cervical cancer, the second-most incidence of female malignant tumors, is estimated by the World Health Organization. The probability that cancer is the first or second cause of death before the age of 70 is about 52.91%, and it is the third cause of death. The probability of four is about 12.79%. Finally, the remaining 59 countries (172 countries in total) rank cancer as the fifth to tenth [15]. The literature [16] proposed a detection method for diagnosing cervical cancer in women or identifying women who are sensitive to cervical cancer. By detecting Brn-3amRNA or Brn-3a protein (quality) in Brn-3a cells, it can then test whether women are sick or not or sensitive to it. Endometrial cancer and ovarian cancer were studied in genome-wide association research [17] to determine the features of the data composition. The literature [18] examines the diagnosis of ICC every six years in women with human immunodeficiency virus in order to determine the frequency of cervical cancer and compare it to women who are not infected with HIV. Women were the focus of a research performed in August. The research looked at whether or not the Pap test was infected with HIV.

Improved Application of Data Mining Algorithm in Intelligent Diagnosis of Cervical Cancer
The proportional difference model is often used in survival data where the mortality of different groups of patients converges over time. If it is assumed that there are two existing patient treatment groups (i = 1) and control groups (i = 2), the PO model assumes that the difference ratio between these two groups is a constant c: Among them, F 1 ðtÞ represents the cumulative distribution function of the treatment group, and F 2 ðtÞ represents the cumulative distribution function of the control group. This model can be extended: we set x to represent the p -dimensional covariate vector and assume that the relationship between the constant c and the covariate is c = exp ðβ ′ xÞ. Among them, β is the p-dimensional coefficient corresponding to the covariate. Therefore, until time t, the following relationship holds between the difference in the case of covariate x ≠ 0 and the difference in the case of covariate x = 0:

2
Computational and Mathematical Methods in Medicine Among them, Fðt | xÞ and F 0 ðtÞ, respectively, represent the unknown cumulative distribution function and the benchmark cumulative distribution function related to the failure time when the covariate x ≠ 0 and the covariate x = 0. Formula (2) can also be equivalently written in the following form [19]: Among them, ∧ 0 ðtÞ = log fF 0 ðtÞ/ð1 − F 0 ðtÞÞg represents the benchmark logarithmic difference function up to time t , which is the nonparametric function part to be approximated by Bernstein polynomials in this article. It is easy to know from its expression: ∧ 0 ðtÞ is monotonous and nondecreasing with respect to t, and there is ∧ 0 ðtÞ = −∞ and ∧ 0 ð tÞ = ∞. The proportional difference model has a reasonable biological explanation: the jth element β j in the coefficient vector β can be interpreted as the unit contribution of the jth covariate to the logarithmic difference in failure time until time t. The biggest advantage of the score (logit> connection function is that the constant difference rate can be obtained from the model. For example, there are two different individuals, and the covariates are x 1 and x 2 , respectively. From (3), it can be seen that the following expression holds: This means that if the covariates between two individuals differ by one unit, then their logarithmic difference ratios differ by coefficient β, which can also be used as a further explanation for equation (3).
In addition, here, we point out the relationship between the logarithmic difference model and logistic regression, which is convenient for the calculation using the normal data enhancement method in the following text. Formula (2) can also be rewritten as follows: Here, G represents the cumulative distribution function of the standard logistic random variable.
Unlike the Cox proportional hazard (PH) model, the PO model must estimate the nonparametric function AoðtÞ while estimating the parameter β, which makes the parameter estimation problem under the PO model very complicated.
This theorem is at the heart of Bayesian analysis. The posterior probability of an occurrence is described by Bayes' theorem based on a priori information. Let us say we know that age has anything to do with the incidence of a certain malignancy, the factor of age can help people more accurately determine the probability of a person having cancer, that is, the role of prior information. The mathematical expression of Bayes' theorem is Among them, A and B both represent events, and the probability of occurrence of event B is PðBÞ ≠ 0. PðA \ BÞ represents the conditional probability of event A occurring under the condition of event B, and PðB \ AÞ represents the conditional probability of event B occurring under the condition of event A. PðAÞ and PðBÞ are the probability of event A and event B occurring independently, also called conditional probability.
Bayesian inference is based on Bayes' theorem. This method updates the posterior probability distribution of the parameters based on the prior probability and therefore incorporates both the prior hypothesis and a random sample to better serve the goal of statistical inference. More specifically, we use πðθÞ to denote the prior distribution of the parameter θ of interest. The prior distribution of parameters can be any distribution in the parameter space, usually given by experience. W represents the observed sample, so the posterior distribution of parameter θ is the distribution of θ = f ðw | θÞπðθÞ under the conditions given by sample W = w, denoted as πðθ | wÞ. For the case of the continuity density function, the density function formula of posterior distribution is [20]: Among them, hðw, θÞ = f ðw | θÞπðθÞ is the joint distribution of w and θ, and there is which is the marginal distribution of sample w. The posterior distribution function calculated in formula (7) combines the information given by the population (that is, the density function f ), the sample and the prior. From the Bayesian point of view, all statistical inferences should be derived from the posterior distribution function πðθ | wÞ, such as calculating the mean value and confidence interval of the parameters.
Model (2) is combined, and the likelihood function expression based on right-censored data is Among them, δ i1 , δ i2 , and δ i3 are all indicative functions. When the value is 1, it means that the ith individual is left censored, right censored, and interval censored and satisfies δ i1 + δ i2 + δ i3 = 1. Data pre-processing includes : missing value processing, unbalanced data processing, continuous attribute discretization, etc; The basic task of feature extraction of data is usually to basically achieve targeting and obtain the corresponding characteristics under certain constraints; Establish a classifier model for the relevant information within the known data;

Computational and Mathematical Methods in Medicine
This article defines the ith individual: Therefore, based on the relationship between the PO model and logistic regression described in equation (5), this paper uses the following data enhancement methods: Among them, KS represents the distribution function of Kolmogorov-Smirnov distribution, and the specific expression is In the above data enhancement method, the logistic distribution of ∈ i can be derived from the scale normal mixed distribution with the Kolmogorov-Smirnov distribution. Therefore, under the transformation of formula (11), the likelihood function in formula (9) can be rewritten as [21] Among them, there is The main part of the above enhanced likelihood function expression (13) is the product of the constrained normal density, which is very beneficial to the Gibbs sampling method given below.
This article defines all parameters a priori in order to enhance the Bayesian approach by ensuring that the parameters not only offer adequate modeling flexibility but also enable efficient posterior computations. More specifically, this article specifies a multivariate normal prior Nðβ 0 , ∑ 0 Þ for the regression parameter β, a normal prior Nðm 0 , v −1 0 Þ for the coefficient γ 0 of the Bernstein approximation polynomial, and an independent exponential prior ExpðηÞ for the remaining Bernstein coefficients γ i ði = 1, ⋯, NÞ. Exponential priors can shrink relatively small coefficients to zero and thus can play a role in penalizing nonzero coefficients. In order to ensure sufficient modeling flexibility at the same The initial values of all parameters are sampled from their prior distributions. The iterative steps of Gibbs sampling in this paper are as follows: (1) For i = 1, ⋯, n, the algorithm extracts the latent variable z i .
(a) If there is δ i1 = 1, the algorithm extracts samples (2) The algorithm extracts samples of γ 0 from the normal distribution NðE 0 , W −1 0 Þ; among them, there is For k = 1, ⋯, N, the algorithm draws samples of γ k . We If there is W k = 0, the algorithm extracts γ k from ExpðηÞ; otherwise, the algorithm extracts γ k from NðE 0 , W −1 k ÞI ðy k >d * k Þ . Among them, there is The algorithm samples β from the normal distribution Nðβ,∑Þ; among them, there is The rejection sampling algorithm l is used to sample each λ i , and the fully conditional distribution of λ i is proportional to the following formula: The algorithm samples η from distribution Gaða η + k, b η + ∑ N k=1 γ k Þ.

Intelligent Diagnosis of Cervical Cancer Based on Data Mining
When traditional data analysis methods deal with similar data, the effect is not ideal, and sometimes, even a lot of information contained in the data is lost, so that the inherent information and laws of the data cannot be mined. As a result, it is critical to understand the data processing and classification processes, as well as to extract the information hidden in the sample data. Figure 1 depicts the processing of the cervical cancer classification issue in this paper. Discretization of continuous attributes is the second step of data processing in this article. The basic process is shown in Figure 2.
The simplest discretization algorithm-binning technique is utilised in this article. Equal frequency binning and equal width binning are two common binning methods. Unsupervised discretization is used in the binning technique. There are many bins created when sorting an attribute data set from top to bottom. The technique of detection has nothing to do with the procedure. It is possible to get an equal-width binned dataset by using the equal-width binning technique, which uses a fixed width for each bin. To split the data into The cervical cancer cell image automatic analysis system can be represented by the frame structure shown in Figure 3.
Image registration is mainly performed before other image processing steps. For example, it observes and com-pares satellite remote sensing images and medical images obtained under different modalities after registration. In addition, for example, it can obtain the topographic map of the vast land, or observe the lesion area of the medical image more conveniently, and register and fuse the CT and MRI images. One image is designated as a reference image, while the other as a floating image, when two pictures need to be registered. To accomplish the objective of matching the coordinates of the reference picture, the latter must employ a comparable transformation model. The first transformation model is generated based on the difference between the two pictures, and the initial transformation image and the reference image are compared for similarity. As long as the anticipated criteria are not being fulfilled, optimization algorithms are used to fine-tune the transformation model parameters. This procedure is repeated until the similarity measure falls within an acceptable range, at which point the registration process is complete. The differences between pictures of various modalities are primarily a limited range of translation and rotation, according to the features of the    Figure 4 shows how to determine a stiff transformation model, let us get into it.
The bilinear interpolation method uses the gray levels of the four neighboring pixels of the pixel to be found to perform linear interpolation in two directions, as shown in Figure 5.
For the removal of shadow areas, this paper considers the classification of colors based on HSI and Lab space, so white light image samples need to be added to participate in diagnosis. The specific algorithm flow is shown in Figure 6.
On the basis of the above research, the intelligent algorithm of this paper is used to verify the performance. Feature extraction is based on a certain pixel, pixel area, or specific target area as a unit and extracts a certain type of characteristic shown by the gray value of the area, usually expressed as a number or a vector. These values or vectors are used to represent these various kinds of pixel point sets in order to make it easier for the computer to comprehend them. Features are a necessary component of any computer classification method. For various pictures, classification algorithms, and classification purposes, the features that need to be extracted are often different. In feature extraction, people tend to extract too many features. This will result in an increase in computations and superfluous characteristics, which will have a detrimental effect on categorization. As a consequence, the final classification outcome will be influenced by the quantity and efficiency of feature extraction and make a significant contribution. They prefer to use basic and refined characteristics for categories that are not as complex. There is a clear contrast in the gray level between the lesion and the surrounding normal tissue, as seen in the multiband picture above. Because the lesion region appears darker, the average gray value of the sample area was used as the first feature to be retrieved in this study. As for the classification algorithm, due to the small number of features, this paper uses the more commonly used K-means clustering algorithm to classify the surface of the cervical tissue.
Since the system is based on a combination of both the local area network and the Internet, the internal systems use the local area network for data exchange, and the system and users use the Internet for information transmission. As a result, inside the LAN, it is preferable to utilise a star topology logic structure (as illustrated in Figure 7). The reason is that the star topology logical structure is an online structure in which the server (switch) is the centre, and all computers are connected to the server (switch). This structure can prevent each computer from affecting each other when a failure occurs and can work independently, and the real wiring of this star-shaped structure will not be very difficult, as long as there is enough connection medium.
In terms of connection of form data, the connection between forms can be divided into two types. One is to use ADODC, the most common connection method for forms, which is the connection method used for form connection. The page link is generated using a relative route as a hyperlink in this method. As long as the link target file exists in the folder, the hyperlink may be generated using this method without thinking about the computer's path. As a result, when the location of the system folder changes, the hyperlink route does not have to be adjusted. Using data transfer between dynamic web pages as a connection mechanism. Programming is used to make this happen. The linking method can be based on the needs of the editor, what data and how much data is sent to that page, and what operations are performed on the data. In fact, this kind of link is a way of linking data. The difficulty of this link method is that because it is derived from programming, the correctness of the compiled program must be ensured. If the correctness is not ensured, the link to the web page will not work correctly.
The manual module process in this system includes evaluation of fresh image information, change of picture information, and review after picture processing, as illustrated in Figure 8.
After constructing the above system model, the performance of the system is verified, and the intelligent diagnosis effect of cervical cancer is evaluated based on the actual situation. This article primarily performs research from three aspects in the system evaluation: picture feature identification, diagnostic and treatment accuracy, and user experience. The results are shown in Table 1 and Figure 9.
From the above research, the cervical cancer intelligent diagnosis system based on the data mining algorithm proposed in this paper has a good auxiliary diagnosis effect, which can be verified in clinical practice in the follow-up.

Conclusion
With the advancement of information technology, computer systems have been more widely utilised in the service sector, particularly in hospitals. The Gynecology Research Office of the hospital is mainly responsible for judging the pictures of cervical cancer cells, classifying the results of judging cells, and judging diseased cancer cells, normal cells, and  Computational and Mathematical Methods in Medicine suspected cells. However, these operations must rely on image processing software for processing. For front and back ends, the data centre uses database software and system development tools to actualize them. The business processing centre will be kept in a different environment. To avoid severe mistakes due to catastrophic failures, the business and data processing environment is protected in this way. Use of image data to segment cervical cancer cell images, picking feature vectors based on cell characteristic information, and development of a classifier using statistical methods are all discussed in this article. According to the results of the tests, this system has a strong automatic recognition effect and a good supplementary diagnostic effect. Therefore, clinical practice in follow-up may verify it.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare that they have no conflicts of interest.