{'fraction_non_alphanumeric': 0.03484995159728945, 'fraction_numerical': 0.009021441370929539, 'mean_word_length': 4.711328079049523, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Scalable and interpretable product recommendations via overlapping co-clustering

We consider the problem of generating interpretable recommendations by identifying overlapping co-clusters of clients and products, based only on positive or implicit feedback. Our approach is applicable on very large datasets because it exhibits almost linear complexity in the input examples and the number of co-clusters. We show, both on real industrial data and on publicly available datasets, that the recommendation accuracy of our algorithm is competitive to that of state-of-art matrix factorization techniques. In addition, our technique has the advantage of offering recommendations that are textually and visually interpretable. Finally, we examine how to implement our technique efficiently on Graphical Processing Units (GPUs).


I. INTRODUCTION
Research on recommender systems is slowly transitioning from its historical focus on prediction accuracy towards a more balanced approach between accuracy and interpretability. Interpretability adds trust, confidence [11], and persuasiveness [50] to the recommendations, and in certain scenarios helps highlight tradeoffs and compromises [25] and even incorporates aspects of social endorsement [46]. Therefore, interpretable recommendations are of interest in many application areas. Here, our goal is to provide interpretable client-product recommendations in a business-to-business (B2B) scenario. In this scenario, the recipient of the recommendation is a salesperson responsible for the client, and not the client itself. A salesperson is responsible for tens or hundreds of clients and to decide whether to pursue a sales opportunity (i.e., recommendation), he or she relies on evaluating the reasoning provided by a generated recommendation.
In a business-to-business, and many other recommender systems, often only positive ratings are present: the products that the clients have already purchased. Negative ratings are unavailable, because absence of a purchase does not necessarily reflect a lack of interest in the item. The problem of generating recommendations based on positive ratings of users only is known as One-Class Collaborative Filtering (OCCF). The lack of negative examples makes the problem challenging, as one has to learn the customer's preferences from what they like, without information on what they dislike. In fact, it can be shown that in an online setup, learning only from what users like requires more examples for making "good" recommendations [15] than learning from both what users like and dislike. This problem also occurs in other important collaborative filtering settings, in particular when only implicit ratings are available. Examples of implicit ratings are the browsing history of a visitor at an e-commerce website or views of videos on a video-sharing platform.
State-of-the-art approches to generate recommendations from positive ratings only are often based on standard matrix factorization. However, they offer low interpretability because "latent factors obtained with mathematical methods applied to the user-item matrix can be hardly interpreted by humans" [35]. Similar observations have been made in several works [3], [36], [51].
Our approach to providing interpretable recommendations from positive examples is based on the detection of co-clusters between users (clients) and items (products) 1 . Co-clusters are groups of both users and items with similar patterns. Users can have several interests, and items might satisfy several needs, so users and items may belong to several co-clusters. So, coclusters may be overlapping. Most importantly, discovery of the overlapping user-item co-clusters offers an interpretable model: Identification of sets of users that are interested in or have bought a set of items not only allows us to infer latent underlying patterns but can also lead to better, textually and visually interpretable recommendations.
Co-clustering principles have been considered before in the context of recommender systems [45], [19], but past work has focused on non-overlapping co-clusters, that is, every item/user pair belongs to only one co-cluster. In contrast, in this work we use overlapping co-clusters to generate interpretable recommendations. Allowing co-clusters to overlap is not only crucial for the recommendation performance, but also enables the discovery of interesting and useful buying patterns. In Figure 1, we provide an example of overlapping co-clusters. A dark square describes a product bought by the user in the past. One can visually identify three potential recommendations indicated by white squares inside the coclusters. The approach we propose identifies overlapping useritem co-clusters and generates recommendations based on a generative model. The algorithm is scalable, on par with state-of-art OCCF methodologies in terms of recommendation accuracy, and at the same time provides interpretable recommendations. The algorithm's computational complexity is very low, specifically it is essentially linear in the problem size. We have also developed highly efficient GPU implementations of the algorithm.
Outline: The paper first reviews related work in Section II. Next, we present our Overlapping co-CLuster Recommendation algorithm, henceforth referred to as "OCuLaR". We start by introducing a generative model, based on which we design a computationally efficient algorithm to produce recommendations and to identify the co-clusters. In Section IV-E, we explain the relation to traditional matrix factorization approaches and in Section VI we discuss how to port our approach onto GPUs. Section VII compares OCuLaR to other state-of-theart recommendation algorithms on a real-world client-product dataset from our institution and publicly available datasets. Our algorithm performs as well as or better than state-ofthe-art matrix factorization techniques, with the added benefit of interpretability. We conclude with a deployment of our algorithm in a B2B recommender system at our organization.

II. RELATED WORK
Collaborative Filtering (CF): Early approaches to recommender systems performed either user-based or item-based collaborative filtering. User-based techniques infer preferences of a given user based on the preferences of similar or likeminded users by, e.g., by recommending products that the nearest neighbors of a user have bought in the past. Similarly, item-based techniques exploit item-to-item similarities to generate recommendations. Item-and user-based techniques yield a reasoning of the sort "similar users have also bought", but are often outperformed by latent factor models, such as matrix factorization approaches [20]. Matrix factorization techniques in their traditional form predict ratings or preferences well, but the latent features make it difficult to explain a recommendation [3], [36], [35]. Recently, [51] investigated a method for explainable factorization by extracting explicit factors (sentiment, keywords, etc) from user reviews; however such an approach is applicable only in the presence of additional textual information for each training example. Such data, however, are rarely available in a B2B setting. Generic methodologies for communicating explanations to the user have been studied in [43], [11], [12].
Conceptually, the algorithm proposed in this work can be regarded as a combination of item-based and user-based approaches. Specifically, we discover similarities that span across both the user and the item space by identifying user-item coclusters. By identifying user-item co-clusters via matrix factorization techniques, we obtain an estimate of the probability that an unknown example is positive, which automatically results in ranked recommendations.
OCCF models are used to predict preferences from implicit feedback. This is a typical recommendation scenario when browsing, buying, or viewing history is available, or in general in setups where the user does not provide any explicit negative feedback. Recommendation techniques based on positive ratings try to learn by considering the implicit feedback either as absolute preferences [29], [17], or as relative preferences [34], [30]. We compare with techniques from both categories in Section VII, and state the formal relation of our approach to matrix factorization techniques in Section IV-E.
Co-clustering: The majority of the literature on co-clustering considers the detection of non-overlapping co-clusters [10], [31]. Notable exceptions are the approach in [4], which introduces a generative model for overlapping co-clusters along with a generic alternating minimization algorithm for fitting the corresponding model to given data, and the approach in [40], which considers the problem of simultaneously clustering documents and terms. That latter approach is based on a generative model and estimation of the corresponding model parameters [40]. In both [4], [40], the focus is on discovering co-clusters, whereas here we are interested in explaining the recommendations produced by our algorithm with co-clusters. Although co-clustering approaches have been used before in the collaborative filtering setting, the majority of those papers [13], [19], [45] is restricted to non-overlapping co-clusters. An exception is [48], which explores the multiclass co-clustering problem in the context of CF.
Community detection: Related to co-clustering is community detection. To see this, observe that the positive examples are the edges in a bipartite graph of users and items (see Figure  1). Then, the co-clusters of users and items correspond to user-item communities in the bipartite graph. One of the bestknown community detection algorithms is based on the notion of modularity [26]. Specifically, the modularity algorithm by Girvan & Newman [14] is one of the most widely-used community detection algorithms and is used in many software packages (gephy, Mathematica, etc). It has the advantage that it can automatically discover the number of communities; however it does not support discovery of overlapping communities. Other work on detecting non-overlapping communities includes [18], [9].
Recently, interest in identification of overlapping communities has grown [28], [2], [22], [49], [32]. Most related to our work is the BIGCLAM algorithm proposed by Yang and Leskovec [49]; the main differences are that in our approach we consider a particular bi-partite graph and use regularization, which turns out to be crucial for recommendation performance. While such approaches share similar concepts, off-the-shelf community detection methodologies are not directly applicable to the one class collaborative filtering problem, as they yield an assignment of users/items to communities, but not a ranked list of recommendations. In fact, they may even fail to reveal the correct community structure. As an illustrative example, in Figure 2 we provide the output of both Modularity and BIGCLAM for the introductory example. Both fail to reveal the correct co-clustering structure, and by recovering incorrect 'community' boundaries they would have identified only one (1) of the three (3) candidate recommendations.

III. B2C VS B2C RECOMMENDATIONS
Our main application of interest, albeit by far not the only one, in the problem of generating interpretable recommendation from one-class examples are B2B recommender systems. A B2B recommender system-in contrast to a Business-to-Consumer (B2C) recommender system (referring to recommender systems employed at sites such as Amazon.com or the Genius recommendations of Apple within their iTunes platform)-is typically not open to the public, but is deployed within an enterprise. Because only a handful of literature is available on B2B recommender systems [27], [52], [1], we briefly point out some differences between a B2C and B2B deployment. This discussion will serve the purpose of highlighting the need for higher interpretability in a B2B recommendation environment.
1) The system users in a B2C system are the consumers, whereas in B2B they are the selling teams or marketing teams of an enterprise. Therefore the recommendations have a different audience. 2) The user entities ("clients") for which the recommendations are created are typically individuals in a B2C setting. In B2B, they are enterprises or companies. This has certain implications on data availability. For example, in B2C we have knowledge about past buying patterns (PBP) or implied interest from items browsed by the users. In B2B, we know the past buying patterns and it is also possible to obtain other type of information about the clients-companies in the form of news events or financial reports published about that company. This eventually can help build a more holistic view of the examined bipartite graph of clients and products.
3) The value of the recommendation in a B2C setting is typically small, ranging from a few dollars (a song or a movie) to a few hundreds of dollars (an electronic gadget). In contrast, in a B2B setting, the recommendation may refer to a large hardware installation or complex and custom software solution and the recommendation may be valued at thousands of dollars, or even millions. It is also interesting to note that in B2C the value of a product/recommendation is (mostly) fixed, or it is publicly available, whereas in B2B it varies based on the complexity of the installation. 4) Finally, the reasoning of the recommendation is short and limited in a B2C scenario, but it has to be detailed and offer sufficient supportive evidence in B2B. This is partially because of the associated value of the recommendation. If a salesperson is to pursue a multithousand dollar sales opportunity, the methodology has to present sufficient reasoning for this action, which is subsequently further evaluated by the responsible seller or marketing team.
In this work we will only focus on the last part, the reasoning of the recommendation. Even though components such as estimating the appropriate value of a B2B recommendation are also non-trivial, they will not be further examined in this work.
In the following section we begin to describe our approach for generating interpretable recommendations based on the client-product purchase graph.

