Repository: scikit-learn/scikit-learn
Pull request: ENH Add inverse_transform to random projection transformers
#### Reference Issues/PRs

Fixes #21687

#### What does this implement/fix? Explain your changes.

Adds a `fit_inverse_transform` parameter to all transformers in  the `sklearn.random_projection` module: `GaussianRandomProjection` and `SparseRandomProjection`. When set to `True`, the pseudo-inverse of the components is computed during `fit()` and stored in `components_pinv_`, and `inverse_transform()` becomes available.

#### Any other comments?

Using the pseudo-inverse makes sense to me, and it seems to work fine, in the sense that `rnd_proj.transform(rnd_proj.inverse_transform(X))` equals `X`. However, this implementation uses `scipy.linalg.pinv()`, which scales very poorly to large datasets, which is a major use case for random projections. Perhaps it would make sense to use another approach if we detect that the `components_` array is large?

And perhaps there's a mathematical way to generate both a random matrix and its inverse more efficiently?

For the `SparseRandomProjection` transformer, computing the pseudo-inverse breaks sparsity. Perhaps there's a way to generate a sparse matrix that is "close enough" rather than using the pseudo-inverse?

In short: if there's a Random Projection expert in the room, please speak up!

That said, it seems to work fine now, so performance improvements could be pushed in follow-up PRs.

Problem statement:
Add inverse_transform() to random projection classes
### Describe the workflow you want to enable

`GaussianRandomProjection` and `SparseRandomProjection` should have an `inverse_transform()` method.

### Describe your proposed solution

Simply add the `inverse_transform()` method, with the appropriate input validation and tests, and return `X @ np.linalg.pinv(self.components_.T)`.

### Describe alternatives you've considered, if relevant

N/A

### Additional context

My understanding is that every transformer that performs a reversible transformation should have an `inverse_transform()` method, unless there's a really good reason not to (e.g., if it's terribly inefficient, or there's a better alternative). Was the `inverse_transform()` left out because it's too inefficient?

Hints from issue discussion:
As computing the inverse is costly, it probably should be done in the `fit()` method, instead of every time we call `inverse_transform()`. But since it's not needed that often, I guess it should only be computed if an `inverse` hyperparameter is set to `True` (and the default would be `False`).
I guess `np.linalg.pinv(self.components_.T)` can be cached/stored the first time `inverse_transform` is called.
I suppose it could be, but I thought the usual approach was to learn stuff in `fit()`, similar to how `probability=True` in SVC makes it learn probas and adds the `predict_proba()` method. WDYT?
true, you could open a PR with that approach.

Commit message: Add missing newline in test

Patch to review:
diff --git a/doc/whats_new/v1.1.rst b/doc/whats_new/v1.1.rst
index cc607aedb3946..84de13fc09f4d 100644
--- a/doc/whats_new/v1.1.rst
+++ b/doc/whats_new/v1.1.rst
@@ -231,6 +231,14 @@ Changelog
   :class:`~sklearn.random_projection.SparseRandomProjection`. :pr:`21330` by
   :user:`Loïc Estève <lesteve>`.
 
+- |Enhancement| Adds a `fit_inverse_transform` parameter to all transformers in
+  the :mod:`~sklearn.random_projection` module:
+  :class:`~sklearn.random_projection.GaussianRandomProjection` and
+  :class:`~sklearn.random_projection.SparseRandomProjection`. When set to True,
+  the pseudo-inverse of the components is computed during `fit` and stored in
+  `components_pinv_`, and `inverse_transform` becomes available.
+  :pr:`21701` by `Aurélien Geron <ageron>`.
+
 Code and Documentation Contributors
 -----------------------------------
 
diff --git a/sklearn/random_projection.py b/sklearn/random_projection.py
index 79ebe46042379..2484398fec6ae 100644
--- a/sklearn/random_projection.py
+++ b/sklearn/random_projection.py
@@ -31,6 +31,7 @@
 from abc import ABCMeta, abstractmethod
 
 import numpy as np
+from scipy import linalg
 import scipy.sparse as sp
 
 from .base import BaseEstimator, TransformerMixin
@@ -38,11 +39,11 @@
 
 from .utils import check_random_state
 from .utils.extmath import safe_sparse_dot
+from .utils.metaestimators import available_if
 from .utils.random import sample_without_replacement
-from .utils.validation import check_is_fitted
+from .utils.validation import check_array, check_is_fitted
 from .exceptions import DataDimensionalityWarning
 
