Distributed PCA or an equivalent

Question

We normally have fairly large datasets to model on, just to give you an idea:

• over 1M features (sparse, average population of features is around 12%);
• over 60M rows.

A lot of modeling algorithms and tools don't scale to such wide datasets.

So we're looking for a dimensionality reduction implementation that runs distributely (i.e. in Spark/Hadoop/ etc). We think to bring number of features down to several thousand.

Since PCA operate on matrix multiplication, which don't distribute very well over a cluster of servers, we're looking at other algorithms or probably, at other implementations of distributed dimensionality reduction.

Anyone ran into similar issues? What do you do to solve this?

There is a Cornell/Stanford Abstract on "Generalized Low-Rank Models" http://web.stanford.edu/~boyd/papers/pdf/glrm.pdf that talks specifically into this:

• page 8 "Parallelizing alternating minimization" tells how it can be distributed;
• also page 9 "Missing data and matrix completion" talks how sparse/ missing data can be handled.

GLRM although seems to be what we are looking for, but we can't find good actual implementations of those ideas.

• Somebody wrote a version for Spark (e.g. https://github.com/rezazadeh/spark/tree/glrm/examples/src/main/scala/org/apache/spark/examples/glrm but it looks more of a proof of a concept rather than a version that we could use in production);
• H2O has a version of GLRM http://docs.h2o.ai/h2o/latest-stable/h2o-docs/data-science/glrm.html but it actually doesn't scale to our dataset size (see above).

Update 7/15/2018: Another Abstract is Fast Randomized SVD from Facebook (read here http://tygert.com/spark.pdf ) and also idea to do low-rank matrix approximation using ALS - http://tygert.com/als.pdf . Although there is no clear way how to use them now - see discussion at https://github.com/facebook/fbpca/issues/6

Any other ideas how to tackle this? Other available GLRM or other distributed dimensionality reduction implementations?

Answer 1

A note about PCA

I'd like to clarify the language I'm going to use. Principal Component Analysis is an extremely wide umbrella term that is used for algorithms that try to capture "the gist" of the data (usually expressed as a matrix). There are two components in those algorithms: what are they trying to capture exactly and model that they employ to capture this information (the general idea is that the compressed version of the data should occupy less space and maybe reveal some patterns in the data).

Given these definitions, SVD is an optimal way to capture L2 norm of the matrix (meaning, we are minimizing $||A-\hat A||$) given low-rank linear model $\hat A = USV^T$ (U and V are matrices with orthonormal columns and S is diagonal). SVD++ or FunkSVD are not exactly the same, they have different loss function (they either do not take zeros into account or subsample them), they have an explicit regularization and they have no guarantees about orthogonality of found singular vectors (which you might not care much about).

Best distributed SVD libs

If you think that capturing the L2 norm of your matrix is a good idea, then I can suggest using one of two libraries that implement Randomized SVD algorithm: SparkRSVD library in Scala for Spark or Dask ML Truncated SVD in Python for, well, Dask. We've written the first and tested it against the second on matrices up to 100 mln x 100 mln rows and columns with very small sparsity (see our blog post here). They are both I think good options, the choice depends mainly on the language/framework you're most comfortable with. Also, we found Dask to be faster on relatively small matrices, but the Spark lib can scale to larger matrices (from what we tested).

Now, why Randomized SVD and not SVD from MLLib, for instance?

If you look under the hood, Spark MLLib implements SVD only for RowMatrix (with the limitation on num_cols <= MAX_INTEGER which is fine) and computes a Gramian matrix ($G=A^T A$) to produce a smaller matrix (in your case it should become almost dense 1 mln x 1 mln) that is easy to work with which should either fit in memory of a single node or at least num_rows * num_factors should. Which might or might not be okay for you depending on what kind of hardware you have or how many latent factors you need to reconstruct.

In contrast, Randomized SVD algorithms are based on the ideas of randomized algorithms and efficient QR decomposition. They do not require candom access to the matrix A, just the ability to multiply by it, and do not require storing the result of those multiplications at the same node, just a $\text{num_factors} \times \text{num_factors}$ matrix.

P.S. Updated the answer to be more precise about what Spark MLLib implements.