IV. OCULAR ALGORITHM
Now we present our Overlapping co-CLuster Recommendation algorithm, or for short "OCuLaR". We assume that we are given a matrix R where the rows correspond, e.g., to users or clients and the columns to items or products. If the (u, i)th element of R takes on the value r ui = 1 this indicates that user u has purchased item i in the past or, more generally, that user u is interested in item i. We consider all values r ui that are not positive (r ui = 1) as unknown (r ui = 0), because they indicate that user u might be interested in i or not. Our goal is to identify those items a user u is likely to be interested in. Put differently, we want to find the positives among the unknowns, given only positive examples.
We assume an underlying model whose parameters are factors associated with the users and items. Those factors are learned, such that the fitted model explains well the given positive examples r ui = 1. The specific choice of our model allows us to design an efficient algorithm to learn the factors. Moreover, the factors encode co-cluster membership and affiliation strength, and are interpretable in that sense.

A. Generative model
We start with the generative model underlying our recommendation approach. It formalizes the following intuition: There exist clusters, groups, or communities of users that are interested in a subset of the items. As users can have several interests, and items might satisfy several needs, each user and item can belong to several co-clusters consisting of users and items. However, a co-cluster must contain at least one user and one item, and can therefore not consist of users or items alone.
Suppose there are K co-clusters (K can be determined from the data, e.g., by cross-validation, as discussed later). Affiliation of a user u and item i with a co-cluster is modeled by the K-dimensional co-cluster affiliation vectors f u and f i , respectively. The entries of f u , f i are constrained to be nonnegative, and [f u ] c = 0 signifies that user u does not belong to co-cluster c. Positive examples are explained by the co-clusters as follows. If user u and item i both lie in co-cluster c, then this co-cluster generates a positive example with probability Assuming that each co-cluster c = 1, ..., K, generates a positive example independently, it follows that A similar generative model also appears in the community detection literature [49]. So far the model cannot explain a positive example by means other than co-cluster affiliation. Although we could incorporate user (b u ), item (b i ), and overall bias (b) by supposing that the probability of an example being positive is given by we found that fitting the corresponding model does not increase the recommendation performance for the datasets considered in Section VII, and therefore will not discuss it further.

