[MRG before #12069] KernelPCA: raise Errors and Warnings according to eigenvalue decomposition numerical/conditioning issues

doc/whats_new/v0.22.rst

@@ -178,6 +178,13 @@ Changelog
 :mod:`sklearn.cross_decomposition`
 ..................................
 
+- |Enhancement| :class:`decomposition.KernelPCA` now properly checks the
+  eigenvalues found by the solver for numerical or conditioning issues. This
+  ensures consistency of results across solvers (different choices for
+  ``eigen_solver``), including approximate solvers such as ``'randomized'`` and
+  ``'lobpcg'`` (see :issue:`12068`).
+  :pr:`12145` by :user:`Sylvain Marié <smarie>`
+
 - |Fix| Fixed a bug where :class:`cross_decomposition.PLSCanonical` and
   :class:`cross_decomposition.PLSRegression` were raising an error when fitted
   with a target matrix `Y` in which the first column was constant.

sklearn/decomposition/_kernel_pca.py

@@ -9,7 +9,8 @@
 
 from ..utils import check_random_state
 from ..utils.extmath import svd_flip
-from ..utils.validation import check_is_fitted, check_array
+from ..utils.validation import (check_is_fitted, check_array,
+                                _check_psd_eigenvalues)
 from ..exceptions import NotFittedError
 from ..base import BaseEstimator, TransformerMixin
 from ..preprocessing import KernelCenterer
@@ -211,6 +212,10 @@ def _fit_transform(self, K):
                                                 maxiter=self.max_iter,
                                                 v0=v0)
 
+        # make sure that the eigenvalues are ok and fix numerical issues
+        self.lambdas_ = _check_psd_eigenvalues(self.lambdas_,
+                                               enable_warnings=False)
+
         # flip eigenvectors' sign to enforce deterministic output
         self.alphas_, _ = svd_flip(self.alphas_,
                                    np.empty_like(self.alphas_).T)

sklearn/decomposition/tests/test_kernel_pca.py

@@ -11,6 +11,7 @@
 from sklearn.pipeline import Pipeline
 from sklearn.model_selection import GridSearchCV
 from sklearn.metrics.pairwise import rbf_kernel
+from sklearn.utils.validation import _check_psd_eigenvalues
 
 
 def test_kernel_pca():
@@ -270,3 +271,19 @@ def test_nested_circles():
     # The data is perfectly linearly separable in that space
     train_score = Perceptron(max_iter=5).fit(X_kpca, y).score(X_kpca, y)
     assert train_score == 1.0
+
+
+def test_kernel_conditioning():
+    """Test that ``_check_psd_eigenvalues`` is correctly called."""
+
+    # create a pathological X leading to small non-zero eigenvalue
+    X = [[5, 1],
+         [5+1e-8, 1e-8],
+         [5+1e-8, 0]]
+    kpca = KernelPCA(kernel="linear", n_components=2,
+                     fit_inverse_transform=True)
+    kpca.fit(X)
+
+    # check that the small non-zero eigenvalue was correctly set to zero
+    assert kpca.lambdas_.min() == 0
+    assert np.all(kpca.lambdas_ == _check_psd_eigenvalues(kpca.lambdas_))

sklearn/exceptions.py

@@ -12,7 +12,8 @@
            'FitFailedWarning',
            'NonBLASDotWarning',
            'SkipTestWarning',
-           'UndefinedMetricWarning']
+           'UndefinedMetricWarning',
+           'PositiveSpectrumWarning']
 
 
 class NotFittedError(ValueError, AttributeError):
