scikit-learn · kstoneriv3 · Apr 24, 2021 · Apr 24, 2021 · Apr 24, 2021 · Apr 24, 2021
diff --git a/doc/whats_new/v1.0.rst b/doc/whats_new/v1.0.rst
@@ -219,6 +219,13 @@ Changelog
  :func:`~sklearn.inspection.permutation_importance`.
  :pr:`19411` by :user:`Simona Maggio <simonamaggio>`.

+:mod:`sklearn.kernel_approximation`
+...................................
+
+- |Enhancement| Add
+  :meth:`~sklearn.kernel_approximation.Nystroem.inverse_transform`.
+  :pr:`19971` by :user:`Kei Ishikawa <kstoneriv3>`.
+
 :mod:`sklearn.linear_model`
 ...........................


diff --git a/sklearn/kernel_approximation.py b/sklearn/kernel_approximation.py
@@ -12,6 +12,7 @@

 import numpy as np
 import scipy.sparse as sp
+from scipy import linalg
 from scipy.linalg import svd
 try:
    from scipy.fft import fft, ifft
@@ -20,6 +21,7 @@

 from .base import BaseEstimator
 from .base import TransformerMixin
+from .exceptions import NotFittedError
 from .utils import check_random_state, as_float_array
 from .utils.extmath import safe_sparse_dot
 from .utils.validation import check_is_fitted
@@ -683,6 +685,9 @@ class Nystroem(TransformerMixin, BaseEstimator):

    Attributes
    ----------
+    basis_indices_: ndarray of shape (n_basis)
+        Indices of ``basis`` in the training set.
+
    components_ : ndarray of shape (n_components, n_features)
        Subset of training points used to construct the feature map.

@@ -730,18 +735,61 @@ class Nystroem(TransformerMixin, BaseEstimator):
    """
    @_deprecate_positional_args
    def __init__(self, kernel="rbf", *, gamma=None, coef0=None, degree=None,
-                 kernel_params=None, n_components=100, random_state=None,
-                 n_jobs=None):
+                 kernel_params=None, n_components=100, alpha=1.0,
+                 fit_inverse_transform=False, random_state=None, n_jobs=None):

        self.kernel = kernel
        self.gamma = gamma
        self.coef0 = coef0
        self.degree = degree
        self.kernel_params = kernel_params
        self.n_components = n_components
+        self.alpha = alpha
+        self.fit_inverse_transform = fit_inverse_transform
        self.random_state = random_state
        self.n_jobs = n_jobs

+    def _fit_inverse_transform(self, X_transformed, X):
+        if hasattr(X, "tocsr"):
+            raise NotImplementedError("Inverse transform not implemented for "
+                                      "sparse matrices!")
+
+        # Let I and T be the indices of inducing points and training points.
+        Kii = pairwise_kernels(X_transformed[self.basis_indices_],  # K_{I, I}
+                               metric=self.kernel,
+                               filter_params=True,
+                               n_jobs=self.n_jobs,
+                               **self._get_kernel_params())
+        Kit = pairwise_kernels(X_transformed[self.basis_indices_],  # K_{I, T}
+                               X_transformed,
+                               metric=self.kernel,
+                               filter_params=True,
+                               n_jobs=self.n_jobs,
+                               **self._get_kernel_params())
+
+        # By Woodbury matrix identity, following formula(typed in latex) holds:
+        #
+        # \tilde{K}_{new, T} (\tilde{K}_{T, T} + \alpha I)^{-1} y
+        # = K_{new, I} K_{I, I}^{-1}K_{I, T}
+        # \left( \frac{1}{\alpha} I -  \frac{1}{\alpha}K_{T, I}(\alpha K_{I, I}
+        # + K_{I, T}K_{T, I})^{-1}K_{I, T} \right )y
+        # = K_{new, I} \left( \alpha \cdot K_{I, I}\right)^{-1}
+        # \left\{K_{I, T}y - K_{I, T} K_{T, I}(\alpha K_{I, I}
+        # + K_{I, T}K_{T, I})^{-1}K_{I, T}y \right\}
+        #
+        # where \tilde{K} implies the kernel matrix approximated by
+        # Nystrom approximation.
+
+        A = Kit @ X
+        B = Kit @ Kit.T
+        # TODO(kstoneriv3): The following line often results in warning like
+        # "Ill-conditioned matrix (rcond=1.80677e-18)". So we use pinv instead.
+        # C = linalg.solve(self.alpha * Kii + B, A, sym_pos=True)
+        C = linalg.pinvh(self.alpha * Kii + B) @ A
+        D = A - B @ C
+        self.dual_coef_ = linalg.pinvh(self.alpha * Kii) @ D
+        self.X_transformed_fit_ = X_transformed[self.basis_indices_]
+
    def fit(self, X, y=None):
        """Fit estimator to data.