B. Fitting the model parameters
Given a matrix R, we fit the model parameters by finding the most likely factors f u , f i to the matrix R by maximizing the likelihood (recall that we assume positive examples to be generated independently across co-clusters and across items and users in co-clusters): Maximizing the likelihood is equivalent to minimizing the negative log-likelihood: To prevent overfitting, we add an 2 penalty, which results in the following optimization problem: where and λ ≥ 0 is a regularization parameter. As will we discuss in more detail in Section IV-E, this optimization problem can be viewed as a variant of non-negative matrix factorization (NMF), specifically NMF with a certain cost function. A common approach to solve an NMF problem is alternating least squares, which iterates between fixing f u , and minimizing with respect to f i , and fixing f i and minimizing with respect to f u , until convergence. This strategy is known as cyclic block coordinate descent or the non-linear Gauss-Seidel and convex in f u (with f i fixed). Therefore, a solution to the subproblems of minimizing Q with fixed f i and minimizing Q with fixed f u can be found, e.g., via gradient descent or Newton's method. As this optimization problem is non-convex, one cannot in general guarantee convergence to a global minimum; however convergence to a stationary point can be ensured. Specifically, provided that λ > 0, Q is strongly convex in f i (with f u fixed) and in f u (with f i fixed). Thus, the subproblems have unique solutions, and therefore, if we solve each subproblem exactly, convergence to a stationary point is ensured by [ However, as noted in the context of matrix factorization [16], solving the subproblems exactly may slow down convergence. Specifically, when f u , f i are far from a stationary point, it is intuitive that there is little reason to allocate computational resources to solve the subproblems exactly. It is therefore often more efficient to solve the subproblem only approximately in each iteration [6], [16].
For the above reasons, we will only approximately solve each subproblem by using a single step of projected gradient descent with backtracking line search, and iteratively update f i and f u by single projected gradient descent steps until convergence. In Section IV-D, we provide further implementation details. Convergence is declared if Q stops decreasing. This results in a very efficient algorithm that is essentially linear in the number of positive examples {(u, i) : r ui = 1}, and the number of co-clusters K. Our simulations have shown that performing only one gradient descent step significantly speeds up the algorithm.
Choice of K and λ: Recall that the number of co-clusters K and the regularization parameter λ are the two model hyperparameters. K and λ can be determined from the data via cross-validation. Specifically, to determine a suitable pair of (K, λ), we train a model on a subset of the given data for different choices of (K, λ), and select the pair for which the corresponding model performs best on the test set. The model can be tuned for an appropriate metric. In our experiments, we measure the recommendation performance in terms of the recall-at-M items [24], which is a typical performance metric in recommender systems. Because this grid-search approach can be costly if one wishes to perform a fine-grained search of the hyper-parameter space, the GPU implementation of the algorithm, reported in Section VI, can help to dramatically reduce the overall grid-search time.

