Randomized online matrix factorization
Sparse matrix factorization is a popular tool to obtain interpretable data decompositions, which are also effective to perform data completion or denoising. Its applicability to large datasets has been addressed with online and randomized methods, that reduce the complexity in one of the matrix dimension, but not in both of them. In this paper, we tackle very large matrices in both dimensions. We propose a new factoriza-tion method that scales gracefully to terabyte-scale datasets, that could not be processed by previous algorithms in a reasonable amount of time. We demonstrate the efficiency of our approach on massive functional Magnetic Resonance Imaging (fMRI) data, and on matrix completion problems for recommender systems, where we obtain significant speed-ups compared to state-of-the art coordinate descent methods.READ FULL TEXT VIEW PDF
We present a matrix-factorization algorithm that scales to input matrice...
Matrix factorization is a key tool in data analysis; its applications in...
Several complex tasks that arise in organizations can be simplified by
If learning methods are to scale to the massive sizes of modern datasets...
Muscle coordination patterns of lingual behaviors are synergies generate...
Multiresolution analysis and matrix factorization are foundational tools...
We present a family of online algorithms for real-time factorization-bas...
Randomized online matrix factorization
Randomized online matrix factorization
In this section, we introduce dictionary learning as a matrix factorization problem, and present stochastic algorithms that observe one column (or a minibatch) at every iteration.
The goal of matrix factorization is to decompose a matrix – typically signals of dimension – as a product of two smaller matrices:
with potential sparsity or structure requirements on and . In statistical signal applications, this is often a dictionary learning problem, enforcing sparse coefficients . In such a case, we call the “dictionary” and the sparse codes. We use this terminology throughout the paper.
Learning the dictionary is typically performed by minimizing a quadratic data-fitting term, with constraints and/or penalties over the code and the dictionary:
where is a convex set of , and a is a penalty over the code, to enforce structure or sparsity. In large and large settings, typical in recommender systems, this problem is solved via block coordinate descent, which boils down to alternating least squares if regularizations on and are quadratic (Hastie et al., 2014).
The constraint set is traditionally a technical constraint ensuring that the coefficients do not vanish, making the effect of the penalty disappear. However, other constraints can also be used to enforce sparsity or structure on the dictionary (see Varoquaux et al., 2013). In our paper, is the Cartesian product of a or norm ball:
where and or . The choice of and typically offers some flexibility in the regularization effect that is desired for a specific problem; for instance, classical dictionary learning uses and , leading to sparse coefficients , whereas our experiments on fMRI uses and , leading to sparse dictionary elements that can be interpreted as brain activation maps.
In stochastic optimization, the number of signals is assumed to be large (or potentially infinite), and the dictionary can be written as a solution of
where the signals
are assumed to be i.i.d. samples from an unknown probability distribution. Based on this formulation,Mairal et al. (2010) have introduced an online dictionary learning approach that draws a single signal at iteration (or a minibatch), and computes its sparse code using the current dictionary according to
Then, the dictionary is updated by approximately minimizing the following surrogate function
which involves the sequence of past signals and the sparse codes that were computed in the past iterations of the algorithm. The function is called a “surrogate” in the sense that it only approximates the objective . In fact, it is possible to show that it converges to a locally tight upper-bound of the objective, and that minimizing at each iteration asymptotically provides a stationary point of the original optimization problem. The underlying principle is that of majorization-minimization, used in a stochastic fashion (Mairal, 2013).
One key to obtain efficient dictionary updates is the observation that the surrogate can be summarized by a few sufficient statistics that are updated at every iteration. In other words, it is possible to describe without explicitly storing the past signals and codes for . Indeed, we may define two matrices and
and the surrogate function is then written:
The gradient of can be computed as
Minimization of is performed using block coordinate descent on the columns of . In practice, the following updates are successively performed by cycling over the dictionary elements for
where Proj denotes the Euclidean projection over the constraint norm constraint . It can be shown that this update corresponds to minimizing with respect to when fixing the other dictionary elements (see Mairal et al., 2010).
Factorization of matrices with missing value have raised a significant interest in signal processing and machine learning, especially as a solution for recommender systems. In the context of dictionary learning, a similar effort has been made bySzabó et al. (2011) to adapt the framework to missing values. Formally, a mask , represented as a binary diagonal matrix in , is associated with every signal , such that the algorithm can only observe the product at iteration instead of a full signal . In this setting, we naturally derive the following objective
where the pairs are drawn from the (unknown) data distribution. Adapting the online algorithm of Mairal et al. (2010) would consist of drawing a sequence of pairs , and building the surrogate
where is the size of the mask and
Unfortunately, this surrogate cannot be summarized by a few sufficient statistics due to the masks : some approximations are required. This is the approach chosen by Szabó et al. (2011). Nevertheless, the complexity of their update rules is linear in the full signal size , which makes them unadapted to the large- regime that we consider.
Using the formalism exposed above, we now consider the problem of factorizing a large matrix in into two factors in and in with the following setting: both and are large (greater than up to several millions), whereas is reasonable (smaller than and often near ), which is not the standard dictionary-learning setting; some entries of may be missing. Our objective is to recover a good dictionary taking into account appropriate regularization.
To achieve our goal, we propose to use an objective akin to (11
), where the masks are now random variablesindependant from the samples. In other words, we want to combine ideas of online dictionary learning with random subsampling, in a principled manner. This leads us to consider an infinite stream of samples , where the signals are i.i.d. samples from the data distribution – that is, a column of selected at random – and “selects” a random subset of observed entries in . This setting can accommodate missing entries, never selected by the mask, and only requires loading a subset of at each iteration.
The main justification for choosing this objective function is that in the large sample regime that we consider, computing the code using only a random subset of the data according to (13
) is a good approximation of the code that may be computed with the full vectorin (5). This of course requires choosing a mask that is large enough; in the fMRI dataset, a subsampling factor of about – that is only of the entries of are observed – resulted in a similar speed-up (see experimental section) to achieve the same accuracy as the original approach without subsampling. This point of view also justifies the natural scaling factor introduced in (11).
An efficient algorithm must address two challenges: (i) performing dictionary updates that do not depend on but only on the mask size; (ii) finding an approximate surrogate function that can be summarized by a few sufficient statistics. We provide a solution to these two issues in the next subsections and present the method in Algorithm 1.
with the same matrix as in (8), which is updated as
and to replace in (8) by the matrix
which is the same as (7) when . Since is a diagonal matrix, is also diagonal and simply “counts” how many times a row has been seen by the algorithm. thus behaves like for large , as in the fully-observed algorithm. By design, only rows of selected by the mask differ from . The update can therefore be achieved in operations:
This only requires keeping in memory the diagonal matrix , and updating the rows of selected by the mask. All operations only depend on the mask size instead of the signal size .
With a surrogate function in hand, we now describe how to update the codes and the dictionary when only partial access to data is possible. The complexity for computing the sparse codes is obviously independent from since (13) consists in solving a reduced
penalized linear regression ofin on in . Thus, we focus here on dictionary update rules.
The naive dictionary update (18) has complexity due to the matrix-vector multiplication for computing . Reducing the single iteration complexity of a factor requires reducing the dimensionality of the dictionary update phase. We propose two strategies to achieve that, both using block coordinate descent, by considering
where is the partial derivative of with respect to the -th column and rows selected by the mask.
The update (18) represents a classical block coordinate descent step involving particular blocks. Following Mairal et al. (2010), we perform one cycle over the columns warm-started on . Formally, the gradient step without projection for the -th component consists of updating the vector
where are the -th columns of , respectively. The update has complexity since it only involves rows of and only entries of have changed.
Block coordinate descent algorithms require orthogonal projections onto the constraint set . In our case, this amounts to the projection step on the unit ball corresponding to the norm in (18). The complexity of such a projection is usually both for and -norms (see Duchi et al., 2008). We consider here two strategies.
When , it is possible to perform the projection implicitly with complexity . The computational trick is to notice that the projection amounts to a simple rescaling operation
which may have low complexity if the dictionary elements are stored in memory as a product where is in and is a rescaling coefficient such that . We code the gradient step (19) followed by -ball projection by the updates
Note that the update of corresponds to the gradient step without projection (19) which costs , whereas the norm of is updated in operations. The computational complexity is thus independent of and the only price to pay is to rescale the dictionary elements on the fly, each time we need access to them.
The case of is slightly different but can be handled in a similar manner, by storing an additional scalar for each dictionary element . More precisely, we store a vector in such that , and a classical result (see Duchi et al., 2008) states that there exists a scalar such that
where is the soft-thresholding operator, applied elementwise to the entries of . Similar to the case , the “lazy” projection consists of tracking the coefficient for each dictionary element and updating it after each gradient step, which only involves coefficients. For such sparse updates followed by a projection onto the -ball, Duchi et al. (2008) proposed an algorithm to find the threshold in
operations. The lazy algorithm involves using particular data structures such as red-black trees and is not easy to implement; this motivated us to investigate another simple heuristic that also performs well in practice.
The heuristic consists in performing the projection by forcing the coefficients outside the mask not to change. This results in the orthogonal projection of each on , which is a subset of the original constraint set .
All the computations require only 4 matrices kept in memory , , , with additional , matrices and vectors for the exact projection case, as summarized in Alg. 1.
Our model is related to the classical -penalized matrix completion model (e.g. Bell & Koren, 2007) we rewrite
With quadratic regularization on and – that is, using and – (11) only differs in that it uses a penalization on instead of a constraint. Srebro et al. (2004) introduced the trace-norm regularization to solve a convex problem equivalent to (23). The major difference is that we adopt a non-convex optimization strategy, thus losing the benefits of convexity, but gaining on the other hand the possibility of using stochastic optimization.
Our algorithm can be slightly modified to use weights that differ from for and , as advocated by Mairal (2013). It also proves beneficial to perform code computation on mini-batches of masked samples. Update of the dictionary is performed on the rows that are seen at least once in the masks .
The proposed algorithm was designed to handle massive datasets: masking data enables streaming a sequence instead of , reducing single-iteration computational complexity and IO stress of a factor , while accessing an accurate description of the data. Hence, we analyze in detail how our algorithm improves performance for sparse decomposition of fMRI datasets. Moreover, as it relies on data masks, our algorithm is well suited for matrix completion, to reconstruct a data stream from the masked stream . We demonstrate the accuracy of our algorithm on explicit recommender systems and show considerable computational speed-ups compared to an efficient coordinate-descent based algorithm.
Matrix factorization has long been used on functional Magnetic Resonance Imaging (McKeown et al., 1998). Data are temporal series of 3D images of brain activity, to decompose in spatial modes capturing regions that activate together. The matrices to decompose are dense and heavily redundant, both spatially and temporally: close voxels and successive records are correlated. Data can be huge: we use the whole HCP dataset (Van Essen et al., 2013), with (2000 records, 1 200 time points) and , totaling 2 TB of dense data.
Interesting dictionaries for neuroimaging capture spatially-localized components, with a few brain regions. This can be obtained by enforcing sparsity on the dictionary: in our formalism, this is achieved with -ball projection for . We set , and . Historically, such decomposition have been obtained with the classical dictionary learning objective on transposed data (Varoquaux et al., 2013): the code holds sparse spatial maps and voxel time-series are streamed. However, given the size of for our dataset, this method is not usable in practice.
Handling such volume of data sets new constraints. First, efficient disk access becomes critical for speed. In our case, learning the dictionary is done by accessing the data in row batches, which is coherent with fMRI data storage: no time is lost seeking data on disk. Second, reducing IO load on the storage is also crucial, as it lifts bottlenecks that appear when many processes access the same storage at the same time, e.g. during cross-validation on within a supervised pipeline. Our approach reduces disk usage by a factor . Finally, parallel methods based on message passing, such as asynchronous coordinate descent, are unlikely to be efficient given the network / disk bandwidth that each process requires to load data. This makes it crucial to design efficient sequential algorithms.
We quantify the effect of random subsampling for sparse matrix factorization, in term of speed and accuracy. A natural performance evaluation is to measure an empirical estimate of the lossdefined in Eq. 4 from unseen data, to rule out any overfitting effect. For this, we evaluate on a test set . Pratically, we sample in a pseudo-random manner: we randomly select a record, from where we select a random batch of rows – we use a batch size of , empirically found to be efficient. We load in memory and perform an iteration of the algorithm. The mask sequence is sampled by breaking random permutation vectors into chunks of size .
Fig. 1(a) compares our algorithm with subsampling ratios in to vanilla online dictionary learning algorithm , plotting trajectories of the test objective against real CPU time. There is no obvious choice of due to the unsupervised nature of the problem: we use and , that bounds the range of providing interpretable dictionaries.
First, we observe the convergence of the objective function for all tested , providing evidence that the approximations made in the derivation of update rules does not break convergence for such . Fig. 1(b) shows the validity of the obtained dictionary relative to the reference output: both objective function and ratio – the relevant value to measure sparsity in our setting – are comparable to the baseline values, up to . For high regularization and , our algorithm tends to yield somewhat sparser solutions ( lower ) than the original algorithm, due to the approximate -projection we perform. Obtained maps still proves as interpretable as with baseline algorithm.
Our algorithm proves much faster than the original one in finding a good dictionary. Single iteration time is indeed reduced by a factor , which enables our algorithm to go over a single epoch times faster than the vanilla algorithm and capture the variability of the dataset earlier. To quantify speed-ups, we plot the empirical objective value of against the number of observed records in Fig. 3. For , increasing little reduces convergence speed per epoch: random subsampling does not shrink much the quantity of information learned at each iteration.
This brings a near speed-up factor: for high and low regularization respectively, our algorithm converges in and hours with subsampling factor , whereas the vanilla online algorithm requires about and hours. Qualitatively, Fig. 2 shows that with the same time budget, the proposed reduction approach with on half of the data gives better results than processing a small fraction of the data without reduction: segmented regions are less noisy and closer to processing the full data.
These results advocates the use of a subsampling rate of in this setting. When sparse matrix decomposition is part of a supervised pipeline with scoring capabilities, it is possible to find efficiently: start by setting it dereasonably high and decrease it geometrically until supervised performance (e.g. in classification) ceases to improve.
We validate the performance of the proposed algorithm on recommender systems for explicit feedback, a well-studied matrix completion problem. We evaluate the scalability of our method on datasets of different dimension: MovieLens 1M, MovieLens 10M, and 140M ratings Netflix dataset.
We compare our algorithm to a coordinate-descent based method (Yu et al., 2012), that provides state-of-the art convergence time performance on our largest dataset. Although stochastic gradient descent methods for matrix factorization can provide slightly better single-run performance (Takács et al., 2009), these are notoriously hard to tune and require a precise grid search to uncover a working schedule of learning rates. In contrast, coordinate descent methods do not require any hyper-parameter setting and are therefore more efficient in practice. We benchmarked various recommender-system codes (MyMediaLite, LibFM, SoftImpute, spira222https://github.com/mblondel/spira), and chose coordinate descent algorithm from spira as it was by far the fastest.
We stream user ratings to our algorithm: is the number of movies and is the number of users. As on Netflix dataset, this increases the benefit of using an online method. We have observed comparable prediction performance streaming item ratings. Past the first epoch, at iteration , every column of can be predicted by the last code that was computed from this column at iteration . At iteration , for all , . Prediction thus only requires an additional matrix computation after the factorization.
Successful prediction should take into account user and item biases. We compute these biases on train data following Hastie et al. (2014) (alternated debiasing). We use them to center the samples that are streamed to the algorithm, and to perform final prediction.
Both baseline and proposed algorithm are implemented in a computationally optimal way, enabling fair comparison based on CPU time. Benchmarks were run using a single 2.7 GHz Xeon CPU, with a 30 components dictionary. For Movielens datasets, we use a random of data for test and the rest for training. We average results on five train/test split for MovieLens in Table 1. On Netflix, the probe dataset is used for testing. Regularization parameter is set by cross-validation on the training set: the training data is split 3 times, keeping of Movielens datasets for evaluation and for Netflix, and grid search is performed on 15 values of between and . We assess the quality of obtained decomposition by measuring the root mean square error (RMSE) between prediction on the test set and ground truth. We use mini-batches of size .
|Dataset||Test RMSE||Convergence time||Speed|
We report the evolution of test RMSE along time in Fig. 4, along with its value at convergence and numerical convergence time in Table 1. Benchmarks are performed on the final run, after selection of parameter .
The two variants of the proposed method converge toward a solution that is at least as good as that of coordinate descent, and slightly better on Movielens 10M and Netflix. Our algorithm brings a substantial performance improvement on medium and large scale datasets. On Netflix, convergence is almost reached in minutes (score under deviation from final RMSE), which makes our method times faster than coordinate descent. Moreover, the relative performance of our algorithm increases with dataset size. Indeed, as datasets grow, less epochs are needed for our algorithm to reach convergence (Fig. 4). This is a significant advantage over coordinate descent, that requires a stable number of cycle on coordinates to reach convergence, regardless of dataset size. The algorithm with partial projection performs slightly better. This can be explained by the extra regularization on brought by this heuristic.
Unlike SGD, and similar to the vanilla online dictionary learning algorithm, our method does not critically suffer from hyper-parameter tuning. We tried weights as described in Sec. 2.3, and observed that a range of yields fast convergence. Theoretically, Mairal (2013) shows that stochastic majorization-minimization converges when . We verify this empirically, and obtain optimal convergence speed for . (Fig. 5). We report results for .
Whether it is sensor data, as fMRI, or e-commerce databases, sample sizes and number of features are rapidly growing, rendering current matrix factorization approaches intractable. We have introduced a online algorithm that leverages random feature subsampling, giving up to 8-fold speed and memory gains on large data. Datasets are getting bigger, and they often come with more redundancies. Such approaches blending online and randomized methods will yield even larger speed-ups on next-generation data.
The research leading to these results was supported by the ANR (MACARON project, ANR-14-CE23-0003-01 – NiConnect project, ANR-11-BINF-0004NiConnect) and has received funding from the European Union Seventh Framework Programme (FP7/2007-2013) under grant agreement no. 604102 (HBP).
Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition, pp. 2865–2872. IEEE, 2011.