Repository: scikit-learn/scikit-learn
Pull request: FIX KernelPCA inverse transform when gamma is not given
#### Reference Issues/PRs

Fixes #26280.

#### What does this implement/fix? Explain your changes.

Quote @jeremiedbb:
> Originally when `gamma` is `None`, its set to `1 / n_features` each time the kernel is called, not once at the beginning of fit. It means that when the kernel is called at `inverse_transform`, the number of features is different and hence `gamma` is different.

This PR intends to set a private attribute `_gamma` at the first `fit` (or `fit_transform`) and use `_gamma` instead of `gamma` in further transforms.

#### Any other comments?

Not sure if this approach is neat enough, so if maintainers want an alternative approach, please let me know!

Problem statement:
KernelPCA inverse transform behaves unexpectly.
### Describe the bug

Hi!
we wanted to use inverse transform in kernel PCA. There is a parameter gamma, which is defined as 1/num_features if gamma=None. However if gamma is provided as gamma=1/num_features, the resultant inverse transform results is different. Why is that? How the parameter gamma is related to inverse transform?

### Steps/Code to Reproduce

```
from sklearn.datasets import load_iris
from sklearn.decomposition import KernelPCA

#load iris dataset (features only)
data = load_iris().data

#number of features in the dataset
n_features = data.shape[1]

#run KernelPCA with default value of gamma (1/n_features)
pca1 = KernelPCA(n_components = 2, random_state = 7, fit_inverse_transform = True, kernel = 'rbf')
pca1.fit(data)

#set the value of gamma to 1/n_features
pca2 = KernelPCA(n_components = 2, gamma = 1/n_features, random_state = 7, fit_inverse_transform = True, kernel = 'rbf')
pca2.fit(data)

#eigenvalues and eigenvectors are the same
print((pca1.eigenvalues_ == pca2.eigenvalues_).all())
print((pca1.eigenvectors_ == pca2.eigenvectors_).all())

#projections of components are the same
latent1 = pca1.transform(data)
latent2 = pca2.transform(data)
print((latent1 == latent2).all())

#but reconstructed data from these projections are different
reconstructed_data1 = pca1.inverse_transform(latent1)
reconstructed_data2 = pca2.inverse_transform(latent2)
print((reconstructed_data1 == reconstructed_data2).all())
```

### Expected Results

True
True
True
True

### Actual Results

True
True
True
False

### Versions

```shell
System:
    python: 3.10.6 (main, Mar 10 2023, 10:55:28) [GCC 11.3.0]
executable: /usr/bin/python3
   machine: Linux-5.15.0-70-generic-x86_64-with-glibc2.35

Python dependencies:
      sklearn: 1.2.2
          pip: 22.0.2
   setuptools: 59.6.0
        numpy: 1.21.5
        scipy: 1.8.0
       Cython: None
       pandas: 2.0.0
   matplotlib: 3.7.1
       joblib: 1.2.0
threadpoolctl: 3.1.0

Built with OpenMP: True

threadpoolctl info:
       user_api: openmp
   internal_api: openmp
         prefix: libgomp
       filepath: /home/petr/.local/lib/python3.10/site-packages/scikit_learn.libs/libgomp-a34b3233.so.1.0.0
        version: None
    num_threads: 8

       user_api: blas
   internal_api: openblas
         prefix: libopenblas
       filepath: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so
        version: 0.3.20
threading_layer: pthreads
   architecture: Haswell
    num_threads: 8
```

Hints from issue discussion:
From what I recall and what the documentation says, the `inverse_transform` is a learnt approximation: https://scikit-learn.org/stable/modules/generated/sklearn.decomposition.KernelPCA.html#sklearn.decomposition.KernelPCA.inverse_transform

It is not surprising to not recover the original data.
> It is not surprising to not recover the original data.

Yes but here it's not a comparison against the original data, it's a comparison between 2 run with the (supposed) same parameters.

Actually I think there's a bug here. 

- When `gamma` is None, its set to `1/n_feature` each time the kernel is called, not once at the beginning of fit. It means that when the kernel is called at `inverse_transform`, the number of features is different and hence gamma is different.
- When gamma is set, the same value is used at fit and at inverse_transform

Hi @jeremiedbb, what about adding a line
```python
self._gamma = 1 / X.shape[1] if self.gamma is None else self.gamma
```
right after validating `X` in the `fit` method, and then in `_get_kernel`, we use `self._gamma` rather than `self.gamma`? This indeed solves this specific issue and passes `pytest sklearn/decomposition/tests`. If you think this is a valid approach and is neat enough, I'm happy to make a PR for this.
Since there is no response, I've opened a PR as above.

Commit message: changelog added

Patch to review:
diff --git a/doc/whats_new/v1.3.rst b/doc/whats_new/v1.3.rst
index 5965cd90c1b89..903a70c9367d9 100644
--- a/doc/whats_new/v1.3.rst
+++ b/doc/whats_new/v1.3.rst
@@ -230,6 +230,11 @@ Changelog
   :class:`decomposition.MiniBatchNMF` which can produce different results than previous
   versions. :pr:`25438` by :user:`Yotam Avidar-Constantini <yotamcons>`.
 
+- |Fix| :class:`decomposition.KernelPCA` now uses the correct `gamma` if `gamma` is
+  `None` every time the kernel is called. It should be equal to `1/num_features` of
+  the data that it is initially fitted on.
+  :pr:`26337` by :user:`Yao Xiao <Charlie-XIAO>`.
+
 :mod:`sklearn.discriminant_analysis`
 ....................................
 