C. Generating interpretable recommendations
Suppose we want to recommend M items to each user. After having fitted the model parameters, we recommend item i to user u if r ui is among the M largest values P [r ui = 1], where i is over all items that user u did not purchase, i.e., over all i with r ui = 0. The probability P [r ui = 1] is large if the user-item pair (u, i) is in one or more user-item coclusters. Thus, along with a recommendation, we can output the corresponding user-item co-clusters that cause P [r ui = 1] or, equivalently, In the B2B setting that we consider it is also important to explicitly mention who are the clients that have purchased a similar bundle of products. Contrary to the B2C setting, where, for privacy reasons, one only mentions what similar clients purchase, in a B2B setting this is not a concern. The salesperson, who is the recipient of the recommendation, can use this information (explicit names of similar clients) to understand better the types of clients/companies that typically require such a solution. Our approach directly provides this information, because each co-cluster consists of specific clients (users) and products, and does not merely describe an average behavior.  To justify recommending Item 4 to User 6, we can therefore give the following, automatic interpretation to the user of the recommender system: Item 4 is recommended to Client 6 with confidence 0.83 because: • Client 6 has purchased Items 1-3. Clients with similar purchase history (e.g., Clients 4-5) also bought Item 4. • Moreover, Client 6 has purchased Items 5-9. Clients with similar purchase history (e.g., Clients 7-9) also bought Item 4.
Naturally, once the co-clusters have been discovered, additional information (derived from the co-clusters) can be attached in the rationale presented. In the experimental Section VII, we discuss the recommendation rationale for an industrial deployment of OCuLaR (see Figure 10).

D. Implementation and complexity
Here we examine in more detail the projected gradient descent approach we use to solve the subproblems and the complexity of the overall optimization algorithm. It is sufficient to discuss minimization of Q with respect to f i , as minimization with respect to f u is equivalent. The following approach for minimizing Q was also used in similar form in [49], [23]. We start by noting that, because of we can minimize Q for each f i individually. The part of Q depending on f i is given by As mentioned, we update the parameter f i by performing a projected gradient descent step. The projected gradient descent algorithm [5,Sec. 2.3] is initialized with a feasible initial factor f 0 i and updates the current solution f k i to f k+1 where (f ) + projects f on its positive part, [(f ) + ] c = max(0, [f ] c ), and the gradient is given by The step size α k is selected using a backtracking line search, also referred to as the Armijo rule along the projection arc [5]. Specifically, α k = β t k , where t k is the smallest positive integer such that where σ, β ∈ (0, 1) are user-set constants. As the computation of both ∇Q(f i ) and Q(f i ) requires u : rui=0 f u , and typically, the number of items for which r ui = 1 is small relative to the total number of items, we precompute u f u before updating all f i , and then compute u : rui=0 f u via This idea is taken from [49], where it was used in the context of community detection. Using the precomputed

E. Relation to Matrix Factorization
Latent factor models are a prevalent approach to recommender systems. Among them, matrix factorization approaches are particularly popular, because of their good performance properties. In this section, we discuss the connection of our approach to standard matrix factorization techniques. Matrix factorization models in their most basic form represent users and items by latent vectors f u , f i in a low-dimensional space R K , and form a rating according to

