{'fraction_non_alphanumeric': 0.03271715274792407, 'fraction_numerical': 0.013349322352195851, 'mean_word_length': 4.257911015940994, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Algorithmic Acceleration of Parallel ALS for Collaborative Filtering: Speeding up Distributed Big Data Recommendation in Spark

Collaborative filtering algorithms are important building blocks in many practical recommendation systems. For example, many large-scale data processing environments include collaborative filtering models for which the Alternating Least Squares (ALS) algorithm is used to compute latent factor matrix decompositions. In this paper, we propose an approach to accelerate the convergence of parallel ALS-based optimization methods for collaborative filtering using a nonlinear conjugate gradient (NCG) wrapper around the ALS iterations. We also provide a parallel implementation of the accelerated ALS-NCG algorithm in the Apache Spark distributed data processing environment, and an efficient line search technique as part of the ALS-NCG implementation that requires only one pass over the data on distributed datasets. In serial numerical experiments on a linux workstation and parallel numerical experiments on a 16 node cluster with 256 computing cores, we demonstrate that the combined ALS-NCG method requires many fewer iterations and less time than standalone ALS to reach movie rankings with high accuracy on the MovieLens 20M dataset. In parallel, ALS-NCG can achieve an acceleration factor of 4 or greater in clock time when an accurate solution is desired; furthermore, the acceleration factor increases as greater numerical precision is required in the solution. In addition, the NCG acceleration mechanism is efficient in parallel and scales linearly with problem size on synthetic datasets with up to nearly 1 billion ratings. The acceleration mechanism is general and may also be applicable to other optimization methods for collaborative filtering.


I. INTRODUCTION AND BACKGROUND
Recommendation systems are designed to analyze available user data to recommend items such as movies, music, or other goods to consumers, and have become an increasingly important part of most successful online businesses. One strategy for building recommendation systems is known as collaborative filtering, whereby items are recommended to users by collecting preferences or taste information from many users (see, e.g., [1], [2]). Collaborative filtering methods provide the basis for many recommendation systems [3] and have been used by online businesses such as Amazon [4], Netflix [5], and Spotify [6].
An important class of collaborative filtering methods are latent factor models, for which low-rank matrix factorizations are often used and have repeatedly demonstrated better accuracy than other methods such as nearest neighbor models and restricted Boltzmann machines [5], [7]. A low-rank latent factor model associates with each user and with each item a vector of rank n f , for which each component measures the tendency of the user or item towards a certain factor or feature. In the context of movies, a latent feature may represent the style or genre (i.e. drama, comedy, or romance), and the magnitude of a given component of a user's feature vector is proportional to the user's proclivity for that feature (or, for an item's feature vector, to the degree to which that feature is manifested by the item). A low-rank matrix factorization procedure takes as its input the user-item ratings matrix R, in which each entry is a numerical rating of an item by a user and for which typically very few entries are known. The procedure then determines the low-rank user (U) and item (M) matrices, where a column in these matrices represents a latent feature vector for a single user or item, respectively, and R ≈ U T M for the known values in R. Once the user and item matrices are computed, they are used to build the recommendation system and predict the unknown ratings. Computing user and item matrices is the first step in building a variety of recommendation systems, so it is important to compute the factorization of R quickly.
The matrix factorization problem is closely related to the singular value decomposition (SVD), but the SVD of a matrix with missing values is undefined. However, since R ≈ U T M, one way to find the user and item matrices is by minimizing the squared difference between the approximated and actual value of the known ratings in R. Minimizing this difference is typically done by one of two algorithms: stochastic gradient descent (SGD) or alternating least squares (ALS) [8], [9]. ALS can be easily parallelized and can efficiently handle models that incorporate implicit data (e.g. the frequency of a user's mouse-clicks or time spent on a website) [10], but it is wellknown that ALS can require a large number of iterations to converge. Thus, in this paper, we propose an approach to significantly accelerate the convergence of the ALS algorithm for computing the user and item matrices. We use a nonlinear optimization algorithm, specifically the nonlinear conjugate gradient (NCG) algorithm [11], as a wrapper around ALS to significantly accelerate the convergence of ALS, and thus refer to this combined algorithm as ALS-NCG. Alternatively, the algorithm can be viewed as a nonlinearly preconditioned NCG method with ALS as the nonlinear preconditioner [12]. Our approach for accelerating the ALS algorithm using the NCG algorithm can be situated in the context of recent research activity on nonlinear preconditioning for nonlinear iterative solvers [13], [14], [12], [15], [16]. Some of the ideas date back as far as the 1960s [17], [18], but they are not well-known and remain under-explored experimentally and theoretically [13].
Parallel versions of collaborative filtering and recommendation are of great interest in the era of big data [19], [20], [21]. For example, the Apache Spark data processing environment [22] contains a parallel implementation of ALS for the collaborative filtering model of [23], [8]. In [19] an advanced distributed SGD method is described in Hadoop environments, followed by work in [20] that considers algorithms based on ALS and SGD in environments that use the Message Passing Interface (MPI). Scalable coordinate descent approaches for collaborative filtering are proposed in [21], also using the MPI framework. In addition to producing algorithmic advances, these papers have shown that the relative performance of the methods considered is often strongly influenced by the performance characteristics of the parallel computing paradigm that is used to implement the algorithms (e.g., Hadoop or MPI). In this paper we implement our proposed ALS-NCG algorithm in parallel using the Apache Spark framework, which is a largescale distributed data processing environment that builds on the principles of scalability and fault tolerance that are instrumental in the success of Hadoop and MapReduce. However, Spark adds crucial new capabilities in terms of in-memory computing and data persistence for iterative methods, and makes it possible to implement more elaborate algorithms with significantly better performance than is feasible in Hadoop. As such, Spark is already being used extensively for advanced big data analytics in the commercial setting [24]. The parallel Spark implementation of ALS for the collaborative filtering model of [23], [8] forms the starting point for applying the parallel acceleration methods proposed in this paper.
Our contributions in this paper are as follows. The specific optimization problem is formulated in Section II, and the accelerated ALS-NCG algorithm is developed in Section III. In Section IV, we study the convergence enhancements of ALS-NCG compared to standalone ALS in small serial tests using subsets of the MovieLens 20M dataset [25]. Section V describes our parallel implementation in Spark 1 , and Section VI contains results from parallel performance tests on a highend computing cluster with 16 nodes and 256 cores, using both the full MovieLens 20M dataset and a large synthetic dataset sampled from the MovieLens 20M data with up to 6 million users and 800 million ratings. We find that ALS-NCG converges significantly faster than ALS in both the serial and distributed Spark settings, demonstrating the overall speedup provided by our algorithmic acceleration. 1 Source code available: https://github.com/mbhynes/als-ncg II. PROBLEM DESCRIPTION The acceleration approach we propose in this paper is applicable to a broad class of optimization methods and collaborative filtering models. For definiteness, we choose a specific latent factor model, the matrix factorization model from [8] and [23], and a specific optimization method, ALS [23]. Given the data we use, the model is presented in terms of users and movies instead of the more generic users and items framework.
Let the matrix of user-movie rankings be represented by R = {r ij } nu×nm where r ij is the rating given to movie j by user i, n u is the number of users, and n m is the number of items. Note that for any user i and movie j, the value of r ij is either a real number or is missing, and in practice very few values are known. For example, the MovieLens 20M dataset [25] with 138,493 users and 27,278 movies contains only 20 million rankings, accounting for less than 1% of the total possible rankings. In the low-rank factorization of R with n f factors, each user i is associated with a vector u i ∈ R n f (i = 1, . . . , n u ), and each movie j is associated with a vector m j ∈ R n f (j = 1, . . . , n m ). The elements of m j measure the degree that movie j possesses each factor or feature, and the elements of u i similarly measures the affinity of user i for each factor or feature. The dot product u T i m j thus captures the interaction between user i and movie j, approximating user i's rating of movie j as r ij ≈ u T i m j . Denoting U = [u i ] ∈ R n f ×nu as the user feature matrix and M = [m j ] ∈ R n f ×nm as the movie feature matrix, our goal is to determine U and M such that R ≈ U T M by minimizing the following squared loss function: where I is the index set of known r ij in R, n ui denotes the number of ratings by user i, and n mj is the number of ratings of movie j. The term λ( i n ui u i 2 + j n mj m j 2 ) is a Tikhonov regularization [26] term commonly included in the loss function to prevent overfitting. The full optimization problem can be stated as

A. Alternating Least Squares Algorithm
The optimization problem in (2) is not convex. However, if we fix one of the unknowns, either U or M, then the optimization problem becomes quadratic and we can solve for the remaining unknown as a least squares problem. Doing this in an alternating fashion is the central idea behind ALS.
Consider the first step of the ALS algorithm in which M is fixed. We can determine the least squares solution to (1) for each u i by setting all the components of the gradient of (1) related to u i to zero: ∂L λ ∂u ki = 0 ∀ i, k where u ki is an element of U. Expanding the terms in the derivative of (1), we have where m kj is an element of M. In vector form, the resultant linear system for any u i is where I is the n f × n f identity matrix, I i is the index set of movies user i has rated, and M Ii represents the sub-matrix of M where columns j ∈ I i are selected. Similarly, R(i, I i ) is a row vector that represents the ith row of R with only the columns in I i included. The explicit solution for u i is then given by where . The analogous solution for the columns of M is found by fixing U, where each m j is given by where Here, I j is the index set of users that have rated movie j, U Ij represents the sub-matrix of U ∈ R n f ×nu where columns i ∈ I j are selected, and R(I j , j) is a column vector that represents the jth column of R with only the rows in I j included.
Algorithm 1 summarizes the ALS algorithm used to solve the optimization problem given in (2). From (3) and (4), we note that each of the columns of U and M may be computed independently; thus ALS may be easily implemented in parallel, as in [23]. However, the ALS algorithm can require many iterations for convergence, and we now propose an acceleration method for ALS that can be applied to both the parallel and serial versions.

B. Accelerated ALS Algorithm
In this section, we develop the accelerated ALS-NCG algorithm to solve the collaborative filtering optimization problem  (1). In practice we found that the nonlinear conjugate gradient algorithm by itself was very slow to converge to a solution of (2), and thus do not consider it as an alternative to the ALS algorithm. Instead we propose using NCG to accelerate ALS, or, said differently, the combined algorithm uses ALS as a nonlinear preconditioner for NCG [12].
The standard NCG algorithm is a line search algorithm in continuous optimization that is an extension of the CG algorithm for linear systems. While the linear CG algorithm is specifically designed to minimize the convex quadratic where A ∈ R n×n is a symmetric positive definite matrix, the NCG algorithm can be applied to general constrained optimization problems of the form min x∈R n f (x). Here, the minimization problem is (2), where the matrices U and M are found by defining the vector x ∈ R n f ×(nu+nm) as with function f (x) = L λ as in (1). The NCG algorithm generates a sequence of iterates x i , i ≥ 1, from the initial guess x 0 using the recurrence relation The parameter α k > 0 is the step length determined by a line search along search direction p k , which is generated by the following rule: where β k+1 is the update parameter and g k = ∇f (x k ) is the gradient of f (x) evaluated at x k . The update parameter β k+1 can take on various different forms. In this paper, we use the variant of β k+1 developed by Polak and Ribière [27]: Note that if a convex quadratic function is optimized using the NCG algorithm with an exact line search, then (7) reduces to the same β k+1 as in the original CG algorithm for linear systems [11]. The preconditioning of the NCG algorithm by the ALS algorithm modifies the expressions for β k+1 and p k+1 as follows, and is summarized in Algorithm 2. Let x k be the iterate generated from one iteration of the ALS algorithm applied to x k , x k = P (x k ), where P represents one iteration of ALS. This iterate is incorporated into the NCG algorithm by defining the preconditioned gradient direction generated by ALS as and replacing g k with g k in (6). Note that −g k is expected to be a descent direction that is an improvement compared to the steepest descent direction −g k . The update parameter β k+1 is redefined as β k+1 with the form (8) is similar to (7), however not every instance of g k has been replaced with g k . The specific form for β k+1 is chosen because if Algorithm 2 were applied to the convex quadratic problem with an exact line search and P were represented by a preconditioning matrix P, then Algorithm 2 would be equivalent to preconditioned CG with preconditioner P. In [12], an overview and an in-depth analysis are given of the different possible forms for β k+1 , however (8) performed best in our numerical experiments. Note that the primary computational cost in Algorithm 2 comes from computing both the ALS iteration, P (x k ), and α k , the step length parameter using a line search.

IV. SERIAL PERFORMANCE OF ALS-NCG
The ALS and ALS-NCG algorithms were implemented in serial in MATLAB, and evaluated using the MovieLens 20M dataset [25]. The entire MovieLens 20M dataset has 138,493 users, 27,278 movies, and just over 20 million ratings, where each user has rated at least 20 movies. To investigate the algorithmic performance as a function of problem size, both the ALS and ALS-NCG algorithms were run on subsets of the MovieLens 20M dataset. In creating subsets, we excluded outlier users with either very many or very few movie ratings relative to the median number of ratings per user. To construct a subset with n u users, the users from the full dataset were sorted in descending order by the values of n ui for each user. Denoting the index of the user with the median number of ratings by c, the set of users from index c − nu 2 to c+( nu 2 −1) were included. Once the users were determined, the same process was used to select the movies, where the ratings per movie were computed only for the chosen users.
All serial experiments were performed on a linux workstation with a quad-core 3.16 GHz processor (Xeon X5460) and 8 GB of RAM. For the ALS-NCG algorithm, the Moré-Thuente line search algorithm from the Poblano toolbox [28] was used to compute α k . The line search parameters were as follows: 10 −4 for the sufficient decrease condition tolerance, 10 −2 for the curvature condition tolerance, and an initial step length of 1 and a maximum of 20 iterations. The stopping criteria for both ALS and ALS-NCG were the maximum number of iterations as well as a desired tolerance value in the gradient norm normalized by the number of variables, 1 N g k for N = n f × (n u + n m ). For both algorithms, a normalized gradient norm of less than 10 −6 was required within at most 10 4 iterations. In addition, ALS-NCG had a maximum number of allowed function evaluations in the line search equal to 10 7 .
The serial tests were performed on ratings matrices of 4 different sizes: n u × n m = 400 × 80, 800 × 160, 1600 × 320 and 3200 × 640, where n u is the number of users and n m is the number of movies. For each ratings matrix, ALS and ALS-NCG were used to solve the optimization problem in (2) with λ = 0.1 and n f = 10. The algorithms were each run using 20 different random starting iterates (which were the same for both algorithms) until one of the stopping criteria was reached. Table I summarizes the timing results for different problem sizes, where the given times are written in the form a ± b where a is the mean time in seconds and b is the standard deviation about the mean. Since computing the gradient is not explicitly required in the ALS algorithm, the computation time for the gradient norm was excluded from the timing results for the ALS algorithm. Runs that did not converge based on the gradient norm tolerance were not included in the mean and standard deviation calculations, however the only run that did not converge to 1 N g k < 10 −6 before reaching the maximum number of iterations was a single 1600 × 320 ALS run. The large standard deviations in the timing measurements stem from the variation in the number of iterations required to reduce the gradient norm. The fourth column of Table I shows the acceleration factor of ALS-NCG, computed as the mean time for convergence of ALS divided by the mean time for convergence of ALS-NCG. We see from this table that ALS-NCG significantly accelerates the ALS algorithm for all problem sizes. Similarly, Table II summarizes the number of iterations required to reach convergence for each algorithm. Again, the results were calculated based on converged runs only.
From Tables I and II it is clear that ALS-NCG accelerates the convergence of the ALS algorithm, using the gradient norm as the measure of convergence. However, since the factor matrices U and M are used to make recommendations, we would also like to examine the convergence of the algorithms in terms of the accuracy of the resultant recommendations. In particular, we are interested in the rankings of the top t movies (e.g. top 20 movies) for each user, and want to explore how these rankings change with increasing number of iterations for both ALS and ALS-NCG. If the rankings of the top t movies for each user no longer change, then the algorithm has likely computed an accurate solution.
To measure the relative difference in rankings we use a metric that is based on the number of pairwise swaps required to convert a vector of movie rankings p 2 into another vector p 1 , but only for the top t movies. We use a modified Kendall-Tau [29] distance to compute the difference between ranking vectors based only on the rankings of the top t items, normalize the distance to range in [0, 1], and subsequently average the distances over all users. To illustrate the ranking metric, consider the following example, where p 1 = [6, 3, 1, 2, 4, 5], and p 2 = [3, 4, 2, 5, 6, 1] are two different rankings of movies for user i, and let t = 2. We begin by finding the top t movies in p 1 . Here, user i ranks movie 6 highest. In p 2 , movie 6 is ranked 5th, and there are 4 pairwise inversions required to place movie 6 in the first component of p 2 , producing a new ranking vector for the second iteration,p 2 = [6, 3, 4, 2, 5, 1]. The 2nd highest ranked movie in p 1 is movie 3. Inp 2 , movie 3 is already ranked 2nd; as such, no further inversions are required. In this case, the total number of pairwise swaps required to match p 2 to p 1 for the top 2 rankings is s i = 4. The distance s i is normalized to 1 if no inversions were needed (i.e. p 1 = p 2 in the top t spots) and 0 if the maximum number of inversions are needed. The maximum number of inversions occurs if p 2 and p 1 are in opposite order, requiring n m − 1 inversions in the first step, n m − 2 inversions in the second step, and so on until n m − t inversions are needed in the t-th step. Thus, the maximum number of inversions is s max = (n m −1)+(n m −2)+. . .+(n m −t) = t 2 (2n m −t−1), yielding the ranking accuracy metric for user i as The total ranking accuracy for an algorithm is taken as the average value of q i across all users relative to ranking of the top t movies for each user produced from the solution obtained after the algorithm converged. The normalized Kendall-Tau ranking metric for the top t rankings described above was used to evaluate the accuracy of ALS and ALS-NCG as a function of running time for t = 20. Tables III and IV summarize the time needed for each   TABLE III   algorithm to reach a specified percentage ranking accuracy for ratings matrices with different sizes. Both algorithms were run with 20 different random starting iterates, and the time to reach a given ranking accuracy is written in the form a±b where a is the mean time across all converged runs and b is the standard deviation about the mean. In the fourth column of Tables III and IV, we have computed the acceleration factor of ALS-NCG as the ratio of the mean time for ALS to the mean time for ALS-NCG. To reach 70% accuracy, ALS is faster for both problem sizes, however for accuracies greater than 70%, ALS-NCG converges in significantly less time than ALS. Thus, if an accurate solution is desired, ALS-NCG is much faster in general than ALS. This is further illustrated in Fig. 1, which shows the average time needed for both algorithms to reach a given ranking accuracy for the 400 × 80 ratings matrix. Here, both ALS and ALS-NCG reach a ranking accuracy of 75% in less than a second, however it then takes ALS approximately 100 s to increase the ranking accuracy from 75% to 100%, while ALS-NCG reaches the final ranking in only 14 s.

A. Apache Spark
Apache Spark is a fault-tolerant, in-memory cluster computing framework designed to supersede MapReduce by maintaining program data in memory as much as possible between distributed operations. The Spark environment is built upon two components: a data abstraction, termed a resilient distributed dataset (RDD) [22], and the task scheduler, which uses a delay scheduling algorithm [30]. We describe the fundamental aspects of RDDs and the scheduler below.
Resilient Distributed Datasets. RDDs are immutable, distributed datasets that are evaluated lazily via their provenance information-that is, their functional relation to other RDDs or datasets in stable storage. To describe an RDD, consider an immutable distributed dataset D of k records with homogeneous type: D = k i d i with d i ∈ D. The distribution of D across a computer network of nodes {v α }, such that d i is stored in memory or on disk on node v α , is termed its partitioning according to a partition function P (d i ) = v α . If D is expressible as a finite sequence of deterministic operations on other datasets D 1 , . . . , D l that are either RDDs or persistent records, then its lineage may be written as a directed acyclic graph L formed with the parent datasets {D i } as the vertices, and the operations along the edges. Thus, an RDD of type D (written RDD [D]) is the tuple (D, P, L).
Physically computing the records {d i } of an RDD is termed its materialization, and is managed by the Spark scheduler program. To allocate computational tasks to the compute nodes, the scheduler traverses an RDD's lineage graph L and divides the required operations into stages of local computations on parent RDD partitions. Suppose that R 0 = ( i x i , P 0 , L 0 ) were an RDD of numeric type RDD [R], and let R 1 = ( i y i , P 1 , L 1 ) be the RDD resulting from the application of function f : R → R to each record of R 0 . To compute {y i }, R 1 has only a single parent in the graph L 1 , and hence the set of tasks to perform is {f (x i )}. This type of operation is termed a map operation. If P 1 = P 0 , L 1 is said to have a narrow dependency on R 0 : each y i may be computed locally from x i , and the scheduler would allocate the task f (x i ) to a node that stores x i .
Stages consist only of local map operations, and are bounded by shuffle operations that require communication and data transfer between the compute nodes. For example, shuffling is necessary to perform reduce operations on RDDs, wherein a scalar value is produced from an associative binary operator applied to each element of the dataset. In implementation, a shuffle is performed by writing the results of the tasks in the preceding stage, {f (x i )}, to a local file buffer. These shuffle files may or may not be written to disk, depending on the operating system's page table, and are fetched by remote nodes as needed in the subsequent stage.
Delay Scheduling and Fault Tolerance. The simple delay scheduling algorithm [30] prioritizes data locality when submitting tasks to the available compute nodes. If v α stores the needed parent partition x i to compute task f (x i ), but is temporarily unavailable due to faults or stochastic delays, rather than submitting the task on another node, the scheduler will wait until v α is free. However, if v α does not become available within a specified maximum delay time (several seconds in practice), the scheduler will resubmit the tasks to a different compute node. However, as x i is not available in memory on the different node, the lineage L 0 of the RDD R 0 must be traversed further, and the tasks required to compute x i from the parent RDDs of R 0 will be submitted for computation in addition to f (x i ). Thus, fault tolerance is achieved in the system through recomputation.

B. ALS Implementation in Spark
The Apache Spark codebase contains a parallel implementation of ALS for the collaborative filtering model of [23],  [8]; see [24]. We briefly outline the main execution and relevant optimizations of the existing ALS implementation. The implementation stores the factor matrices U and M as single precision column block matrices, with each block as an RDD partition. A U column block stores factor vectors for a subset of users, and an M column block stores factor vectors for a subset of movies. The ratings matrix is stored twice: both R and R T are stored in separate RDDs, partitioned in row blocks (i.e., R is partitioned by users and R T by movies). Row blocks of R and R T are stored in a compressed sparse matrix format. The U and R RDDs have the same partitioning, such that a node that stores a partition of U also stores the partition of R such that the ratings for each user in its partitions of U are locally available. When updating a U block according to (3), the required ratings are available in a local R block, but the movie factor vectors in M corresponding to the movies rated by the users in a local U block must be shuffled across the network. These movie factors are fetched from different nodes, and, as explained below, an optimized routing table strategy is used from [24] that avoids sending duplicate information. Similarly, updating a block of M according to (4) uses ratings data stored in a local R T block, but requires shuffling of U factor vectors using a second routing table.
Block Partitioning. All RDDs are partitioned into n b partitions, where n b is an integer multiple of the number of available compute cores in practice. For example, M is divided into column blocks M j b with block (movie) index j b ∈ {0, . . . , n b − 1} by hash partitioning the movie factor vectors such that m j ∈ R j b if j ≡ j b (mod n b ) as in Fig. 2. Similarly, U is hash partitioned into column blocks U i b with block (user) index i b ∈ {0, . . . , n b − 1}. The RDDs for M and U can be taken as type where the blocks are tracked by the indices j b and i b . R is partitioned by rows (users) into blocks with type RDD [(i b , R i b )] with the same partitioning as the RDD representing U (and similarly for the R T and M RDDs). By sharing the same user-based partitioning scheme, the blocks R i b and U i b are normally located on the same compute node, except when faults occur. The same applies to R T j b and M j b due to the movie-based partitioning scheme.
Routing  Fig. 3. Schematic of the use of the routing table Tm(m j ) in the Spark ALS shuffle. In (a) the blocks M j b are filtered using Tm(m j ) for each destination R i b and shuffled to the respective blocks in (b), where arrows between the shaded backgrounds represent network data transfer between different compute nodes. In (c), when updating block U i b , the ratings information is locally available in R i b . data shuffling in ALS. Suppose we want to update the user factor block U i b according to (3). The required ratings data, R i b , is stored locally, but a global shuffle is required to obtain all movie factor vectors in M that correspond to the movies rated by the users in U i b . To optimize the data transfer, a routing table T m (m j ) is constructed by determining, for each of the movie factor blocks M j b , which factor vectors have to be sent to each partition of R i b (that may reside on other nodes). In Fig. 3 (a), the blocks M j b are filtered using T m (m j ) such that a given m j ∈ M j b is written to the buffer destined j b , only once regardless of how many u i in U i b have ratings for movie j, and only if there is at least one u i in U i b that has rated movie j. This is shown by the hatching of each m j vector in Fig. 3 (a); for instance, the first column in M 1 is written only to M [1] 1 and has one set of hatching lines, but the last column is written to both M and correspondingly has two sets of hatching lines. Once the buffers are constructed, they are shuffled to the partitions of R, as in Fig. 3 (b) such that both the movie factors and ratings are locally available to compute the new U i b block, as in Fig. 3 (c). The routing table formulation for shuffling U, with mapping T u (u i ), is analogous. Note that the routing tables are constructed before the while loop in Algorithm 1, and hence do not need recomputation in each iteration. Fig. 3 (c), a compute node that stores R i b , will obtain n b buffered arrays of filtered movie factors. Once the factors have been shuffled, A j is computed for each u i as j∈Ii m j m T j +n ui I, and v j as j∈Ii r ij m j using the Basic Linear Algebra Subprograms (BLAS) library [31]. The resulting linear system for u i is then solved via the Cholesky decomposition using LAPACK routines [32], giving the computation to update U an asymptotic complexity of O(n u n 3 f ) since n u linear systems must be solved. Solving for m j is an identical operation with the appropriate routing table and has O(n m n 3 f ) complexity.

C. ALS-NCG Implementation in Spark
We now discuss our contributions in parallelizing ALS-NCG for Spark. Since the calculation of α k in Algorithm 2 requires a line search, the main technical challenges we address are how to compute the loss function in (1) and its gradient in an efficient way in Spark, obtaining good parallel performance. To this end, we formulate a backtracking line search procedure that dramatically reduces the cost of multiple function evaluations, and we take advantage of the routing table mechanism to obtain fast communication. We also extend the ALS implementation in Spark to support the additional NCG vector operations required in Algorithm 2 using BLAS routines.
Vector Storage. The additional vectorsx, g,ḡ, and p were each split into two separate RDDs, such that blocks corresponding to the components of u i (see (5)) were stored in one RDD and partitioned in the same way as U with block index i b . Analogously, blocks corresponding to components of m j were stored in another RDD, partitioned in the same way as M with block index j b . This ensured that all vector blocks were aligned component-wise for vector operations; furthermore, the p blocks could also be shuffled efficiently using the routing tables in the line search (see below).
Vector Operations. RDDs have a standard operation termed a join, in which two RDDs R 1 and R 2 representing the datasets of tuples i (k i , a i ) and i (k i , b i ) with type (K, A) and (K, B), respectively, are combined to produce an RDD R 3 of type RDD [ (K, (A, B))], where R 3 = R 1 .join(R 2 ) represents the dataset l (k l , (a l , b l )) of combined tuples and k l is a key common to both R 1 and R 2 . Parallel vector operations between RDD representations of blocked vectors x and y were implemented by joining the RDDs by their block index i b and calling BLAS level 1 interfaces on the vectors within the resultant tuples of aligned vector blocks, Since the RDD implementations of vectors had the same partitioning schemes, this operation was local to a compute node, and hence very inexpensive. One caveat, however, is that BLAS subprograms generally perform modifications to vectors in place, overwriting their previous components. For fault tolerance, RDDs must be immutable; as such, whenever BLAS operations were required on the records of an RDD, an entirely new vector was allocated in memory and the RDD's contents copied to it. Algorithm 3 shows the operations required for the vector addition ax+ by using the BLAS.axpby routine (note that the result overwrites y in place). The inner product of two block vector RDDs and norm of a single block vector RDD were implemented in a similar manner to vector addition, with an additional reduce operation to sum the scalar componentwise dot products. These vector operations were used to computeβ k+1 in (8).
Line Search & Loss Function Evaluation. We present a computationally cheap way to implement a backtracking line search in Spark for minimizing (1). A backtracking line search minimizes a function f (x) along a descent direction p by decreasing the step size in each iteration by a factor of Algorithm 3: RDD Block Vector BLAS.axpby ij in the summation for each coefficient only require level 1 BLAS operations between the block vectors x ui , p ui , x mj , and p mj for known (i, j) ∈ I. Thus, coefficients of Q(α) were computed at the beginning of the line search with a single shuffle operation using the routing table T m (m j ) to match vector pairs with dot products contributing to the coefficients. Here, T m (m j ) was chosen since there is far less communication required to shuffle {m j }, as n u ≫ n m . Since each iteration of the line search was very fast after computing the coefficients of Q(α), we used relatively large values of τ = 0.9, c = 0.5, and α 0 = 10 in Algorithm 4 that searched intensively along direction p.

Gradient Evaluation.
We computed g k with respect to a block for u i using only BLAS level 1 operations as 2λn ui u i + 2 j∈Ii m j (u T i m j − r ij ), with an analogous operation with respect to each m j . As this computation requires matching the u i and m j factors, the routing tables T u (u i ) and T m (m j ) were used to shuffle u i and m j consecutively. Evaluating g k can be performed with O(n f (n u i n ui + n m j n mj )) operations, but requires two shuffles. As such, the gradient computation required as much communication as a single iteration of ALS in which both U and M are updated.

VI. PARALLEL PERFORMANCE OF ALS-NCG
Our comparison tests of the ALS and ALS-NCG algorithms in Spark were performed on a computing cluster composed of 16 homogeneous compute nodes, 1 storage node hosting a network filesystem, and 1 head node. The nodes were interconnected by a 10 Gb ethernet managed switch (Pow-erConnect 8164). Each compute node was a 64 bit rack server (PowerEdge R620) running Ubuntu 14.04, with linux kernel 3.13. The compute nodes all had two 8-core 2.60 GHz chips (Xeon E5-2670) and 256 GB of SDRAM. The head node had the same processors as the compute nodes, but had 512 GB of RAM. The single storage node (PowerEdge R720) contained two 2 GHz processors, each with 6 cores (Xeon E5-2620), 64 GB of memory, and 12 hard disk drives of 4 TB capacity and 7200 RPM nominal speed. Finally, compute nodes were equipped with 6 ext4-formatted local SCSI 10k RPM hard disk drives, each with a 600 GB capacity.
Our Apache Spark assembly was built from a snapshot of the 1.3 release using Oracle's Java 7 distribution, Scala 2.10, and Hadoop 1.0.4. Input files to Spark programs were stored on the storage node in plain text. The SCSI hard drives on the compute nodes' local filesystems were used as Spark spilling and scratch directories, and the Spark checkpoint directory for persisting RDDs was specified in the network filesystem hosted by the storage node, and accessible to all nodes in the cluster. Shuffle files were consolidated into larger files, as recommended for ext4 filesystems [33]. In our experiments, the Spark master was executed on the head node, and a single instance of a Spark executor was created on each compute node. It was empirically found that the ideal number of cores to make available to the Spark driver was 16 per node, or the number of physical cores for a total of 256 available cores. The value of n b was set to the number of cores in all experiments.
To compare the performance of ALS and ALS-NCG in Spark, the two implementations were tested on the Movie-Lens 20M dataset with λ = 0.01 and n f = 100. In both algorithms, the RDDs were checkpointed to persistent storage every 10 iterations, since it is a widely known issue in the Spark community that RDDs with very long lineage graphs cause stack overflow errors when the scheduler recursively traverses their lineage [34]. Since RDDs are materialized via lazy evaluation, to obtain timing measurements, actions were triggered to physically compute the partitions of each block vector RDD at the end of each iteration in Algorithm 1 and 2. For each experimental run of ALS and ALS-NCG, two experiments with the same initial user and movie factors were performed: in one, the gradient norm in each iteration was computed and printed (incurring additional operations); in the other experiment, no additional computations were performed such that the elapsed times for each iteration were correctly measured. Fig. 4 shows the convergence in gradient norm, normalized by the degrees of freedom N = n f × (n u + n m ), for six separate runs of ALS and ALS-NCG on 8 compute nodes for the MovieLens 20M dataset. The subplots (a) and (b) show  N g k over 100 iterations and 25 minutes, respectively. This time frame was chosen since it took just over 20 minutes for ALS-NCG to complete 100 iterations; note that in Fig. 4 (b), 200 iterations of ALS are shown, since with this problem size it took approximately twice as long to run a single iteration of ALS-NCG. The shaded regions in Fig. 4 show the standard deviation about the mean value for 1 N g k across all runs, computed for non-overlapping windows of 3 iterations for subplot (a) and 30 s for subplot (b). Even within the uncertainty bounds, ALS-NCG requires much less time and many fewer iterations than ALS to reach accurate values of 1 N g k (e.g. below 10 −3 ).
The operations that we have implemented in ALS-NCG that are additional to standard ALS in Spark have computational complexity that is linear in problem size. To verify the expected linear scaling, experiments with a constant number of nodes and increasing problem size up to 800 million ratings were performed on 16 compute nodes. Synthetic ratings matrices were constructed by sampling the MovieLens 20M dataset such that the synthetic dataset had the same statistical sparsity, realistically simulating the data transfer patterns in each iteration. To do this, first the dimension n u of the sampled dataset was fixed, and for each user n ui was sampled from the empirical probability distribution p(n ui |R) of how many movies each user ranked, computed empirically from the MovieLens 20M dataset. The n ui movies were then sampled from the empirical likelihood of sampling the jth movie, p(m j |R), and the resultant rating value was sampled from the distribution of numerical values for all ratings. The choice of scaling up the users (for a fixed set of movies) was made Linear scaling of computation time per iteration with increasing nu on a synthetic dataset for ALS and ALS-NCG on 16 compute nodes (256 cores) for up to 6M users, corresponding to 800M ratings.
to model the situation in which an industry's user base grows far more rapidly than its items. Fig. 5 shows the linear scaling in computation time for both ALS and ALS-NCG. The values shown are average times per iteration over 50 iterations, for n u from 1 to 6 million, corresponding to the range from 133 to 800 million ratings. The error bars show the uncertainty in this measurement, where the standard deviation takes into account that iterations with and without checkpointing come from two different populations with different average times. The uncertainty in the time per iteration for ALS-NCG is larger due to the greater overhead of memory management and garbage collection by the Java Virtual Machine required with more RDDs. While each iteration of ALS-NCG takes longer due to the additional line search and gradient computations, we note that many fewer iterations are required to converge.
Finally, we compute the relative speedup that was attainable on the large synthetic datasets. For the value of 1 N g k in each iteration of ALS-NCG, we determined how many iterations of regular ALS were required to achieve an equal or lesser value gradient norm. Due to the local variation in 1 N g k (as in Fig. 4), a moving average filter over every two iterations was applied to the ALS-NCG gradient norm values. The total time required for ALS and ALS-NCG to reach a given gradient norm was then estimated from the average times per iteration in Fig. 5. The ratios of these total times for ALS and ALS-NCG are shown in Fig. 6 as the relative speedup factor for the 1M, 3M, and 6M users ratings matrices. When an accurate solution is desired, ALS-NCG often achieves faster convergence by a factor of 3 to 5, with the acceleration factor increasing with greater desired accuracy in the solution.
VII. CONCLUSION In this paper, we have demonstrated how the addition of a nonlinear conjugate gradient wrapper can accelerate the convergence of ALS-based collaborative filtering algorithms when accurate solutions are desired. Furthermore, our parallel ALS-NCG implementation can significantly speed up big data recommendation in the Apache Spark distributed computing environment, and the proposed NCG acceleration can be naturally extended to Hadoop or MPI environments. It may also be applicable to other optimization methods for collaborative filtering. Though we have focused on a simple latent factor model, our acceleration can be used with any collaborative filtering model that uses ALS. We expect that our acceleration approach will be especially useful for advanced collaborative filtering models that achieve low root mean square error (RMSE), since these models require solving the optimization problem accurately, and that is precisely where accelerated ALS-NCG shows the most benefit over standalone ALS.