{'fraction_non_alphanumeric': 0.03724920802534319, 'fraction_numerical': 0.01438753959873284, 'mean_word_length': 4.6488219504920965, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Relative Pairwise Relationship Constrained Non-negative Matrix Factorisation

Non-negative Matrix Factorisation (NMF) has been extensively used in machine learning and data analytics applications. Most existing variations of NMF only consider how each row/column vector of factorised matrices should be shaped, and ignore the relationship among pairwise rows or columns. In many cases, such pairwise relationship enables better factorisation, for example, image clustering and recommender systems. In this paper, we propose an algorithm named, Relative Pairwise Relationship constrained Non-negative Matrix Factorisation (RPR-NMF), which places constraints over relative pairwise distances amongst features by imposing penalties in a triplet form. Two distance measures, squared Euclidean distance and Symmetric divergence, are used, and exponential and hinge loss penalties are adopted for the two measures respectively. It is well known that the so-called"multiplicative update rules"result in a much faster convergence than gradient descend for matrix factorisation. However, applying such update rules to RPR-NMF and also proving its convergence is not straightforward. Thus, we use reasonable approximations to relax the complexity brought by the penalties, which are practically verified. Experiments on both synthetic datasets and real datasets demonstrate that our algorithms have advantages on gaining close approximation, satisfying a high proportion of expected constraints, and achieving superior performance compared with other algorithms.

A seminal approach in NMF is the so-called "multiplicative update rules" which guarantees both the convergence of the algorithm and the non-negativity of factorised matrices [10]. Though the "multiplicative update rules" were proved not converging to a stationary point numerically [11], and they are not strictly well-defined because of possible zero entries [12], it practically produces satisfactory results, especially for large scale data, which makes it a popular solution for NMF. However, the original NMF only imposes the non-negativity constraints on both of the factorising and factorised matrices, which in practice may not be enough to satisfy additional requirements. Thus, researchers in this area have been proposing new algorithms under this framework to cater for incremental improvements, variations, and/or application oriented constraints [13], [14], [15], [16], [17], [18].
A main sub category of NMF is Constrained Non-negative Matrix Factorisation (CNMF), which imposes constraints based on variables as regularisation terms [19]. The most commonly used regularisations for NMF are L1 norm and L2 norm, the former increases the sparseness of the factorised matrices, while the latter makes the results smooth to prevent overfitting. However, these constraints are only imposed on each of the rows or columns and has not considered the relationship among pairwise rows or columns. Such pairwise relationships exists wildly in many systems, especially in those systems where the factorised matrices are representing features. A typical instance is the matrix factorisation technique used in recommender systems. In a recommender system, the factorising matrix is usually the rating matrix whose entries denote the ratings given by the corresponding user (row) to the corresponding item (column), and the factorised matrices are usually regarded as the user feature matrix (the left factorised matrix) and the item feature matrix (the right factorised matrix). Fig. 1 shows an example of a simple movie recommender system. In this example, after factorisation by NMF, the distance between "Star Wars" and "Titanic" becomes less than the distance between "Star Wars" and "Star Trek". However, as we all know, the movie "Star Wars" should be closer to the movie "Star Trek" -both are within the scientific fiction genre, unlike "Titanic", which is a love story. Thus intuitively, it will result in a better factorisation and make better recommendations if we can incorporate such human-aware relative relationships into the NMF model.
There have been a few studies that attempt to consider the above pairwise relationship in the area of matrix factorisation. Although none of these have properly addressed our problem, two existing methods are still worth noting here: one is Graph Regularised Non-negative Matrix Factorisation (GNMF) [20], the other is Label Constrained Non-negative Matrix Factorisation (LCNMF) [21].
GNMF constructs a weight matrix of the graph from the observed data, and then applies the weights (similarities) on the factorised low-dimensional data representation as regularisations. arXiv:1803.02218v1 [cs.LG] 5 Mar 2018 Fig. 1: A simple movie recommender system. The factorisation algorithm used here is NMF with Euclidean measure proposed in [10]. Original missing ratings are denoted by 0 and recovered ratings are showed in brackets. As shown under the right factorised matrix, the distance between feature vectors of movies "Star Wars" and "Titanic" is less than the distance between movies "Star Wars" and "Star Trek".
It is designed as a dimensionality reduction method, and it works well on image clustering applications where the data points in different classes are distinctively different. However, in many other cases, such as where the data points are not spread and where the matrix factorisation is not used to reduce dimensionality (like in recommender systems), GNMF cannot guarantee the relative relationship denoted by similarities retained as expected after factorisation. Besides, the setting of similarities is sensitive to the factorisation results, especially when the similarities cannot be simply set zeros and ones, such as when there exists chain constraints.
LCNMF was proposed to cater for scenarios where partially labeled grouping data was made available. Its key idea is, if two feature vectors are labeled into the same class, they are assumed to have the same feature representation in the latent space. This approach has addressed the need of applications on image clustering, however, such a setting is far too restrictive in general: "Star Wars" and "Star Trek" could be very similar to each other, but setting their features identical is unacceptable and impractical.
In this paper, we propose a novel matrix factorisation algorithm, called RPR-NMF. Rather than using explicit similarities or previously known labels, RPR-NMF imposes penalties for relative pairwise relationships (RPRs) in a triplet form. The penalties are not limited to be within [0, 1] as for similarities or to be binary values as for labels. Both of the squared Euclidean distance and the symmetric divergence measure are used in the objective of RPR-NMF, and the penalties are in exponential and hinge loss forms respectively. The update rules for RPR-NMF conform to the well-known "multiplicative update rules" in which the proofs of convergence are essential for the whole algorithm. Due to the complexity of proof brought by the imposed penalties, we approximate partial terms in the proof part and have verified its practical benefits through numerous experiments.
Compared with the existing methods GNMF and LCNMF, RPR-NMF can guarantee more pairwise relationships retrained after factorisation. Fig. 2 gives a demonstration of the RPRs among four data points after running GNMF, LCNMF, and our proposed algorithm RPR-NMF using Euclidean measure respectively. In this example, GNMF failed on retaining one RPR (points 2 and 4 should be closer than points 3 and 4 ), LCNMF projected all points onto one, while RPR-NMF retains all RPRs after factorisation.
The main contributions of this paper are: 1. We propose a novel algorithm named RPR-NMF, which utilises relative pairwise relationship among rows or columns of factorised matrices to achieve a better factorisation with high constraint satisfied rate and close approximation simultaneously.
Different from GNMF and LCNMF, RPR-NMF can guarantee the expected RPRs retained after factorisation and does not limit the feature vectors to be identical.
2. In the part of method for RPR-NMF, we use different forms of penalties for Euclidean measure and Divergence measure based on our observations on method convergence through numerous experiments. For the Euclidean measure, we incorporate the RPRs in natural exponential functions, while for the Divergence measure, we use hinge loss function.
3. The solution of RPR-NMF conforms to that of "multiplicative rules", in which we provide complete and sufficient proofs for both of the distance measures with the help of relaxation on partial terms. Such relaxation has been verified reasonable and practical through our experiments.
4. The complexity analysis shows that RPR-NMF does not increase much of processing time by introducing penalty terms when comparing to NMF. Synthetic and real datasets experiments both demonstrate that RPR-NMF have advantages on close approximation, high constraint satisfied rate and outstanding application performance.
The rest of this paper is organised as follows: in the next section, we review related literature to our work. The Method section contains the details of our algorithm and the proof of convergence, followed by the Experiments section in which we evaluate our algorithm as well as the Conclusions.

RELATED WORK
NMF was first proposed to solve the following optimisation problem: given a non-negative matrix V , find non-negative matrix factors W and H such that V ≈ W H. Since the seminal work of [10] which has proposed the so-called, "multiplicative update rules" for non-negative matrix factorisation, a number of approaches followed suit dealing with various NMF issues from different aspects. Our work falls under the category of the Constrained Non-negative Matrix Factorisation (CNMF) which was first proposed in [19]. In general, it has the following representation: for W 0 and H 0, where · F is the Frobenius norm, α and β are regularisation coefficients. The functions of J 1 (W ) and J 2 (H) are penalty terms used to enforce certain constraints on the solution of Eq.(1). In their work, the penalties are set J 1 (W ) = W 2 F and J 2 (H) = H 2 F in order to enforce the smoothness in W and H respectively.
Many studies followed the above CNMF framework. [22] used an adaptive potential function as penalties to characterise the piecewise smoothness of spectral data. [23] imposed three additional constraints on the NMF basis to reveal local features. [24] imposed both L1 and L2 norms to control the shape of base matrix and increase the sparseness of the coefficient matrix so as to enhance the clustering performance on multiple manifolds.
Besides the above studies based on CNMF that consider how to constrain the value of each vector in factorised matrices, there are two algorithms that take the relationship among vectors of factorised matrix into consideration: Graph Regularised Nonnegative Matrix Factorisation (GNMF) proposed by [20], and Labelled Constrained Non-negative Matrix Factorisation (LCNMF) proposed by [21].
GNMF, also a CNMF algorithm, is to utilise the relationship among factorised rows or columns for a dimensionality reduction issue. It is in an effort to ensure that similarities between data points at the original space are also retained after they are transformed to the low dimensional subspace through factorisation. Its objective function with Euclidean measure is as following: where L is a graph Laplacian matrix obtained from the similarity matrix and λ is a regularisation coefficient. As showed in the objective function, GNMF tries to minimise two parts simultaneously: the squared errors between the product of factorised matrices and the factorising matrix, and the similarity matrix. The ideal solution is to minimise them at the same time. However, it often happens that if GNMF minimised the squared error, the RPRs implied by similarities might not be satisfied. As shown in Fig. 2, the RPR among points 2 , 3 and 4 was opposite to that when they were in the original 3D space. Our proposed RPR-NMF works in a different way: it imposes penalties with respect to the expected RPRs, which forces the factorised feature vectors to keep as many of the expected RPRs as possible. LCNMF also utilises the relationship among factorised vectors. Instead of imposing regularised constraints, it represents the relationship by altering the factorising structure. Thus it is not a CNMF method. LCNMF uses partial label information as hard constraints and turns the original NMF task into a semi-supervised problem: they represent the right factorised matrix by a product of a class matrix and a reduced feature matrix where the class matrix contains binary entries to divide data into a predefined number of classes. Its objective function with Euclidean measure is as: where Y is the reduced feature matrix and B is the class matrix. The method however, assumes that if two data points have the same label, their corresponding feature vectors must be identical, as showed in Fig. 2 where all four data points are labelled the same and are projected onto one point after factorisation. Such constraints are too restrictive under many general settings. For example, if two movies are by the same director, in the same genre and even feature the same actors, setting their features identical is to ignore any difference between them, which is what our method aims to mitigate.

METHOD
In this section, we introduce a new factorisation algorithm when RPR constraints are in place, called RPR-NMF. Note that both GNMF and LCNMF only impose constraints on the right factorised matrix because they were proposed as data dimensionality reduction methods. For generality, RPR-NMF imposes constraints on both factorised matrices, and it is trivial to only impose constraints on one factorised matrix. Consider a dataset represented by a non-negative N × M matrix V . This matrix is then approximately factorised into an N × K matrix W and a K × M matrix H, where K is usually set to be smaller than both N and M , and is commonly referred to as the latent dimension. The RPR constraints placed on the factorised matrices can be defined as two sets of integer indexed triples: Specifically, each triple represents the relative relationship between two pairs of vectors with one sharing vector. In our work, if W was the matrix in question and l th triple specified the distance between vector q and r to be less than the distance between vectors q and s, the relationship could be denoted as dis(W q l : , W r l : ) < dis(W q l : , W s l : ) where W q: is the q th row vector of matrix W , and dis(x, y) measures the distance between vectors x and y. We follow the most commonly used two distance measures in our work, which are the squared Euclidean distance and the Divergence (when the variables are two unit vectors/distributions, it becomes KL-Divergence) Since the Divergence of two vectors is not symmetric (D(x||y) = D(y||x)), when characterising the imposed RPR constraints, we use the Symmetric Divergence defined as Then the constraints are incorporated as penalty terms in the objective function. In our work, the penalty format for Euclidean measure is of an addition of natural exponential functions, while that for Divergence measure is of a hinge loss function. The reason for not using the same format of penalties is that the Divergence measure with exponential penalties cannot guarantee a high proportion of satisfied constrains after factorisation, and that the Euclidean measure with hinge loss penalties cannot steadily converge. Thus with two independent coefficients λ W and λ H , we define the objective function using Euclidean distance as [exp(E(W q l : , W r l : )) + exp(−E(W q l : , W s l : ))] and the objective function using Divergence as where l W and l H are the numbers of constraints.

Solving Objective Functions
To solve the above objective functions, we need to derive the update rules for W and H. As a matter of fact, the penalties are not convex even when fixing one of the matrix factors. However, we found it is still feasible to obtain the update rules by constructing and solving the corresponding Lagrange functions. Once we obtained the updating rules, we could minimise the objective functions by iteratively updating W and H.

Updating rules for Euclidean measure
As for the objective function in Eq.(9), we first construct a Lagrange function with non-negative constraints W ij 0 and H ij 0: [exp(E(H :q l , H :r l )) + exp(−E(H :q l , H :s l ))].
(11) The partial derivative with respect to W ab is: where Let the partial derivative vanish and considering the nonnegative constraints (α ab W ab = 0 under K.K.T conditions), we obtain: thus the update rule for W ab is formulated as following: Similarly, we have the update rule for H ab as As for the update rules, we have the following theorem: (9) is non-increasing under the update rules in Eq. (15) and Eq. (16) with appropriate penalty coefficients λ W and λ H .
We provide the prove of the convergence for the above updating rules and Theorem 1 in Section 3.2.1.

Updating rules for Divergence measure
As for the objective function in Eq.(10), we first construct a Lagrange function with non-negative constraints W ij 0 and H ij 0: The partial derivative with respect to W ab is: where Let the partial derivative equal to zero as well as considering the non-negative constraints (α ab W ab = 0 under K.K.T conditions), we obtain: thus the update rule for W ab is formulated as following: Similarly, we have the update rule for H ab as Notice that the value of P row (W ab ) and P col (H ab ) may be negative during updating. Thus if the denominator in an iteration is less than zero, then we abandon the penalty parts. Besides, we dynamically change the penalty coefficients to ensure the convergence since the update rules are obtained by approximation (see proof of convergence in Section 3.2.2). As for the update rules, we have the following theorem: The objective function F 2 in Eq. (10) is nonincreasing under the update rules in Eq.(23) and Eq. (24) with appropriate penalty coefficients λ W and λ H .
We provide the prove of the convergence for the above updating rules and Theorem 2 in Section 3.2.2.

Proofs of Convergence and Theorems
We construct an auxiliary function G(x, x ) to help prove Theorem 1 & 2, which satisfies the conditions G(x, x ) F(x) and G(x, x) = F(x), and guarantees F(x) to be non-increasing under the following update: We illustrate our proofs for W ab and the proofs for H ab can be derived in a similar fashion. Let F ab (W ) denote part of F(W ) concerning W ab and so as G ab (W , W t ).

Convergence of Euclidean updating rules and Theorem 1
As for updating rules in Eq.(15) & (16), we have is an auxiliary function for F ab (W ).
Comparing Eq. (26) and Eq. (27), in order to prove Rewrite F(W ) as an addition of two functions (omitting the penalties for H since it is irrelevant with W ) Then Eq. (28) becomes the only remaining work is to prove Recall that the function V(W ) is a summation of additions of two natural exponential functions, and the exponents are quadratic functions of W . Calculation of their derivatives is intricate and causes difficulties to prove the above inequality. Here we give a proof by applying approximations to the natural exponential functions.
For the general natural exponential function f (x) = e x with support [0, +∞), we can use a piecewise function to approximate it where each sub-function is a linear function of x. The relationship between f (x) and the piecewise function is demonstrated in Fig.  3. Similarly, we can approximate g(x) = e −x with a piecewise function whose n th sub-function is (e −n − e −(n−1) )x − (n − 1)e −n + ne −(n−1) , x ∈ [n − 1, n).
(34) If the support [0, +∞) is divided into infinite intervals, the piecewise functions will be infinitely close to f (x) and g(x). (35) Its first and second order derivatives w.r.t W ab are And the positive part of the first derivative is (38) Then we have (39) Thus Eq. (28) holds. Now we can prove the convergence of Theorem 1: Proof of Theorem 1. According to Eq.(25) and Eq.(26), we get  (41) This is an auxiliary function for F ab (W ).
F ab (W ), we utilize the Jensen Inequality to obtain which holds for all non-negative α k that sum to unity. Setting This is the only different part between F ab (W ) and Now we can prove the convergence of Theorem 2: Proof of Theorem 2. According to Eq.(25) and Eq.(41), we get Notice that W also exists in the penalty regularisation terms, and it is difficult to calculate their corresponding derivatives with respect to W . For simplicity and efficiency, we substitute W in penalties with W t so that they become irrelevant to W . This approximation of the local optima of G ab (W , W t ) is numerically proved feasible through our experiments.
As G ab (W , W t ) is an auxiliary function, F ab (W ) is nonincreasing under this update rule.
The complete algorithms of RPR-NMF are demonstrated in Algorithm 1 and Algorithm 2. for a = 1 to N do 6: Calculate C + row (W ab ) and C − row (W ab ) as Eq. (13) 7: Calculate C + col (H ab ) and C − col (H ab ) as Eq. (13) 12:

EXPERIMENTS
In this section, we conducted experiments on both synthetic datasets and real datasets to demonstrate the performance of our proposed algorithm RPR-NMF. The baseline algorithms we chose for comparison are: NMF [10], GNMF [20], and LCNMF [21]. Each of the algorithms have two versions: one with Euclidean measure (denoted by suffix " euc" in figures and tables), the other with Divergence measure (denoted by suffix " div"). As we stated in the previous section that GNMF and LCNMF are both designed as dimensionality reduction methods, the weight matrix for GNMF only performs as constraints imposed on the right factorised matrix and the label matrix for LCNMF also only affects the right factorised matrix. Thus in our experiments on synthetic datasets as well as datasets for image clustering, RPR-NMF only impose RPR constraints on the right factorised matrix. However, when it comes down to recommender systems, we impose constraints on both factorised matrices which are considered as users' and items' features.
The metrics we used to evaluate the performance of algorithms are: 1. Mean Squared Loss (MSL) for algorithms using Euclidean measure and Mean Divergence (MD) for algorithms using Divergence measure, which evaluate the approximation performance and are defined as Algorithm 2 Relative Pairwise Relationship constrained Nonnegative Matrix Factorisation (RPR-NMF) using Divergence measure Input: V ∈ R N ×M 0 (factorising matrix), K ∈ N (latent dimension), L W ∈ N l W ×3 (pairwise relationship constraints on W ), L H ∈ N l H ×3 (pairwise relationship constraints on H), λ W ∈ R + , λ H ∈ R + (penalty coefficients) Output: W ∈ R N ×K 0 (the left factorised matrix), H ∈ R K×M 0 (the right factorised matrix) 1 where l * W and l * H are the numbers of satisfied constraints. 3. Clustering Accuracy (ACC), which we calculate by constructing a cost matrix, solving the matrix by Munkres Assign Algorithm [25], and mapping one clustering to the other; and Normalised Mutual Information (NMI) [26], which is to calculate the entropy information correlation of clusterings. They are used to evaluate the performance of the clustering results in image clustering experiments. 4. Root Mean Squared Error (RMSE) [27], which evaluates the difference between recovering ratings and original ratings (the latter are removed before factorisation in cross validation); and F1 score [28], which tells how much proportion of correct recommendations the recovered rating matrix gives. They are used to evaluate the performance on recommendation systems and defined as where M is a binary matrix in which ones denote the chosen entries and zeros denote not chosen entries in a cross validation experiment. As for each user, we use its average rating as a threshold. Ratings greater than the threshold suggest the corresponding items are recommended. Thus we obtain two recommendation binary vectors before and after factorisation, which can be used to calculate the value of Precision and Recall. We first introduce two algorithms used to convert the additional information. Then we validate the effectiveness of RPR-NMF and analyse the computing complexity through two synthetic experiments, followed by experiments on real datasets for applications of image clustering and recommender systems. The statistics of the datasets we used are presented in TABLE 1.
All the code of RPR-NMF as well as preprocessed datasets can be downloaded from: https://github.com/shawn-jiang/RPRNMF.

Additional Information Conversion
Since GNMF, LCNMF and RPR-NMF all need additional information besides the factorising matrix itself, it is necessary to make it fair for all the three algorithms to equally access available additional information. GNMF needs a weight matrix which describes how close the data vectors should be, LCNMF requires which data vectors should have the same label, while RPR-NMF demands a list of RPR constraints among factorised vectors. These three additional information can be converted from each to the others. However, their intensities are different: labels of LCNMF are the strongest, followed by the similarities of GNMF, and the RPR constraints of RPR-NMF are the weakest. It will cause loss of information if we convert stronger constraints to weaker ones, but not if it is done the other way around. Thus here we introduce two algorithms which convert the RPR constraints for RPR-NMF to the weight matrix for GNMF (Algorithm 3) and the label matrix for LCNMF (Algorithm 4) respectively.

Effectiveness Validation & Complexity Analysis
We conducted two experiments on synthetic datasets to verify the effectiveness and study the computing complexity of RPR-NMF comparing to the baseline algorithms.
It is worth noting that, the RPR constraints in the experiments are randomly generated but there are chains among them. A chain Put the two nodes into the network (if they have not been there) with an directed edge from node 1 to node 2 7: end for 8: Recursively calculate the depth of each node (nodes with no outgoing edges are set depth 1, their direct parents are set depth 2, and so on). 9: t ← (maxs − mins)/(max(depth) − 1) 10: for each node n do 11: i, j ← two points of n of constraints is like dis(a, b) < dis(b, c) < dis(c, d), which is a 3-chain of constraints. As we are going to show in the results, RPR-NMF can satisfy chain constraints while others are incapable.
In the first experiment, we explored the influence of the number of RPR constraints with respect to MSL/MD, CSR and processing time. First, we randomly generated a 100 × 20 matrix W 0 and a 20 × 100 matrix H 0 , and got the factorising matrix V by their product. Then we generate 1 to 10 groups of 5-chain constraints (i.e. 5 to 50 constraints) from H 0 . K is set 20 and λ H is set 1. For each experiment with different number of constraints, we repeated it 10 times and present the average results in Fig. 4.
As the number of constraints increases, (a) NMF and RPR-NMF algorithms obtain similar and the lowest errors (average 1.963×10 −4 MSL and 2.197×10 −4 MD for NMF, 1.984×10 −4 MSL and 2.195 × 10 −4 MD for RPR-NMF), while the errors of GNMF and LCNMF are increasing and much higher than those of NMF and RPR-NMF; (b) RPR-NMF using Divergence measure achieves the highest CSRs (average 99.11%) followed by its Euclidean version (88.72%), and NMF, GNMF, LCNMF obtain 80.16%, 64.10%, 16.14% average CSR respectively; (c) the processing time of RPR-NMF and the Euclidean version of NMF is increasing while others' is nearly unchanged or slightly decreasing, and GNMF are the fastest while LCNMF are the slowest when comparing algorithms solely with Euclidean measure or Divergence measure.
The second experiment focused on how the size of factorising matrix affects the performance of different algorithms. The size of the factorising matrix V varies from 20 to 200 with step 20. Each V is also obtained by the product of randomly generated W 0 and H 0 . For each matrix size, the latent dimension K is set N/5, and the number of constraints equals to K. The constraints are also of 5-chains. λ H is set 1.The average results of 10 repeated experiments are presented in Fig. 5.
As the size of factorising matrix increases, (a) NMF and RPR-NMF algorithms obtain similar and the lowest errors again with slight rises, while the errors of GNMF and LCNMF are decreasing but still much higher than those of NMF and RPR-NMF; (b) the performance on CSR is similar to that in the first experiment, where RPR-NMF using Divergence measure achieves the highest score; (c) the processing time of all algorithms is increasing, and GNMF are still the fastest while LCNMF are still the slowest respectively for Euclidean measure and Divergence measure.
The results on CSR also suggest that the weight matrix of GNMF and the labels of LCNMF are not helpful on satisfying such RPR constraints. When there are only independent constraints, GNMF and LCNMF can satisfy them by setting proper weights and labels (e.g. for a constraint dis(a, b) < dis(b, c), GNMF sets weight 1 for (a, b) and 0 for (b, c), and LCNMF sets (a, b) the same label). However, when the constraints are not independent (such as in chains), GNMF cannot simply set the weights 1s and 0s, and LCNMF can only guarantee satisfying one constraint of a group of chain constraints (e.g. for a 3-chain of constraints dis(a, b) < dis(b, c) < dis(c, d), the weight of (b, c) in GNMF has to be greater than 0 and less than 1, while in LCNMF, setting (a, b) the same label can satisfy the first constraint, but the only way to satisfy the second constraint is to set (b, c) the same label, which will contradict with the first constraint). Thus, GNMF and LCNMF cannot guarantee most of the constraints are satisfied after factorisation (recall the example in Fig. 2).
The above two synthetic experiments both demonstrate that RPR-NMF algorithms have the advantages of getting accurate factorisation, satisfying relative constraints and being applied to large scale datasets comparing with other algorithms.

Parameter Selection
The parameters introduced by RPR-NMF are the penalty coefficients. We conducted experiments with varying values of these parameters under the following settings: N = 100, M = 100, K = 20, l W = l H = N/10. λ W and λ H vary from 0.4 to 4 with step 0.4 and from 20 to 100 with step 20. The average results of 10 repetitions are presented in Fig. 6.
As showed in the figure, the coefficients have nearly no influence on the CSR for both algorithms. However, the MSL/MD of RPR-NMF using Euclidean measure increases when the parameters become larger, while the Divergence version is stable all the time. This suggests that RPR-NMF is not sensitive to its parameters which can be set properly without much effort.

Performance Analysis in Image Clustering
We validate the performance of algorithms in image clustering on two public datasets with ground truth: AT&T ORL and CMU PIE.

AT&T ORL Dataset
The AT&T ORL database consists of 400 images for 40 classes with 10 different facial images in each class [29]. The images were taken at different times, lighting and facial expressions. The faces are in an upright position in frontal view, with a slight leftright rotation. Each image is preprocessed into a 32 × 32 matrix with 256 grey levels [21]. Thus the size of factorising matrix is 1, 024 × 400.
We adopt the following steps in this experiment: i. Randomly choose K classes and mix up images from these classes to form the factorising matrix; ii. In the K classes, randomly select 2 images in each class and set them more similar to images chosen in other classes, which forms the list of RPR constraints for RPR-NMF; transform the constraints into a weight matrix and a label matrix; the latent dimension is set as the number of chosen classes K; iii. Run algorithms to obtain the right factorised matrix H; utilise K-means method on H to get the clustering results; iv. Calculate the ACC and NMI for each algorithm. The number of class K varies from 2 to 10 with step 2. The penalty coefficients for RPR-NMF using Euclidean measure are   [30]. We followed the pre-processed PIE dataset used in [20] which contains 2, 856 images for 68 different people with 42 images for each person. The images are processed into 32 × 32 matrices denoting the grey level of pixels. Thus the size of factorising matrix is 1, 024 × 2, 856.
We adopt similar experimental steps as we did for AT&T ORL dataset with a few changes on the extraction of RPR constraints. For this experiment, we randomly select 4 images instead of 2 in each cluster. Moreover, considering that the similarity can exist not only among intra-class images, but also inter-class images (e.g. the image of a dog is more similar to that of another dog, rather than a cat), we extract constraints in both ways. The number of class K varies from 10 to 60 with step 10. The penalty coefficients for RPR-NMF using Euclidean measure are set 20 and for its Divergence version are set 2. The results are presented in TABLE  3.
According to the table, we can see that GNMF using Euclidean measure and NMF using Divergence measure achieve the best average approximation (135.5 and 1.242). As for the other three evaluation metrics, RPR-NMF using Euclidean measure achieves the best average CSR (87.03%) and NMI (73.76%) while RPR-NMF using Divergence measure achieves the best average ACC (66.02%). Notice that the CSR of LCNMF are all zeros because inter-class and intra-class constraints lead to cyclic chain constraints. Among the algorithms using Euclidean measure, RPR-NMF outperforms NMF, GNMF and LCNMF by 6.40%, 11.94%, 14.67% on ACC and by 2.66%, 5.37%, 16.02% on NMI respectively. For the algorithms using Divergence measure, RPR-NMF outperforms NMF, GNMF and LCNMF by 7.59%, 24.85%, 18.74% on ACC. However, its performance on NMI is a bit lower than NMF, while it outperforms GNMF and LCNMF by 16.25%, 14.23% on NMI. The average improvement of RPR-NMF algorithms is 14.03% on ACC and 9.02% on NMI.

Performance Analysis in Recommender Systems
As for the performance analysis in recommender systems, we compare our algorithm RPR-NMF with NMF, GNMF, and LC-NMF on Movielens 1M dataset. For the reason that the rating matrices in recommender systems have missing values which are usually denoted by zeros, all the algorithms have to be modified with a MASK matrix for incomplete factorising matrix as proposed in [31]. To our best effort, we implemented the modified version for all algorithms except for the GNMF using Divergence measure. In the Appendix of [20], the authors only mentioned how to deal with incomplete factorising matrix for GNMF using Euclidean measure. Thus we did not compare GNMF using Divergence measure in this part.
Specifically, the pairwise relationship constraints in this part are extracted from meta information and are imposed on both factorised matrices: we utilised users' gender, age and occupation as well as movies' genre to obtain RPRs for the Movielens dataset.

Movielens 1M Dataset
The Movielens 1M dataset is a well-known stable baseline dataset in recommender systems [32]. It has 1, 000, 209 ratings (1 to 5) from 6, 040 users on 3, 883 movies. Indeed some movies have no ratings, thus after removing these movies, we derived a preprocessed 6, 040 × 3, 706 rating matrix. Three groups of cross validation experiments are conducted on this dataset: (1) 300 constraints, K = 20; (2) 500 constraints and K = 50; (3) 1000 constraints and K = 100. The penalty coefficients for RPR-NMF using Euclidean measure are set 200, 10 and 1 respectively, and the coefficients for the Divergence version of RPR-NMF are set 0.1, 0.01, 0.001 respectively. For each group, we conducted 5 cross validation experiments and the average results are showed in TABLE 4. As presented in the table, NMF achieves the lowest MSL and RPR-NMF achieves the lowest MD; both of the algorithms obtain close approximation errors while the other two methods, GNMF and LCNMF, have a higher rate of errors. RPR-NMF also achieves the highest CSR on both its versions. As for the recommendation evaluation, GNMF using Euclidean measure achieves the lowest RMSE while RPR-NMF using Euclidean measure achieves the highest F1 score. It is worth noting that, in this experiment, the RPRs are randomly selected through meta information, and some of the constraints may not represent the true ratings' pattern. Thus the RPR-NMF algorithms seem not greatly outperform existing methods. However, their overall performance is still better compared to other alternatives.

CONCLUSIONS
In this paper, we proposed a novel matrix factorisation algorithm called RPR-NMF, to effectively utilise the relative pairwise relationship among rows or columns of factorised matrices. Both of the Euclidean and Divergence measures are used in the objective function. RPR-NMF imposes penalties for each relative pairwise relationship constraint in a form of addition of natural exponential functions for Euclidean measure and hinge loss for Divergence measure. Complete and sufficient proofs of convergence are also provided to ensure that RPR-NMF conforms to the "multiplicative update rules". Numerical analysis shows that the proposed algorithm achieves superior performance on both the overall loss and the accuracy of satisfied constraints compared with the other algorithms. Experiments on synthetic and real datasets for image clustering and recommender systems demonstrate the effectiveness of RPR-NMF algorithms which outperform baseline methods on several evaluation criteria.