-
 __all__ = [
     "SparseRandomProjection",
     "GaussianRandomProjection",
@@ -302,11 +303,18 @@ class BaseRandomProjection(
 
     @abstractmethod
     def __init__(
-        self, n_components="auto", *, eps=0.1, dense_output=False, random_state=None
+        self,
+        n_components="auto",
+        *,
+        eps=0.1,
+        dense_output=False,
+        fit_inverse_transform=False,
+        random_state=None,
     ):
         self.n_components = n_components
         self.eps = eps
         self.dense_output = dense_output
+        self.fit_inverse_transform = fit_inverse_transform
         self.random_state = random_state
 
     @abstractmethod
@@ -395,6 +403,12 @@ def fit(self, X, y=None):
             " not the proper shape."
         )
 
+        if self.fit_inverse_transform:
+            components = self.components_
+            if sp.issparse(components):
+                components = components.toarray()
+            self.components_pinv_ = linalg.pinv(components, check_finite=False)
+
         return self
 
     def transform(self, X):
@@ -431,6 +445,26 @@ def _n_features_out(self):
         """
         return self.n_components
 
+    @available_if(lambda self: self.fit_inverse_transform)
+    def inverse_transform(self, X):
+        """Project data back to its original space.
+
+        Returns an array X_original whose transform would be X.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix}, shape (n_samples, n_components)
+            Data to be transformed back.
+
+        Returns
+        -------
+        X_original : {ndarray, sparse matrix} of shape (n_samples, n_features)
+            Transformed data.
+        """
+        check_is_fitted(self)
+        X = check_array(X)
+        return X @ self.components_pinv_.T
+
 
 class GaussianRandomProjection(BaseRandomProjection):
     """Reduce dimensionality through Gaussian random projection.
@@ -463,6 +497,11 @@ class GaussianRandomProjection(BaseRandomProjection):
         Smaller values lead to better embedding and higher number of
         dimensions (n_components) in the target projection space.
 
+    fit_inverse_transform : bool, default=False
+        Learn the inverse transform by computing the pseudo-inverse of the
+        components during fit. Note that computing the pseudo-inverse does not
+        scale well to large matrices.
+
     random_state : int, RandomState instance or None, default=None
         Controls the pseudo random number generator used to generate the
         projection matrix at fit time.
@@ -505,11 +544,19 @@ class GaussianRandomProjection(BaseRandomProjection):
     (25, 2759)
     """
 
-    def __init__(self, n_components="auto", *, eps=0.1, random_state=None):
+    def __init__(
+        self,
+        n_components="auto",
+        *,
+        eps=0.1,
+        fit_inverse_transform=False,
+        random_state=None,
+    ):
         super().__init__(
             n_components=n_components,
             eps=eps,
             dense_output=True,
+            fit_inverse_transform=fit_inverse_transform,
             random_state=random_state,
         )
 
@@ -599,6 +646,13 @@ class SparseRandomProjection(BaseRandomProjection):
         If False, the projected data uses a sparse representation if
         the input is sparse.
 
+    fit_inverse_transform : bool, default=False
+        Learn the inverse transform by computing the pseudo-inverse of the
+        components during fit. Note that the components will be converted to a
+        dense array, and the resulting pseudo-inverse is also a dense array.
+        Moreover, computing the pseudo-inverse does not scale well to large
+        matrices.
+
     random_state : int, RandomState instance or None, default=None
         Controls the pseudo random number generator used to generate the
         projection matrix at fit time.
@@ -665,12 +719,14 @@ def __init__(
         density="auto",
         eps=0.1,
         dense_output=False,
+        fit_inverse_transform=False,
         random_state=None,
     ):
         super().__init__(
             n_components=n_components,
             eps=eps,
             dense_output=dense_output,
+            fit_inverse_transform=fit_inverse_transform,
             random_state=random_state,
         )
 
diff --git a/sklearn/tests/test_random_projection.py b/sklearn/tests/test_random_projection.py
index 1e894d906a3ad..fdb777853b559 100644
--- a/sklearn/tests/test_random_projection.py
+++ b/sklearn/tests/test_random_projection.py
@@ -373,3 +373,31 @@ def test_random_projection_feature_names_out(random_projection_cls):
     )
 
     assert_array_equal(names_out, expected_names_out)
+
+
+###############################################################################
+# tests on inverse transform
+###############################################################################
+
+
+@pytest.mark.parametrize("random_projection_cls", all_RandomProjection)
+def test_inverse_transform(random_projection_cls):
+    random_projection = random_projection_cls(
+        n_components=2, fit_inverse_transform=True
+    )
+    projected = random_projection.fit_transform(data)
+    projected_back = random_projection.inverse_transform(projected)
+    assert projected_back.shape == data.shape
+    projected_again = random_projection.transform(projected_back)
+    assert_array_almost_equal(projected, projected_again)
+    assert random_projection.components_pinv_.shape == (n_features, 2)
+
+
+@pytest.mark.parametrize("random_projection_cls", all_RandomProjection)
+def test_no_inverse_transform_by_default(random_projection_cls):
+    random_projection = random_projection_cls(n_components=2)
+    projected = random_projection.fit_transform(data)
+    with pytest.raises(AttributeError):
+        random_projection.inverse_transform(projected)
+    with pytest.raises(AttributeError):
+        random_projection.components_pinv_


Code review comments:
diff_hunk: @@ -395,6 +403,12 @@ def fit(self, X, y=None):
             " not the proper shape."
         )
 
+        if self.fit_inverse_transform:
+            components = self.components_
+            if sp.issparse(components):
+                components = components.toarray()
path: sklearn/random_projection.py
text: @user1:
Doing this silently seems dangerous from a memory usage point of view. I would rather use a local pinv implementation inspired from scipy's but based on `scipy.sparse.linalg.svds`.

This could even be contributed upstream as `scipy.sparse.linalg.pinv` or even by upgrading `scipy.linalg.pinv` to accept sparse arrays once https://github.com/scipy/scipy/pull/14822 is merged.

@user1:
I you don't feel like doing this as par of this PR, we do a first PR to add `fit_inverse_transform` to `GaussianRandomProjection` only for now and deal with `SparseRandomProjection` in a follow-up PR.

@author:
I understand. I think I'll go for the second option, and just add `fit_inverse_transform` to `GaussianRandomProjection` for now.

@author:
Actually, scratch that, I think I'll have a go at using a local `pinv` based on `scipy.sparse.linalg.svds`, as you suggest.

diff_hunk: @@ -431,6 +445,26 @@ def _n_features_out(self):
         """
         return self.n_components
 