A common approach to fit the latent factors based on a given set of examples {r
where is a loss function and f i 2 2 , f u 2 2 are regularization terms to prevent overfitting. A common choice for non-binary ratings is the quadratic loss function (r ui , For the one-class collaborative filtering problem, quadratic loss is not directly applicable, as the question remains on how to deal with the unknowns, for which r ui = 0. Performing summation in (7) only over the positive examples (r ui = 1) is not sensible, as it results in a trivial solution of f u = f i , for all u and all i to (7). A different approach is to treat the unknowns as negative ratings. However, this might bias recommendations as some of the unknown examples are actually positive ratings. To resolve this issue, it has been proposed in [29] to give different weights to the error terms corresponding to positive and unknown ratings in the objective function, specifically to use the cost function where b < 1 is a weight assigned to unknown ratings. A common approach to solve the corresponding optimization problem is weighted alternating least squares (wALS) [29]. Our approach is equivalent to choosing the loss function This results in a large penalty if f i , f u is small for a positive example r ui = 1, and a moderate penalty if f i , f u is large for an unknown example r ui = 0. Our approach is therefore similar in spirit to that of giving different weights to positive (r ui = 1) and unknown (r ui = 0) examples.
Interpretability: Matrix factorization approaches, such as the wALS algorithm (with loss function (8)), yield good empirical performance, as we will see in the experiments. Their main disadvantage for our scenario is that the latent space is typically not easy to interpret. This drawback of MF techniques is attested in various studies [21], [3], [36], [44]. This statement also applies to standard non-negative matrix factorization (NMF) techniques, where the factors are constraint to be non-negative. Our approach also uses factorization principles, however, we confine the factors to explicitly model user and item participation, which is key for not compromising the model interpretability.

V. RELATIVE OCULAR ALGORITHM (R-OCULAR)
As we saw in the related work section, OCCF problems can be viewed as learning of either absolute or relative preferences of users. The OCuLaR algorithm belongs to the first category because it tries to learn the absolute ratings of the users. Specifically, as we saw, OCuLaR can be viewed as a nonnegative matrix factorization approach with a particular loss function (i.e., (9)). The loss function assigns a large penalty to positive examples that are not well explained by the model (i.e., f i , f u small for r ui = 1), and only a moderate penalty to unknown examples that are not well explained by the model (i.e., f i , f u large for r ui = 0).
In this section, we examine how OCuLaR could be adapted to treat the positive examples as relative preferences. The notion of relative preferences in OCCF was first explored by Rendle et al. [34], who proposed to predict the relative preferences (u prefers i over j) rather than the absolute rankings. The underlying idea for this case is that there is an associated latent personalized ranking > u for each user u, where i > u j signifies that user u prefers item i over item j. To this end, in a first step, one can construct a "training" item ranking set for each user from the set of positive examples, denoted by S = {(u, i) : r ui = 1}. The underlying assumption for this constuction is that, if r ui = 1, then user u prefers item i item over all items j with unknown rating, i.e., r uj = 0. Specifically, the training item ranking data set D S is defined as: D S = {(u, i, j) : r ui = 1 and r uj = 0}.
Note that negative examples are accounted for implicitly because if (u, i, j) ∈ D S , then j is not prefered over i. Triplets not in D S correspond to triplets for which no direct preference information is available; it is those triplets for which we want to learn the preferences. As shown in [34], assuming that the users act independently of each other and that > u is independent across u, the model likelihood can be maximized by maximizing To adapt the notion of relative preferences for OCuLaR, suppose that the probability of a user preferring item i over item j is given by To see the formal relation to the model in Section IV-A, simply note that, with this choice, with P [r ui = 1] as defined in (1). Maximizing (10) is equivalent to minimizing the logarithm of (10), i.e., minimizing where we defined w u = |{i : rui=0}| |{i : rui=1}| . Note that the RHS above differs only in the factors w u from the negative log-likelihood (2) in Section IV-B. The factors w u assign a large weight to positive examples of users that have few positive examples associated with them.
We denote the algorithm obtained by substituting the loglikelihood − log L in Section IV-B with the RHS above as the relative OCuLaR (R-OCuLaR) algorithm. Its implementation is essentially equivalent to that of the (original) OCuLaR algorithm; in fact, it has exactly the same complexity.

VI. USING MASSIVELY PARALLEL PROCESSORS
The OCuLaR algorithm can be dramatically accelerated by exploiting its inherent parallelism. In fact, one can easily map and implement the costly training phase onto a graphics processing unit (GPU). This results in two benefits: • Leverage the reduction in training time to train the model and produce recommendations in significantly reduced time. For the real industrial dataset used in our experiments, both training of the model and generation of recommendations can be completed in mere seconds when using a GPU. • Perform a more fine-grained grid search over the algorithm's hyper-parameters to obtain a higher precision/recall. OCuLaR requires the learning of K (number of co-clusters) and λ by using a cross-validated grid search. This phase can be costly, so in practice one can only search across a certain range of values and keep the pair that results in the best performance under the desirable accuracy metric. Achieving very fast learning using GPUs essentially allows the methodology to examine a broader range of hyper-parameter value pairs that help improve the overall accuracy. In order to achieve the potential acceleration offered by GPU devices, communication of large amounts of data between the host and the device must be avoided, especially during iterative computations. Furthermore, arithmetic computations must be carefully designed to optimally make use of the underlying computational architecture and memory hierarchy. Below, we describe in more detail how OCuLaR can be mapped onto an efficient GPU implementation.

A. GPU implementation
We begin by describing how data moves between the host and the GPU during the training of the OCuLaR algorithm. Firstly, the entire training data (in a sparse format) is copied from the host memory into the GPU main memory along with a set of initial values for the f i and f u vectors. The iterative training algorithm described in Section IV-D is then performed by launching a sequence of kernel functions. Throughout the training, all data remains on the GPU and communication back and forth between the host and the device is limited to a small amount of control logic. Once a predetermined set of iterations has been executed, the learned values of the f i and f u are copied back from the GPU to the host memory and the training is completed.
The arithmetic computations that comprise the training are split into a number of separate kernel functions. As an illustration, we provide a detailed description of how the kernel that computes the gradient vectors for all items is structured to efficiently use the underlying hardware. We begin by noting that equation (6) can be expressed: where C is a constant independent of the item index. As in the reference CPU-only implementation, an initial computation is performed to calculate the sum of all f u vectors and the gradient vectors are initialized in GPU memory as C + 2λf i . A kernel function is then called which launches a thread block for every positive rating in the training data. Since this number is typically very large, this mapping is well suited to the massively parallel GPU architecture in which a number of streaming multiprocessors are capable of executing multiple thread blocks concurrently.
We will now describe the computation performed by a single thread block in more detail, making reference to Figure 4. The thread block first fetches the corresponding f i vector and f u vector from GPU main memory and computes the inner product between the two. Following the approach of [37, p. 76], an individual thread within the block handles computation of only part of the inner product. The memory access patterns are carefully arranged so that threads that are executed concurrently (as a "warp") access contiguous regions of memory, allowing reads to be effectively coalesced. The partial results are then stored in high-speed shared memory and a reduction is performed to obtain the final value of the inner product (P ). A single thread within the block then computes the scalar α. Finally, the thread block multiplies this scalar by the f u vector and updates the corresponding item gradient in GPU main memory using atomic operations. In this manner, the sum in equation (11) is computed entirely asynchronously for all items. Once all the thread blocks have finished execution, the correct value of the gradient vectors exist in GPU main memory and the training can proceed. The other essential computations, such as the evaluation of the likelihood function are mapped to the GPU in a similar fashion. Memory: The memory footprint of the GPU-based OCuLaR implementation scales as: O (max (|{(u, i) : r ui = 1}|, n u K, n i K)) , where n u and n i are the number of users and the number of items respectively. This property allows for training on very large datasets despite the relatively limited memory capacity of modern GPU devices. For example, around 2.7GB of GPU memory is required to train on the Netflix dataset (assuming K = 200) and thus the problem easily fits within the main memory of an inexpensive GPU (typically up to 12GB). In contrast, a previous attempt to implement an alternatingleast-square based matrix factorization approach on GPUs determined that the memory requirements for the same dataset exceeded 12GB (for the equivalent of K=100) [42].

VII. EXPERIMENTS
Here, we compare the prediction accuracy of OCuLaR to that of other interpretable and non-interpretable one-class CF algorithms. We compare with interpretable user-based and item-based collaborative filtering approaches, and noninterpretable state-of-art one-class recommendation algorithms based on matrix factorization. Both OCuLaR and R-OCuLaR typically outperform or are on par with the existing oneclass recommendation algorithms in terms of recommendation performance. We also discuss the choice of the input parameters for our approach. We show that OCuLaR exhibits linear scalability, demonstrate the acceleration achieved by the GPU implementation, and conclude with an industrial deployment.

A. Datasets
We use four datasets. First, we consider a real-world dataset from our institution. It consists of the buying history of 80,000 clients with whom our institution interacts and 3,000 products or services offered by our institution. The clients in this case are not individuals but companies, and the recommender system operates in a B2B setting. We call this dataset B2B-DB. We also consider three public datasets. The first was extracted from the CiteULike website, which assists users in creating collections of scientific articles [47]. The dataset consists of 5,551 users and 16,980 articles. Each user has a number of articles in their collection, which are considered as positive examples. Based on the positive examples, the goal is to generate new article recommendations. The second dataset is the Movielens 1 million dataset, which consists of 1 million ratings from 6,000 users on 4,000 movies. The final dataset is the Netflix dataset, consisting of about 100 millions ratings that 480,189 users gave to 17,770 movies.
In both the Movielens and the Netflix dataset, the users provide ratings between 1 and 5 stars. As we consider a one-class collaborative filtering task, we adopt the convention from many previous works (e.g., [30], [41]) to only consider ratings greater than or equal to 3 as positive examples and ignore all other ratings. Therefore, the task now is equivalent to predicting whether the user will give a rating greater than 3 (i.e., is likely to enjoy the movie).
B. Recommendation performance 1) Evaluation Metrics: We measure performance in terms of recall at M items (recall@M ), and mean average precision at M items (MAP@M). In the one-class setting, recall is a more sensible measure than precision, because an example being unknown (r ui = 0) does not mean that user u would not rate item i positively [39]. Given an ordered (by relevance) list of M recommendations for user u, denoted by i 1 , ..., i M , the recall@M items for user u is defined as The overall recall@M is obtained as the average over recall@M (u).
MAP is commonly used in information retrieval for evaluating a ranked or ordered list of items, and is considered a good measure of performance when a short list of the most relevant items is shown to a user. MAP@M items is the mean Note that because Prec(m) ∈ [0, 1], we have that AP@M ≤ 1.

2) Comparison with baselines:
We compare the OCuLaR algorithm with various interpretable and non-interpretable baseline algorithms: • User-based collaborative filtering using a cosine similarity metric (e.g., [38]). Such an algorithm is interpretable because a recommendation can be justified with a reasoning of the type: "item i is recommended because the similar users u 1 , .., u k also bought item i". • Item-based collaborative filtering using cosine similarity (e.g., [8]). This algorithm is also interpretable. It can be accompanied with a recommendation such as: "item i is recommended because user u bought the similar items i 1 , ..., i k ". • Weighted Alternating Least Squares (wALS), a stateof-the-art one-class matrix factorization approach [29]. wALS minimizes the loss function by alternatingly optimizing the user and the item factor via least-squares. It offers good prediction performance but the reasoning is not directly interpretable, as mentioned earlier. • Bayesian personalized ranking (BPR), a state-of-the-art matrix factorization approach that converts the set of positive examples into a set of relative preferences [34] (see also Section V). The recommendations provided by BPR are also not directly interpretable. For BPR we used the python/theano implementation from https: //github.com/bbcrd/theano-bpr.git. All approaches require the setting of hyper-parameters. For each technique we test a number of hyper-parameters and report only the best results. For OCuLaR and R-OCuLaR, we performed a grid search over K and λ (within the range 100 and 200). As we show later, this range does not include the parameters leading to the 'optimal' results, but it is sufficient to showcase a baseline performance without excessive tuning. For the user-and item-based collaborative filtering approaches, we performed a grid search over the number of nearest neighbors. For wALS, we choose the weight in the loss function (8) as b = 0.01, the regularization parameter as λ = 0.01, and performed a grid search over the dimension of the latent vectors. Finally, for BPR we searched over the dimension of the latent factors K, and the regularization parameter λ. We computed the recall@M and MAP@M by splitting the datasets into a training and a test dataset, with a splitting ratio of training/test of 75/25, and averaging over 10 problem instances. We plot the recall@M and MAP@M for the Movielens dataset for varying M in Figure 5. We see that OCuLaR and R-OCuLaR are consistently better or at least as good as the other recommendation techniques. We summarize the results for the other datasets 2 in Table I. Across all datasets the OCuLaR variants are either the best or the second-best performing algorithm (together with wALS). OCuLaR has the advantage of providing interpretable recommendations, an aspect that is compromised when using other OCCF approaches such as wALS or BPR. OCuLaR is also significantly better than the user-and item-based algorithms, its interpretable competitors. This is not surprising because the user-and item-based approaches both consider only similarities in either the user or the item space, whereas OCuLaR can discover more complex structures in the joint user-item space.