diff --git a/sklearn/decomposition/_kernel_pca.py b/sklearn/decomposition/_kernel_pca.py
index 72f20d99bd548..4e5e009a43f74 100644
--- a/sklearn/decomposition/_kernel_pca.py
+++ b/sklearn/decomposition/_kernel_pca.py
@@ -300,7 +300,7 @@ def _get_kernel(self, X, Y=None):
         if callable(self.kernel):
             params = self.kernel_params or {}
         else:
-            params = {"gamma": self.gamma, "degree": self.degree, "coef0": self.coef0}
+            params = {"gamma": self._gamma, "degree": self.degree, "coef0": self.coef0}
         return pairwise_kernels(
             X, Y, metric=self.kernel, filter_params=True, n_jobs=self.n_jobs, **params
         )
@@ -419,6 +419,7 @@ def fit(self, X, y=None):
         if self.fit_inverse_transform and self.kernel == "precomputed":
             raise ValueError("Cannot fit_inverse_transform with a precomputed kernel.")
         X = self._validate_data(X, accept_sparse="csr", copy=self.copy_X)
+        self._gamma = 1 / X.shape[1] if self.gamma is None else self.gamma
         self._centerer = KernelCenterer()
         K = self._get_kernel(X)
         self._fit_transform(K)
diff --git a/sklearn/decomposition/tests/test_kernel_pca.py b/sklearn/decomposition/tests/test_kernel_pca.py
index c6b41ab169810..12bea0b9ddd47 100644
--- a/sklearn/decomposition/tests/test_kernel_pca.py
+++ b/sklearn/decomposition/tests/test_kernel_pca.py
@@ -523,3 +523,28 @@ def test_kernel_pca_feature_names_out():
 
     names = kpca.get_feature_names_out()
     assert_array_equal([f"kernelpca{i}" for i in range(2)], names)
+
+
+def test_kernel_pca_inverse_correct_gamma():
+    """Check that gamma is set correctly when not provided.
+
+    Non-regression test for #26280
+    """
+    rng = np.random.RandomState(0)
+    X = rng.random_sample((5, 4))
+
+    kpca1 = KernelPCA(
+        n_components=2, random_state=rng, fit_inverse_transform=True, kernel="rbf"
+    ).fit(X)
+    kpca2 = KernelPCA(
+        n_components=2,
+        gamma=1 / X.shape[1],
+        random_state=rng,
+        fit_inverse_transform=True,
+        kernel="rbf",
+    ).fit(X)
+
+    reconstructed_X1 = kpca1.inverse_transform(kpca1.transform(X))
+    reconstructed_X2 = kpca2.inverse_transform(kpca1.transform(X))
+
+    assert (reconstructed_X1 == reconstructed_X2).all()


Code review comments:
diff_hunk: @@ -523,3 +523,28 @@ def test_kernel_pca_feature_names_out():
 
     names = kpca.get_feature_names_out()
     assert_array_equal([f"kernelpca{i}" for i in range(2)], names)
+
+
+def test_kernel_pca_inverse_correct_gamma():
+    """Check that gamma is set correctly when not provided.
+
+    Non-regression test for #26280
+    """
+    rng = np.random.RandomState(0)
+    X = rng.random_sample((5, 4))
+
+    kpca1 = KernelPCA(
+        n_components=2, random_state=rng, fit_inverse_transform=True, kernel="rbf"
+    ).fit(X)
path: sklearn/decomposition/tests/test_kernel_pca.py
text: @user1:
Just a small nitpick, but I think it would be nice to set `gamma=None` explicitly just to make it clearer at a glance what is being compared. Additionally you could also set the shared arguments as a dictionary that you unpack in both, just to minimize clutter.

diff_hunk: @@ -523,3 +523,28 @@ def test_kernel_pca_feature_names_out():
 
     names = kpca.get_feature_names_out()
     assert_array_equal([f"kernelpca{i}" for i in range(2)], names)
+
+
+def test_kernel_pca_inverse_correct_gamma():
+    """Check that gamma is set correctly when not provided.
+
+    Non-regression test for #26280
+    """
+    rng = np.random.RandomState(0)
+    X = rng.random_sample((5, 4))
+
+    kpca1 = KernelPCA(
+        n_components=2, random_state=rng, fit_inverse_transform=True, kernel="rbf"
+    ).fit(X)
+    kpca2 = KernelPCA(
+        n_components=2,
+        gamma=1 / X.shape[1],
+        random_state=rng,
+        fit_inverse_transform=True,
+        kernel="rbf",
+    ).fit(X)
+
+    reconstructed_X1 = kpca1.inverse_transform(kpca1.transform(X))
+    reconstructed_X2 = kpca2.inverse_transform(kpca1.transform(X))
+
+    assert (reconstructed_X1 == reconstructed_X2).all()
path: sklearn/decomposition/tests/test_kernel_pca.py
text: @user1:
The `assert_array_equal` is provided by `numpy` and provides a more detailed error message along with being a bit more robust to certain edge-cases present in `numpy` arrays.

```suggestion
    X1_recon = kpca1.inverse_transform(kpca1.transform(X))
    X2_recon = kpca2.inverse_transform(kpca1.transform(X))

    assert_array_equal(X1_recon, X2_recon)
```