BUG: Fix GED tolerances and improve cov validation (#13346)

Co-authored-by: Eric Larson <larson.eric.d@gmail.com>

Author: Genuster

mne/decoding/_ged.py

@@ -54,28 +54,39 @@ def _smart_ged(S, R, restr_mat=None, R_func=None):
     return evals, evecs
 
 
-def _is_cov_symm_pos_semidef(
-    cov, rtol=1e-10, atol=1e-11, eval_tol=1e-15, check_pos_semidef=True
-):
+def _is_cov_symm(cov, rtol=1e-7, atol=None):
+    if atol is None:
+        atol = 1e-7 * np.max(np.abs(cov))
     is_symm = scipy.linalg.issymmetric(cov, rtol=rtol, atol=atol)
-    if not is_symm:
-        return False
+    return is_symm
 
-    if check_pos_semidef:
-        # numerically slightly negative evals are considered 0
-        is_pos_semidef = np.all(scipy.linalg.eigvalsh(cov) >= -eval_tol)
-        return is_pos_semidef
 
-    return True
+def _get_cov_def(cov, eval_tol=None):
+    """Get definiteness of symmetric cov matrix.
 
+    All evals in (-eval_tol, eval_tol) will be considered zero,
+    while all evals smaller than -eval_tol will be considered
+    negative.
+    """
+    evals = scipy.linalg.eigvalsh(cov)
+    if eval_tol is None:
+        eval_tol = 1e-7 * np.max(np.abs(evals))
+    if np.all(evals > eval_tol):
+        return "pos_def"
+    elif np.all(evals >= -eval_tol):
+        return "pos_semidef"
+    else:
+        return "indef"
+
+
+def _is_cov_pos_semidef(cov, eval_tol=None):
+    cov_def = _get_cov_def(cov, eval_tol=eval_tol)
+    return cov_def in ("pos_def", "pos_semidef")
 
-def _is_cov_pos_def(cov, eval_tol=1e-15):
-    is_symm = _is_cov_symm_pos_semidef(cov, check_pos_semidef=False)
-    if not is_symm:
-        return False
-    # numerically slightly positive evals are considered 0
-    is_pos_def = np.all(scipy.linalg.eigvalsh(cov) > eval_tol)
-    return is_pos_def
+
+def _is_cov_pos_def(cov, eval_tol=None):
+    cov_def = _get_cov_def(cov, eval_tol=eval_tol)
+    return cov_def == "pos_def"
 
 
 def _smart_ajd(covs, restr_mat=None, weights=None):
@@ -90,8 +101,8 @@ def _smart_ajd(covs, restr_mat=None, weights=None):
     from .csp import _ajd_pham
 
     if restr_mat is None:
-        is_all_pos_def = all([_is_cov_pos_def(cov) for cov in covs])
-        if not is_all_pos_def:
+        are_all_pos_def = all([_is_cov_pos_def(cov) for cov in covs])
+        if not are_all_pos_def:
             raise ValueError(
                 "If C_ref is not provided by covariance estimator, "
                 "all the covs should be positive definite"
@@ -100,6 +111,12 @@ def _smart_ajd(covs, restr_mat=None, weights=None):
         return evecs
 
     else:
+        are_all_pos_semidef = all([_is_cov_pos_semidef(cov) for cov in covs])
+        if not are_all_pos_semidef:
+            raise ValueError(
+                "All the covs should be positive semi-definite for "
+                "approximate joint diagonalization"
+            )
         covs = np.array([restr_mat @ cov @ restr_mat.T for cov in covs], float)
         evecs_restr, D = _ajd_pham(covs)
         evecs = _normalize_eigenvectors(evecs_restr.T, covs, weights)

mne/decoding/base.py

@@ -26,7 +26,13 @@
 
 from ..parallel import parallel_func
 from ..utils import _check_option, _pl, _validate_type, logger, pinv, verbose, warn
-from ._ged import _handle_restr_mat, _is_cov_symm_pos_semidef, _smart_ajd, _smart_ged
+from ._ged import (
+    _handle_restr_mat,
+    _is_cov_pos_semidef,
+    _is_cov_symm,
+    _smart_ajd,
+    _smart_ged,
+)
 from ._mod_ged import _no_op_mod
 from .transformer import MNETransformerMixin
 
@@ -133,23 +139,21 @@ def fit(self, X, y=None):
         covs, C_ref, info, rank, kwargs = self.cov_callable(X, y)
         covs = np.stack(covs)
         self._validate_covariances(covs)
-        self._validate_covariances([C_ref])
+        if C_ref is not None:
+            self._validate_covariances([C_ref])
         mod_ged_callable = (
             self.mod_ged_callable if self.mod_ged_callable is not None else _no_op_mod
         )
+        restr_mat = _handle_restr_mat(C_ref, self.restr_type, info, rank)
 
         if self.dec_type == "single":
             if len(covs) > 2:
-                weights = (
-                    kwargs["sample_weights"] if "sample_weights" in kwargs else None
-                )
-                restr_mat = _handle_restr_mat(C_ref, self.restr_type, info, rank)
+                weights = kwargs.get("sample_weights", None)
                 evecs = _smart_ajd(covs, restr_mat, weights=weights)
                 evals = None
             else:
                 S = covs[0]
                 R = covs[1]
-                restr_mat = _handle_restr_mat(C_ref, self.restr_type, info, rank)
                 evals, evecs = _smart_ged(S, R, restr_mat, R_func=self.R_func)
 
             evals, evecs, self.sorter_ = mod_ged_callable(evals, evecs, covs, **kwargs)
@@ -160,7 +164,6 @@ def fit(self, X, y=None):
         elif self.dec_type == "multi":
             self.classes_ = np.unique(y)
             R = covs[-1]
-            restr_mat = _handle_restr_mat(C_ref, self.restr_type, info, rank)
             all_evals, all_evecs = list(), list()
             all_patterns, all_sorters = list(), list()
             for i in range(len(self.classes_)):
@@ -251,18 +254,66 @@ def _validate_ged_params(self):
         )
 
     def _validate_covariances(self, covs):
-        for cov in covs:
-            if cov is None:
-                continue
-            # XXX: A lot of mne.decoding classes use mne.cov._regularized_covariance.
-            # Depending on the data it sometimes returns negative semidefinite matrices.
-            # So adding the validation of positive semidefinitiveness
-            # will require overhauling covariance estimation first.
-            is_cov = _is_cov_symm_pos_semidef(cov, check_pos_semidef=False)
-            if not is_cov:
+        error_template = (
+            "{matrix} is not {prop}, but required to be for {decomp}. "
+            "Check your cov_callable"
+        )
+        if len(covs) == 1:
+            C_ref = covs[0]
+            is_C_ref_symm = _is_cov_symm(C_ref)
+            if not is_C_ref_symm:
+                raise ValueError(
+                    error_template.format(
+                        matrix="C_ref covariance",
+                        prop="symmetric",
+                        decomp="decomposition",
+                    )
+                )
+        elif self.dec_type == "single" and len(covs) > 2:
+            # make only lenient symmetric check here.
+            # positive semidefiniteness/definiteness will be
+            # checked inside _smart_ajd
+            for ci, cov in enumerate(covs):
+                if not _is_cov_symm(cov):
+                    raise ValueError(
+                        error_template.format(
+                            matrix=f"cov[{ci}]",
+                            prop="symmetric",
+                            decomp="approximate joint diagonalization",
+                        )
+                    )
+        else:
+            if len(covs) == 2:
+                S_covs = [covs[0]]
+                R = covs[1]
+            elif self.dec_type == "multi":
+                S_covs = covs[:-1]
+                R = covs[-1]
+
+            are_all_S_symm = all([_is_cov_symm(S) for S in S_covs])
+            if not are_all_S_symm:
+                raise ValueError(
+                    error_template.format(
+                        matrix="S covariance",
+                        prop="symmetric",
+                        decomp="generalized eigendecomposition",
+                    )
+                )
+            if not _is_cov_symm(R):
+                raise ValueError(
+                    error_template.format(
+                        matrix="R covariance",
+                        prop="symmetric",
+                        decomp="generalized eigendecomposition",
+                    )
+                )
+            if not _is_cov_pos_semidef(R):
                 raise ValueError(
-                    "One of covariances is not symmetric (or positive semidefinite), "
-                    "check your cov_callable"
+                    error_template.format(
+                        matrix="R covariance",
+                        prop="positive semi-definite",
+                        decomp="generalized eigendecomposition",
+                    )
                 )
 
     def __sklearn_tags__(self):

mne/decoding/tests/test_ged.py

@@ -19,10 +19,10 @@
 from mne._fiff.proj import make_eeg_average_ref_proj
 from mne.cov import Covariance, _regularized_covariance
 from mne.decoding._ged import (
+    _get_cov_def,
     _get_restr_mat,
     _handle_restr_mat,
-    _is_cov_pos_def,
-    _is_cov_symm_pos_semidef,
+    _is_cov_symm,
     _smart_ajd,
     _smart_ged,
 )
@@ -345,34 +345,27 @@ def test__handle_restr_mat_invalid_restr_type():
 
 def test_cov_validators():
     """Test that covariance validators indeed validate."""
-    asymm = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
+    asymm_indef = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
     sing_pos_semidef = np.array([[1, 2, 3], [2, 4, 6], [3, 6, 9]])
     pos_def = np.array([[5, 1, 1], [1, 6, 2], [1, 2, 7]])
 
-    assert not _is_cov_symm_pos_semidef(asymm)
-    assert _is_cov_symm_pos_semidef(sing_pos_semidef)
-    assert _is_cov_symm_pos_semidef(pos_def)
+    assert not _is_cov_symm(asymm_indef)
+    assert _get_cov_def(asymm_indef) == "indef"
+    assert _get_cov_def(sing_pos_semidef) == "pos_semidef"
+    assert _get_cov_def(pos_def) == "pos_def"
 
-    assert not _is_cov_pos_def(asymm)
-    assert not _is_cov_pos_def(sing_pos_semidef)
-    assert _is_cov_pos_def(pos_def)
 
-
-def test__is_cov_pos_def():
-    """Test _is_cov_pos_def works."""
-    asymm = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
-    sing_pos_semidef = np.array([[1, 2, 3], [2, 4, 6], [3, 6, 9]])
-    pos_def = np.array([[5, 1, 1], [1, 6, 2], [1, 2, 7]])
-    assert not _is_cov_pos_def(asymm)
-    assert not _is_cov_pos_def(sing_pos_semidef)
-    assert _is_cov_pos_def(pos_def)
-
-
-def test__smart_ajd_when_restr_mat_is_none():
-    """Test _smart_ajd raises ValueError when restr_mat is None."""
+def test__smart_ajd_raises():
+    """Test _smart_ajd raises proper ValueErrors."""
+    asymm_indef = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
     sing_pos_semidef = np.array([[1, 2, 3], [2, 4, 6], [3, 6, 9]])
     pos_def1 = np.array([[5, 1, 1], [1, 6, 2], [1, 2, 7]])
     pos_def2 = np.array([[10, 1, 2], [1, 12, 3], [2, 3, 15]])
+
+    bad_covs = np.stack([sing_pos_semidef, asymm_indef, pos_def1])
+    with pytest.raises(ValueError, match="positive semi-definite"):
+        _smart_ajd(bad_covs, restr_mat=pos_def2, weights=None)
+
     bad_covs = np.stack([sing_pos_semidef, pos_def1, pos_def2])
     with pytest.raises(ValueError, match="positive definite"):
         _smart_ajd(bad_covs, restr_mat=None, weights=None)