C. OCuLaR parameters
We briefly discuss how to set the parameters for OCuLaR and their impact. Recall that OCuLaR expects the number of co-clusters K and the regularization parameter λ. Values for K and λ are chosen such that the recommendation performance is optimized for a particular metric, as determined via crossvalidation. Figure 6 shows the impact of the parameter values on the recommendation performance and on various co-cluster properties on the MovieLens dataset. The top panel shows the recall@50 items. The graph demonstrates that either too little (λ = 0) or too much regularization (λ = 100) can hurt the recommendation accuracy. K may be selected in such a way to ensure that the size of the co-clusters is neither too big nor too small, and also that each user or item does not belong to too many co-clusters.The size of the co-clusters can depend on application-specific criteria, such as, for example, that a co-cluster should contain at least 100 users. For the specific example, a value of K in the range between 100-200 would be adequate to ensure good prediction and avoid excessively big, thus not dense, co-clusters.

D. Scalability
In this section, we investigate the scalability of the OCuLaR algorithm. We consider the full Netflix dataset, consisting of 100, 480, 507 ratings from 480, 189 users on 17, 770 movie titles. As before, we take the ratings ≥ 3 as positive examples, i.e., r ui = 1. As analyzed in Section IV-D, the training time required by the OCuLaR algorithm is essentially linear in the number of positive examples |{(u, i) : r ui = 1}|, and linear in the number of co-clusters K. In Figure 7 we plot the running time per iteration for increasing fractions of the Netflix dataset (i.e., non-zero entries), chosen uniformly from the whole Netflix dataset. We see that the training time is indeed linear in the number of positive examples |{(u, i) : r ui = 1}| and linear in the number of co-clusters K.

