In his well known paper [1], Achlioptas showed that Random Projections could be performed with a sparse projection matrix, whose nonzero entries are either $1$ or $-1$.

I have noticed that scikit-learn's implementation [2] takes advantage of the sparsity of the projection matrix to accelerate the projections by using scipy's sparse matrix multiplication. However, as described in the original paper, since the projection matrix only contains entries from $\{0, 1, -1\}$ the projection reduces to aggregate evaluation, i.e., summations and subtractions (but no multiplications). Therefore, computing the projection as a multiplication of sparse matrices is not taking full advantage of the simplicity of the projection matrix.

Formally, we want to compute the following matrix multiplication:

$X^\prime_{n\times k} = X_{n\times d}~R_{d\times k}$

where $X$ is an arbitrary data matrix and $R$ is the projection matrix, which is a sparse matrix where all non-zero elements are either $1$ or $-1$.

¿The question is, would it be possible to implement this with numpy in a more efficient manner, given the nature of the projection matrix?

[1] Achlioptas, D. (2001, May). Database-friendly random projections. In Proceedings of the twentieth ACM SIGMOD-SIGACT-SIGART symposium on Principles of database systems (pp. 274-281). ACM.

[2] https://github.com/scikit-learn/scikit-learn/blob/a24c8b46/sklearn/random_projection.py#L502


Yes, of course, technically you could. However, what would the gain be ? The projection itself (especially if the matrix is sparse) should be quite short compared to any other treatment you will apply (running a regression...).

Besides, these operations are often implemented at a very low level (probably calling sparse BLAS routines, themselves being highly optimized). A small investigation shows that the multiplication comes from : from .utils.extmath import safe_sparse_dot which is (with some dimension checks) np.dot(a, b), from numpy. Now, if you need some (potentially huge) speedups numpy can call ATLAS or BLAS routines (see per example http://markus-beuckelmann.de/blog/boosting-numpy-blas.html).

Rewriting routines at that level is possible but after a look at the git repo of OpenBLAS, this seems a Pyrrhic victory...

  • $\begingroup$ In the case that one of the matrices to be multiplied by safe_sparse_dot is sparse, scipy's sparse multiplication routines are called instead of numpy's dot. $\endgroup$ – Daniel López May 16 '18 at 7:18

In case this is useful for anyone, I found an easy way to implement the operation described in the question efficiently. In particular, I created a numpy ufunc with the help of numba:

from numba import guvectorize

@guvectorize(['void(float64[:], intp[:,:], intp[:,:], float64[:], float64[:])'], '(d),(y,l1),(y,l2),(k)->(k)',
             target='parallel', nopython=True, cache=True)

def databaseFriendlyRandomProjection(X, Rones, Rmones, zeros, out):

        for j in range(Rones.shape[1]):
            out[Rones[1,j]] += X[Rones[0,j]]

        for j in range(Rmones.shape[1]):
            out[Rmones[1,j]] -= X[Rmones[0,j]]

When the function is called, $X$ is the $n\times d$ data matrix, $Rones$ is a $2\times ?$ matrix which stores the indexes of all $1$s in the sparse projection matrix, $Rmones$ is a $2\times ?$ matrix which stores the indexes of all $-1$s in the sparse projection matrix, and $zeros$, $out$ are references to the $n\times k$ output matrix which needs to be initialized whith zeros by the user of the function.

Note that the function is designed to process a single data sample, and the broadcasting mechanism in numpy enables it to run in parallel for multiple samples.

In my experiments, this runs a bit faster than scipy's sparse matrix multiplication, especially if the sparsity of the projection matrix is bellow 50%. It also runs notably faster than numpy's dot (given that the projection matrix be sparse).

I'm sure this can be further optimized by someone with more experience with numba though, so I am open to suggestions.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.