@@ -778,11 +826,17 @@ def fit(self, X, y=None):
                                        **self._get_kernel_params())

        # sqrt of kernel matrix on basis vectors
-        U, S, V = svd(basis_kernel)
+        U, S, V = svd(basis_kernel)  # TODO(kstoneriv3): Why not linalg.eigh()?
        S = np.maximum(S, 1e-12)
        self.normalization_ = np.dot(U / np.sqrt(S), V)
        self.components_ = basis
-        self.component_indices_ = inds
+        self.component_indices_ = basis_inds
+        self.basis_indices_ = basis_inds
+
+        if self.fit_inverse_transform:
+            X_transformed = self.transform(X)
+            self._fit_inverse_transform(X_transformed, X)
+
        return self

    def transform(self, X):
@@ -812,6 +866,38 @@ def transform(self, X):
                                    **kernel_params)
        return np.dot(embedded, self.normalization_.T)

+    def inverse_transform(self, X):
+        """Transform X back to original space.
+
+        This method approximates the inverse transformation using
+        a learned pre-image. The pre-image is learned by kernel ridge
+        regression of the original data on their low-dimensional representation
+        vectors. For the efficiency of kernel ridge regression, the kernel for
+        regression is a Nystrom approximation of the original kernel.
+
+        Parameters
+        ----------
+        X : {array-like, sparse matrix} of shape (n_samples, n_components)
+
+        Returns
+        -------
+        X_new : ndarray of shape (n_samples, n_features)
+
+        References
+        ----------
+        "Learning to Find Pre-Images", G BakIr et al, 2004.
+        """
+        if not self.fit_inverse_transform:
+            raise NotFittedError("The fit_inverse_transform parameter was not"
+                                 " set to True when instantiating and hence "
+                                 "the inverse transform is not available.")
+        K = pairwise_kernels(X, self.X_transformed_fit_,
+                             metric=self.kernel,
+                             filter_params=True,
+                             n_jobs=self.n_jobs,
+                             **self._get_kernel_params())
+        return np.dot(K, self.dual_coef_)
+
    def _get_kernel_params(self):
        params = self.kernel_params
        if params is None:

diff --git a/sklearn/tests/test_kernel_approximation.py b/sklearn/tests/test_kernel_approximation.py
@@ -7,12 +7,12 @@
 from sklearn.utils._testing import assert_array_equal
 from sklearn.utils._testing import assert_array_almost_equal

-from sklearn.metrics.pairwise import kernel_metrics
 from sklearn.kernel_approximation import RBFSampler
 from sklearn.kernel_approximation import AdditiveChi2Sampler
 from sklearn.kernel_approximation import SkewedChi2Sampler
 from sklearn.kernel_approximation import Nystroem
 from sklearn.kernel_approximation import PolynomialCountSketch
+from sklearn.metrics.pairwise import kernel_metrics
 from sklearn.metrics.pairwise import polynomial_kernel, rbf_kernel, chi2_kernel

 # generate data
@@ -332,3 +332,21 @@ def test_nystroem_precomputed_kernel():
                      **param)
        with pytest.raises(ValueError, match=msg):
            ny.fit(K)
+
+
+def test_nystroem_inverse_transform_reconstruction():
+    # Test if the reconstruction is a good approximation.
+    n = 100
+    d = 2
+    rng = np.random.RandomState(0)
+    X = rng.randn(n * d).reshape(n, d)
+    Y = rng.randn(n * d).reshape(n, d)
+
+    nystroem = Nystroem(
+        n_components=40, kernel='rbf', fit_inverse_transform=True,
+        alpha=3e-3, gamma=4e-1, random_state=0
+    )
+    nystroem.fit(X)
+    Y_trans = nystroem.fit_transform(Y)
+    Y_reconst = nystroem.inverse_transform(Y_trans)
+    assert np.linalg.norm(Y - Y_reconst) / np.linalg.norm(Y) < 1e-1