E. GPU implementation
Speedup: Here we compare the performance of the CPU implementation with the GPU implementation for the Netflix dataset. The CPU implementation was written in C++ and utilizes the boost library to accelerate certain parts of the computations (such as the evaluation of inner products). It was executed on an Intel Xeon CPU with a clock speed of 2.40GHz. The GPU implementation was written in CUDA C/C++ and executed on a GeForce GTX TITAN X device. In Figure 8 we plot the complete time progression of the training phase, in which OCuLaR minimizes the likelihood of its loss function. The GPU implementation can achieve the same training accuracy over 50 times faster than the CPU implementation. Grid search: The GPU implementation can also help accelerate the grid-search process for learning the algorithm hyper-parameters, namely the regularization parameter λ and the number of co-clusters K. This process can be further accelerated by scaling out across a GPU-enabled cluster. In Figure 9 we show the result of a very fine grid-search over 625 different parameter pairs for the IBM-B2B dataset. The parameter pairs were distributed using Apache Spark across a cluster of 8 machines, each fitted with a NVIDIA Quadro M4000 GPU. Each Spark worker performed training for the parameters assigned to it by calling compiled CUDA code via the Python-C API. The entire grid-search using 8 GPUs took around 8 minutes. On a single CPU this would have required more than 2 days. In Figure 9, the two rectangles show the range of 'optimal' hyper-parameters lies, and the grid-search range used in the CPU-only experiment (Fig. 5). Therefore, the recall reported in the previous experiments could have been improved even further had one used a more exhaustive hyper-parameter grid-search method. This shows the benefit of having a very fast GPU implementation.  Fig. 9. Grid search for hyper-parameters (K, λ) of the algorithm. The heatmap shows the recall@50 products for the IBM-B2B dataset. 'Hot' areas resulted in higher recall. One can benefit from the accelarated execution on a GPU to perform a more fine-grained search of the hyper-parameter space.