@@ -171,3 +172,15 @@ class UndefinedMetricWarning(UserWarning):
     .. versionchanged:: 0.18
        Moved from sklearn.base.
     """
+
+
+class PositiveSpectrumWarning(UserWarning):
+    """Warning raised when the eigenvalues of a PSD matrix have issues
+
+    This warning is typically raised by ``_check_psd_eigenvalues`` when the
+    eigenvalues of a positive semidefinite (PSD) matrix such as a gram matrix
+    (kernel) present significant negative eigenvalues, or bad conditioning i.e.
+    very small non-zero eigenvalues compared to the largest eigenvalue.
+
+    .. versionadded:: 0.22
+    """

sklearn/utils/tests/test_validation.py

@@ -41,13 +41,15 @@
     check_non_negative,
     _num_samples,
     check_scalar,
+    _check_psd_eigenvalues,
     _deprecate_positional_args,
     _check_sample_weight,
     _allclose_dense_sparse,
     FLOAT_DTYPES)
+
 import sklearn
 
-from sklearn.exceptions import NotFittedError
+from sklearn.exceptions import NotFittedError, PositiveSpectrumWarning
 from sklearn.exceptions import DataConversionWarning
 
 from sklearn.utils._testing import assert_raise_message
@@ -937,6 +939,81 @@ def test_check_scalar_invalid(x, target_name, target_type, min_val, max_val,
     assert type(raised_error.value) == type(err_msg)
 
 
+_psd_cases_valid = {
+    'nominal': ((1, 2), np.array([1, 2]), None, ""),
+    'nominal_np_array': (np.array([1, 2]), np.array([1, 2]), None, ""),
+    'insignificant_imag': ((5, 5e-5j), np.array([5, 0]),
+                           PositiveSpectrumWarning,
+                           "There are imaginary parts in eigenvalues "
+                           "\\(1e\\-05 of the maximum real part"),
+    'insignificant neg': ((5, -5e-5), np.array([5, 0]),
+                          PositiveSpectrumWarning, ""),
+    'insignificant neg float32': (np.array([1, -1e-6], dtype=np.float32),
+                                  np.array([1, 0], dtype=np.float32),
+                                  PositiveSpectrumWarning,
+                                  "There are negative eigenvalues \\(1e\\-06 "
+                                  "of the maximum positive"),
+    'insignificant neg float64': (np.array([1, -1e-10], dtype=np.float64),
+                                  np.array([1, 0], dtype=np.float64),
+                                  PositiveSpectrumWarning,
+                                  "There are negative eigenvalues \\(1e\\-10 "
+                                  "of the maximum positive"),
+    'insignificant pos': ((5, 4e-12), np.array([5, 0]),
+                          PositiveSpectrumWarning,
+                          "the largest eigenvalue is more than 1e\\+12 "
+                          "times the smallest"),
+}
+
+
+@pytest.mark.parametrize("lambdas, expected_lambdas, w_type, w_msg",
+                         list(_psd_cases_valid.values()),
+                         ids=list(_psd_cases_valid.keys()))
+@pytest.mark.parametrize("enable_warnings", [True, False])
+def test_check_psd_eigenvalues_valid(lambdas, expected_lambdas, w_type, w_msg,
+                                     enable_warnings):
+    # Test that ``_check_psd_eigenvalues`` returns the right output for valid
+    # input, possibly raising the right warning
+
+    if not enable_warnings:
+        w_type = None
+        w_msg = ""
+
+    with pytest.warns(w_type, match=w_msg) as w:
+        assert_array_equal(
+            _check_psd_eigenvalues(lambdas, enable_warnings=enable_warnings),
+            expected_lambdas
+        )
+    if w_type is None or not enable_warnings:
+        assert not w
+
+
+_psd_cases_invalid = {
+    'significant_imag': ((5, 5j), ValueError,
+                         "There are significant imaginary parts in eigenv"),
+    'all negative': ((-5, -1), ValueError,
+                     "All eigenvalues are negative \\(maximum is -1"),
+    'significant neg': ((5, -1), ValueError,
+                        "There are significant negative eigenvalues"),
+    'significant neg float32': (np.array([3e-4, -2e-6], dtype=np.float32),
+                                ValueError,
+                                "There are significant negative eigenvalues"),
+    'significant neg float64': (np.array([1e-5, -2e-10], dtype=np.float64),
+                                ValueError,
+                                "There are significant negative eigenvalues"),
+}
+
+
+@pytest.mark.parametrize("lambdas, err_type, err_msg",
+                         list(_psd_cases_invalid.values()),
+                         ids=list(_psd_cases_invalid.keys()))
+def test_check_psd_eigenvalues_invalid(lambdas, err_type, err_msg):
+    # Test that ``_check_psd_eigenvalues`` raises the right error for invalid
+    # input
+
+    with pytest.raises(err_type, match=err_msg):
+        _check_psd_eigenvalues(lambdas)
+
+
 def test_check_sample_weight():
     # check array order
     sample_weight = np.ones(10)[::2]

sklearn/utils/validation.py

@@ -6,6 +6,7 @@
 #          Lars Buitinck
 #          Alexandre Gramfort
 #          Nicolas Tresegnie
+#          Sylvain Marie
 # License: BSD 3 clause
 
 from functools import wraps
@@ -22,7 +23,7 @@
 
 from .fixes import _object_dtype_isnan
 from .. import get_config as _get_config
-from ..exceptions import NonBLASDotWarning
+from ..exceptions import NonBLASDotWarning, PositiveSpectrumWarning
 from ..exceptions import NotFittedError
 from ..exceptions import DataConversionWarning
 
@@ -1020,6 +1021,169 @@ def check_scalar(x, name, target_type, min_val=None, max_val=None):
         raise ValueError('`{}`= {}, must be <= {}.'.format(name, x, max_val))
 
 
+def _check_psd_eigenvalues(lambdas, enable_warnings=False):
+    """Check the eigenvalues of a positive semidefinite (PSD) matrix.
+
+    Checks the provided array of PSD matrix eigenvalues for numerical or
+    conditioning issues and returns a fixed validated version. This method
+    should typically be used if the PSD matrix is user-provided (e.g. a
+    Gram matrix) or computed using a user-provided dissimilarity metric
+    (e.g. kernel function), or if the decomposition process uses approximation
+    methods (randomized SVD, etc.).
+
+    It checks for three things:
+
+    - that there are no significant imaginary parts in eigenvalues (more than
+      1e-5 times the maximum real part). If this check fails, it raises a
+      ``ValueError``. Otherwise all non-significant imaginary parts that may
+      remain are set to zero. This operation is traced with a
+      ``PositiveSpectrumWarning`` when ``enable_warnings=True``.
+
+    - that eigenvalues are not all negative. If this check fails, it raises a
+      ``ValueError``
+
+    - that there are no significant negative eigenvalues with absolute value
+      more than 1e-10 (1e-6) and more than 1e-5 (5e-3) times the largest
+      positive eigenvalue in double (simple) precision. If this check fails,
+      it raises a ``ValueError``. Otherwise all negative eigenvalues that may
+      remain are set to zero. This operation is traced with a
+      ``PositiveSpectrumWarning`` when ``enable_warnings=True``.
+
+    Finally, all the positive eigenvalues that are too small (with a value
+    smaller than the maximum eigenvalue divided by 1e12) are set to zero.
+    This operation is traced with a ``PositiveSpectrumWarning`` when
+    ``enable_warnings=True``.
+
+    Parameters
+    ----------
+    lambdas : array-like of shape (n_eigenvalues,)
+        Array of eigenvalues to check / fix.
+
+    enable_warnings : bool, default=False
+        When this is set to ``True``, a ``PositiveSpectrumWarning`` will be
+        raised when there are imaginary parts, negative eigenvalues, or
+        extremely small non-zero eigenvalues. Otherwise no warning will be
+        raised. In both cases, imaginary parts, negative eigenvalues, and
+        extremely small non-zero eigenvalues will be set to zero.
+
+    Returns
+    -------
+    lambdas_fixed : ndarray of shape (n_eigenvalues,)
+        A fixed validated copy of the array of eigenvalues.
+
+    Examples
+    --------
+    >>> _check_psd_eigenvalues([1, 2])      # nominal case
+    array([1, 2])
+    >>> _check_psd_eigenvalues([5, 5j])     # significant imag part
+    Traceback (most recent call last):
+        ...
+    ValueError: There are significant imaginary parts in eigenvalues (1
+        of the maximum real part). Either the matrix is not PSD, or there was
+        an issue while computing the eigendecomposition of the matrix.
+    >>> _check_psd_eigenvalues([5, 5e-5j])  # insignificant imag part
+    array([5., 0.])
+    >>> _check_psd_eigenvalues([-5, -1])    # all negative
+    Traceback (most recent call last):
+        ...
+    ValueError: All eigenvalues are negative (maximum is -1). Either the
+        matrix is not PSD, or there was an issue while computing the
+        eigendecomposition of the matrix.
+    >>> _check_psd_eigenvalues([5, -1])     # significant negative
+    Traceback (most recent call last):
+        ...
+    ValueError: There are significant negative eigenvalues (0.2 of the
+        maximum positive). Either the matrix is not PSD, or there was an issue
+        while computing the eigendecomposition of the matrix.
+    >>> _check_psd_eigenvalues([5, -5e-5])  # insignificant negative
+    array([5., 0.])
+    >>> _check_psd_eigenvalues([5, 4e-12])  # bad conditioning (too small)
+    array([5., 0.])
+
+    """
+
+    lambdas = np.array(lambdas)
+    is_double_precision = lambdas.dtype == np.float64
+
+    # note: the minimum value available is
+    #  - single-precision: np.finfo('float32').eps = 1.2e-07
+    #  - double-precision: np.finfo('float64').eps = 2.2e-16
+
+    # the various thresholds used for validation
+    # we may wish to change the value according to precision.
+    significant_imag_ratio = 1e-5
+    significant_neg_ratio = 1e-5 if is_double_precision else 5e-3
+    significant_neg_value = 1e-10 if is_double_precision else 1e-6
+    small_pos_ratio = 1e-12
+
+    # Check that there are no significant imaginary parts
+    if not np.isreal(lambdas).all():
+        max_imag_abs = np.abs(np.imag(lambdas)).max()
+        max_real_abs = np.abs(np.real(lambdas)).max()
+        if max_imag_abs > significant_imag_ratio * max_real_abs:
+            raise ValueError(
+                "There are significant imaginary parts in eigenvalues (%g "
+                "of the maximum real part). Either the matrix is not PSD, or "
+                "there was an issue while computing the eigendecomposition "
+                "of the matrix."
+                % (max_imag_abs / max_real_abs))
+
+        # warn about imaginary parts being removed
+        if enable_warnings:
+            warnings.warn("There are imaginary parts in eigenvalues (%g "
+                          "of the maximum real part). Either the matrix is not"
+                          " PSD, or there was an issue while computing the "
+                          "eigendecomposition of the matrix. Only the real "
+                          "parts will be kept."
+                          % (max_imag_abs / max_real_abs),
+                          PositiveSpectrumWarning)
+
+    # Remove all imaginary parts (even if zero)
+    lambdas = np.real(lambdas)
+
+    # Check that there are no significant negative eigenvalues
+    max_eig = lambdas.max()
+    if max_eig < 0:
+        raise ValueError("All eigenvalues are negative (maximum is %g). "
+                         "Either the matrix is not PSD, or there was an "
+                         "issue while computing the eigendecomposition of "
+                         "the matrix." % max_eig)
+
+    else:
+        min_eig = lambdas.min()
+        if (min_eig < -significant_neg_ratio * max_eig
+                and min_eig < -significant_neg_value):
+            raise ValueError("There are significant negative eigenvalues (%g"
+                             " of the maximum positive). Either the matrix is "
+                             "not PSD, or there was an issue while computing "
+                             "the eigendecomposition of the matrix."
+                             % (-min_eig / max_eig))
+        elif min_eig < 0:
+            # Remove all negative values and warn about it
+            if enable_warnings:
+                warnings.warn("There are negative eigenvalues (%g of the "
+                              "maximum positive). Either the matrix is not "
+                              "PSD, or there was an issue while computing the"
+                              " eigendecomposition of the matrix. Negative "
+                              "eigenvalues will be replaced with 0."
+                              % (-min_eig / max_eig),
+                              PositiveSpectrumWarning)
+            lambdas[lambdas < 0] = 0
+
+    # Check for conditioning (small positive non-zeros)
+    too_small_lambdas = (0 < lambdas) & (lambdas < small_pos_ratio * max_eig)
+    if too_small_lambdas.any():
+        if enable_warnings:
+            warnings.warn("Badly conditioned PSD matrix spectrum: the largest "
+                          "eigenvalue is more than %g times the smallest. "
+                          "Small eigenvalues will be replaced with 0."
+                          "" % (1 / small_pos_ratio),
+                          PositiveSpectrumWarning)
+        lambdas[too_small_lambdas] = 0
+
+    return lambdas
+
+
 def _check_sample_weight(sample_weight, X, dtype=None):
     """Validate sample weights.