[ENH] Adjusted SSD algorithm to support non-full rank data

mne/decoding/ssd.py

@@ -1,15 +1,21 @@
 # Author: Denis A. Engemann <denis.engemann@gmail.com>
 #         Victoria Peterson <victoriapeterson09@gmail.com>
+#         Thomas S. Binns <t.s.binns@outlook.com>
 # License: BSD-3-Clause
 
 import numpy as np
+import scipy as sp
 
-from ..filter import filter_data
-from ..cov import _regularized_covariance
 from . import TransformerMixin, BaseEstimator
-from ..time_frequency import psd_array_welch
-from ..utils import _time_mask, fill_doc, _validate_type, _check_option
+from ..cov import _regularized_covariance, Covariance
+from ..defaults import _handle_default
+from ..filter import filter_data
 from ..io.pick import _get_channel_types, _picks_to_idx
+from ..rank import compute_rank
+from ..time_frequency import psd_array_welch
+from ..utils import (
+    fill_doc, logger, _check_option, _time_mask, _validate_type,
+    _verbose_safe_false)
 
 
 @fill_doc
@@ -157,7 +163,6 @@ def fit(self, X, y=None):
         self : instance of SSD
             Returns the modified instance.
         """
-        from scipy.linalg import eigh
         self._check_X(X)
         X_aux = X[..., self.picks_, :]
 
@@ -172,17 +177,26 @@ def fit(self, X, y=None):
 
         cov_signal = _regularized_covariance(
             X_signal, reg=self.reg, method_params=self.cov_method_params,
-            rank=self.rank, info=self.info)
+            rank="full", info=self.info)
         cov_noise = _regularized_covariance(
             X_noise, reg=self.reg, method_params=self.cov_method_params,
-            rank=self.rank, info=self.info)
+            rank="full", info=self.info)
+
+        cov_signal_red, cov_noise_red, dim_red = (_dimensionality_reduction(
+            cov_signal, cov_noise, self.info, self.rank))
 
-        eigvals_, eigvects_ = eigh(cov_signal, cov_noise)
+        eigvals_, eigvects_ = sp.linalg.eigh(
+            cov_signal_red, cov_noise_red)
         # sort in descending order
         ix = np.argsort(eigvals_)[::-1]
-        self.eigvals_ = eigvals_[ix]
-        self.filters_ = eigvects_[:, ix]
-        self.patterns_ = np.linalg.pinv(self.filters_)
+        self.eigvals_ = np.real(eigvals_[ix])
+        filters_ = eigvects_[:, ix]
+        # restores dimensionality
+        self.filters_ = np.matmul(dim_red, filters_)
+        self.patterns_ = np.linalg.solve(
+            np.matmul(self.filters_.T, np.matmul(cov_signal, self.filters_)).T,
+            np.matmul(cov_signal, self.filters_).T)
+
         # We assume that ordering by spectral ratio is more important
         # than the initial ordering. This ordering should be also learned when
         # fitting.
@@ -191,6 +205,7 @@ def fit(self, X, y=None):
         if self.sort_by_spectral_ratio:
             _, sorter_spec = self.get_spectral_ratio(ssd_sources=X_ssd)
         self.sorter_spec = sorter_spec
+        logger.info("Done.")
         return self
 
     def transform(self, X):
@@ -290,3 +305,36 @@ def apply(self, X):
         pick_patterns = self.patterns_[self.sorter_spec][:self.n_components].T
         X = pick_patterns @ X_ssd
         return X
+
+
+def _dimensionality_reduction(cov_signal, cov_noise, info, rank):
+    """Perform dimensionality reduction on the covariance matrices."""
+    n_channels = cov_signal.shape[0]
+    rank_signal = list(compute_rank(
+        Covariance(cov_signal, info.ch_names, list(), list(), 0,
+                   verbose=_verbose_safe_false()),
+        rank, _handle_default('scalings_cov_rank', None), info).values())[0]
+    rank_noise = list(compute_rank(
+        Covariance(cov_signal, info.ch_names, list(), list(), 0,
+                   verbose=_verbose_safe_false()),
+        rank, _handle_default('scalings_cov_rank', None), info).values())[0]
+    rank = np.min([rank_signal, rank_noise])  # should be identical
+
+    if rank < n_channels:
+        eigvals, eigvects = sp.linalg.eigh(cov_signal)
+        # sort in descending order
+        ix = np.argsort(eigvals)[::-1]
+        eigvals = eigvals[ix]
+        eigvects = eigvects[:, ix]
+        # compute dimensionality reduction transformation matrix
+        dim_red = np.matmul(
+            eigvects[:, :rank], np.eye(rank) * (eigvals[:rank]**-0.5))
+        logger.info(
+            "Reducing covariance rank from %i -> %i" % (n_channels, rank,))
+    else:
+        dim_red = np.eye(n_channels)
+        logger.info("Preserving covariance rank (%i)" % (rank,))
+
+    cov_signal_red = np.matmul(dim_red.T, np.matmul(cov_signal, dim_red))
+    cov_noise_red = np.matmul(dim_red.T, np.matmul(cov_noise, dim_red))
+    return cov_signal_red, cov_noise_red, dim_red

mne/decoding/tests/test_ssd.py

@@ -1,5 +1,6 @@
 # Author: Denis A. Engemann <denis.engemann@gmail.com>
 #         Victoria Peterson <victoriapeterson09@gmail.com>
+#         Thomas S. Binns <t.s.binns@outlook.com>
 # License: BSD-3-Clause
 
 import numpy as np
@@ -322,3 +323,23 @@ def test_return_filtered():
     psd_out, freqs = psd_array_welch(out[0], sfreq=250, n_fft=250)
     freqs_up = int(freqs[psd_out > 0.5][0]), int(freqs[psd_out > 0.5][-1])
     assert (freqs_up != freqs_sig)
+
+
+def test_non_full_rank_data():
+    """Test that the method works with non-full rank data."""
+    n_channels = 10
+    X, _, _ = simulate_data(SNR=0.9, freqs_sig=[4, 13], n_channels=n_channels)
+    info = create_info(ch_names=n_channels, sfreq=250, ch_types='eeg')
+
+    filt_params_signal = dict(l_freq=freqs_sig[0], h_freq=freqs_sig[1],
+                              l_trans_bandwidth=1, h_trans_bandwidth=1)
+    filt_params_noise = dict(l_freq=freqs_noise[0], h_freq=freqs_noise[1],
+                             l_trans_bandwidth=1, h_trans_bandwidth=1)
+
+    # Make data non-full rank
+    rank = 5
+    X[rank:] = X[:rank]  # an extreme example, but a valid one
+    assert np.linalg.matrix_rank(X) == rank
+
+    ssd = SSD(info, filt_params_signal, filt_params_noise)
+    ssd.fit(X)