VIII. DEPLOYMENT AND USER FEEDBACK
We used our algorithm in a B2B recommender system of our institution. These recommendations are not offered directly to the clients of our institution, but rather to our sales teams. The salesperson responsible for an account examines the recommendations provided and decides whether to act on the recommendation and approach the client. This decision is based on the reasoning provided and on the salesperson's own experience through past interaction with the client.
The textual output of the recommendation, which conveys the corresponding rationale is shown in Figure 10. For reasons of anonymity, the names of the companies/clients belonging to each co-cluster have been omitted.
In the example, we see that the service "Custom Cloud" is recommended for "Client 1" with a confidence of 65.4%. The reasoning explains that Client 1 belongs to three co-clusters that have bought the same service, as well as products similar Fig. 10. Example from industrial deployment (client names are suppressed). Product "Custom Cloud" is recommended to Client 1 because the client is affiliated with three co-clusters. Co-cluster 1 shows an affinity of Client 1 with several airlines, and in Co-cluster 3 with telco companies.
to the ones that Client 1 has already purchased. Based on the co-clusters discovered, the interface also presents a price estimate of the potential business deal, based on historical purchases of the same product by the related clients belonging to the co-clusters discovered.
During the deployment of the OCuLaR algorithm in our organization, we interacted with many sellers that used the platform. Sellers expressed satisfaction about the reasoned aspect of the recommendations. An interesting comment that we received is that the tool could also constitute an educational platform [7] for young sellers, because the detailed reasoning can teach them the currently discovered buying patterns.

IX. CONCLUSION
A large body of work on recommender systems and machine learning has focused primarily on the accuracy of prediction rather than on interpretability. This work explicitly addresses the aspect of interpretability by enabling the detection of overlapping co-clusters that can be easily visualized and transcribed into a textual description. The methodology presented is interpretable, scalable, and does not sacrifice accuracy of prediction. We demonstrated also an efficient GPU implementation which allows the process to be further accelerated by more than 50 times and helps to better explore the hyperparameter space of the algorithm.
Finally, even though the focus of this work was on recommender systems, we feel that the algorithm presented can be used for solving large co-clustering problems in other disciplines as well, including community discovery in social networks [28], or for the analysis of gene expression data [33].