The test failure in sklearn/tests/test_calibration.py looks unrelated to this PR.

Acccording to this SO answer:

numpy.linalg.pinv() approximates the Moore-Penrose pseudo inverse using an SVD (the LAPACK method dgesdd to be precise), whereas scipy.linalg.pinv() solves a model linear system in the least squares sense to approximate the pseudo inverse (using dgelss).

They perform differently, both in terms of speed and memory usage, depending on the size of the input. So perhaps a follow-up PR should give the user the option to select the algorithm to use?

This looks pretty good an clean, thanks @ageron . The only thing I'd add is a section to the relevant user guide (the relevant .rst) to explain the inverse transform, and how it's calculated, and how it's enabled.

numpy.linalg.pinv() approximates the Moore-Penrose pseudo inverse using an SVD (the LAPACK method dgesdd to be precise), whereas scipy.linalg.pinv() solves a model linear system in the least squares sense to approximate the pseudo inverse (using dgelss).

Are you sure this is still valid with recent version of scipy? https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.pinv.html

It's not possible to select the LAPACK driver though.

Thanks @adrinjalali , @glemaitre , and @ogrisel for taking such a thorough look at this PR, and for the interesting feedback.

@ogrisel , I tested both np.linalg.pinv() and scipy.linalg.pinv() on an array of shape [5_000, 20_000], and SciPy took 25.5 minutes while NumPy took only 11.6 minutes (on Colab). So NumPy's implementation was twice as fast! SciPy uses its decomp_svd() function, while NumPy uses its svd() function. I haven't looked further yet, so I'm still not sure why NumPy's implementation is twice as fast as SciPy's.

I'll work on the pinv() implementation for sparse arrays soon, but in the meantime I made the rest of the changes you suggested, including:

• Added a section in the User Guide in random_projection.rst
• Renamed components_pinv_ to inverse_components_. I didn't like the name either!
• Added a note to say that even if X is sparse, X_original is dense, which may use a lot of RAM.
• Fixed: X : {array-like, sparse matrix} of shape (n_samples, n_components)
• fit_transform() now accepts both dense and sparse arrays
• Added testing for both dense and sparse inputs in test_random_projection.py
• Removed the ######### header in test_random_projection.py
• Added the test that random_projection.components_ @ random_projection.inverse_components_ is the identity matrix

Ok, SparseRandomProjection.fit() now computes the pseudo-inverse (when fit_inverse_transform is set to True) without converting the components_ to a dense array.

The pseudo-inverse implementation is based on scipy.sparse.linalg.svds(), as you suggested @ogrisel . I had to call it twice, though: once for the first half of the components, and once for the second half. That's because scipy.sparse.linalg.svds() does not support getting all the components at once.

The rest of the pseudo-inverse implementation was taken in large part from scipy.linalg.pinv(), but I simplified it a lot, removing every parameter except the sparse matrix. I just wanted to keep things as basic as possible, since this is meant to be a temporary workaround until SciPy's pinv() supports sparse matrices.

Done

Yikes, something's weird: I'm getting no errors in my tests locally, but it's failing in CI. I'll have to investigate this. I'm thinking it might be an error that happens occasionally, if the first k//2 components returned by the first call to scipy.sparse.linalg.svds() are not compatible with the k - k//2 components returned by the second call. If that's the case, I don't really see a way to fix that. Any ideas? If there's no solution, I'll just remove inverse_transform() for the sparse case. All this hard work for nothing! 😭

Edit: I've run some tests and I can confirm that my implementation of svd() based on calling svds() twice only worked about 90% of the time, it depended on the dataset. Since the random seed was fixed, I was always testing on the same dataset, so I never ran into this issue. I wrote a new implementation that fixes this issue, see below.

Okay, I've rewritten _svd_for_sparse_matrix(), and everything seems to work fine now. Instead of calling svds() twice, the code just calls it once. To workaround the fact that svds can only return min(a.shape) - 1 components (missing one), the code now adds an extra row or column (or both) full of zeros and crops U or Vt appropriately.

Note that I removed the test assert_array_almost_equal(random_projection.components_ @ random_projection.inverse_components_, np.eye(random_projection.n_components)). Indeed, a matrix multiplied by its pseudo-inverse is not always equal to the identity matrix, for example:

>>> import numpy as np
>>> from scipy.linalg import pinv
>>> a = np.array([[1, 2, 0], [2, 4, 0]])
>>> a @ pinv(a)
array([[0.2, 0.4],
       [0.4, 0.8]])

However, I have added a test that inverse_components_ is equal to scipy.linalg.pinv(components_), or equal to scipy.linalg.pinv(components_.toarray()) if components_ is sparse.

Also, the test now runs many times with different random states and different shapes: with n_rows < n_cols, or n_rows > n_cols, or n_rows == n_cols. I had to silence an unrelated warning in the case where n_components > n_cols.

Thanks @ageron !

Since the inverted components are always dense, why do we bother to implement a pinv for sparse ? Densify the array first and call numpy's pinv could (maybe ?) be more efficient. I understand that the goal was to mitigate the memory usage but we are creating the dense inverse components array anyway.

My pleasure @jeremiedbb !

I think there were two reasons: (1) as you pointed out, it saves saves memory usage (probably halves it), and (2) someone (perhaps @ogrisel ? I can't recall) told me that it may be useful in other places until scipy offers the same functionality, which I believe is on their roadmap.

I think there were two reasons: (1) as you pointed out, it saves saves memory usage (probably halves it), and (2) someone (perhaps @ogrisel ? I can't recall) told me that it may be useful in other places until scipy offers the same functionality, which I believe is on their roadmap.

Yes, the worst case is temporarily having 2 dense arrays like components instead of 1. After an irl discussion with @ogrisel, we agreed that it's acceptable if it allows to avoid the burden of maintaining our own version of a pseudo inverse for sparse matrices.

Tbh I don't have much time right now, I'm trying to finish writing the 3rd edition of my book, so if a kind soul wants to investigate this point, that would be great.

If you don't have much time I can directly push these changes if you want.

Hi @jeremiedbb ,
Thanks for reviewing, and for offering to push the last changes, I really appreciate it. I understand the trade-off about the sparse pinv. There was a bit of a gambler's fallacy on my part: I spent so much time working on the sparse pinv that I was relunctant to just drop it! But not worries, please feel free to delete it.

I spent so much time working on the sparse pinv that I was relunctant to just drop it!

I understand the feeling :)

Thanks @ageron !

Thanks @ogrisel, @jeremiedbb and @glemaitre for the thorough review, I'm impressed by how pro and helpful you all are, it's no wonder Scikit-Learn is so good! 👍