+    @available_if(lambda self: self.fit_inverse_transform)
+    def inverse_transform(self, X):
+        """Project data back to its original space.
+
+        Returns an array X_original whose transform would be X.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix}, shape (n_samples, n_components)
+            Data to be transformed back.
+
+        Returns
+        -------
+        X_original : {ndarray, sparse matrix} of shape (n_samples, n_features)
+            Transformed data.
+        """
+        check_is_fitted(self)
+        X = check_array(X)
path: sklearn/random_projection.py
text: @user1:
We should probably accept sparse data here.

diff_hunk: @@ -373,3 +373,31 @@ def test_random_projection_feature_names_out(random_projection_cls):
     )
 
     assert_array_equal(names_out, expected_names_out)
+
+
+###############################################################################
+# tests on inverse transform
+###############################################################################
+
+
+@pytest.mark.parametrize("random_projection_cls", all_RandomProjection)
+def test_inverse_transform(random_projection_cls):
path: sklearn/tests/test_random_projection.py
text: @user2:
Please parametrize this test to also run on sparse data:

```python
@user1.mark.parametrize("data", [data, data_csr])
```

diff_hunk: @@ -373,3 +373,31 @@ def test_random_projection_feature_names_out(random_projection_cls):
     )
 
     assert_array_equal(names_out, expected_names_out)
+
+
+###############################################################################
+# tests on inverse transform
+###############################################################################
+
+
+@pytest.mark.parametrize("random_projection_cls", all_RandomProjection)
+def test_inverse_transform(random_projection_cls):
+    random_projection = random_projection_cls(
+        n_components=2, fit_inverse_transform=True
+    )
+    projected = random_projection.fit_transform(data)
+    projected_back = random_projection.inverse_transform(projected)
+    assert projected_back.shape == data.shape
+    projected_again = random_projection.transform(projected_back)
+    assert_array_almost_equal(projected, projected_again)
+    assert random_projection.components_pinv_.shape == (n_features, 2)
path: sklearn/tests/test_random_projection.py
text: @user1:
```suggestion
    assert random_projection.components_pinv_.shape == (n_features, 2)
    assert_all_close(
        random_projection.components_ @ random_projection.components_pinv_,
        np.eye(random_projection.n_components)
    )
```

diff_hunk: @@ -395,6 +403,12 @@ def fit(self, X, y=None):
             " not the proper shape."
         )
 
+        if self.fit_inverse_transform:
+            components = self.components_
+            if sp.issparse(components):
+                components = components.toarray()
+            self.components_pinv_ = linalg.pinv(components, check_finite=False)
path: sklearn/random_projection.py
text: @user1:
Is it useful to expose publicly the pseudo-inverse or can we make it kinda of private (i.e. `self._components_pinv`

@user1:
You mentioned that `pinv` does not scale with large matrices. I have 2 questions:

- Should we use `pinv` or `pinv2` depending on the shape of the `X`?
- Is there a way to approximate the inverse? I am recalling the following [PR](https://github.com/scikit-learn/scikit-learn/pull/19971) where we wanted to do so with an approximation in Nystroem (I did not look at the PR thought). I don't know if we could have a trick here to have such an approximation here?

@user2:
I think it's useful to expose publicly. Maybe we could change the name to `inverse_components_` though to make the name less dependent on an implementation detail.

@user2:
> Should we use pinv or pinv2 depending on the shape of the X?

`pinv2` is deprecated.

diff_hunk: @@ -431,6 +445,26 @@ def _n_features_out(self):
         """
         return self.n_components
 
+    @available_if(lambda self: self.fit_inverse_transform)
+    def inverse_transform(self, X):
+        """Project data back to its original space.
+
+        Returns an array X_original whose transform would be X.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix}, shape (n_samples, n_components)
path: sklearn/random_projection.py
text: @user1:
```suggestion
        X : {array-like, sparse matrix} of shape (n_samples, n_components)
```

diff_hunk: @@ -373,3 +373,31 @@ def test_random_projection_feature_names_out(random_projection_cls):
     )
 
     assert_array_equal(names_out, expected_names_out)
+
+
+###############################################################################
+# tests on inverse transform
+###############################################################################
path: sklearn/tests/test_random_projection.py
text: @user1:
We would not put this comment.