ged.rst

orphan:

Generalized eigendecomposition in decoding

This section describes the mathematical formulation and application of Generalized Eigendecomposition (GED), often used in spatial filtering and source separation algorithms, such as :class:`mne.decoding.CSP`, :class:`mne.decoding.SPoC`, :class:`mne.decoding.SSD` and :class:`mne.decoding.XdawnTransformer`.

The core principle of GED is to find a set of channel weights (spatial filter) that maximizes the ratio of signal power between two data features. These features are defined by the researcher and are represented by two covariance matrices: a "signal" matrix S and a "reference" matrix R. For example, S could be the covariance of data from a task time interval, and S could be the covariance from a baseline time interval. For more details see :footcite:`Cohen2022`.

Algebraic formulation of GED

A few definitions first: Let n \in \mathbb{N}^+ be a number of channels. Let \text{Symm}_n(\mathbb{R}) \subset M_n(\mathbb{R}) be a vector space of real symmetric matrices. Let S^n_+, S^n_{++} \subset \text{Symm}_n(\mathbb{R}) be sets of real positive semidefinite and positive definite matrices, respectively. Let S, R \in S^n_+ be covariance matrices estimated from electrophysiological data X_S \in M_{n \times t_S}(\mathbb{R}) and X_R \in M_{n \times t_R}(\mathbb{R}).

GED (or simultaneous diagonalization by congruence) of S and R is possible when R is full rank (and thus R \in S^n_{++}):

SW = RWD,

where W \in M_n(\mathbb{R}) is an invertible matrix of eigenvectors of (S, R) and D is a diagonal matrix of eigenvalues \lambda_i.

Each eigenvector \mathbf{w} \in W is a spatial filter that solves an optimization problem of the form:

\operatorname{argmax}_{\mathbf{w}} \frac{\mathbf{w}^t S \mathbf{w}}{\mathbf{w}^t R \mathbf{w}}

That is, using spatial filters W on time-series X \in M_{n \times t}(\mathbb{R}):

\mathbf{A} = W^t X,

results in "activation" time-series A of the estimated "sources", such that the ratio of their variances, \frac{\text{Var}(\mathbf{w}^T X_S)}{\text{Var}(\mathbf{w}^T X_R)} = \frac{\mathbf{w}^T S \mathbf{w}}{\mathbf{w}^T R \mathbf{w}}, is sequentially maximized spatial filters \mathbf{w}_i, sorted according to \lambda_i.

GED in the principal subspace

Unfortunately, R might not be full rank depending on the data X_R (for example due to average reference, removed PCA/ICA components, etc.). In such cases, GED can be performed on S and R in the principal subspace Q = \operatorname{Im}(C_{ref}) \subset \mathbb{R}^n of some reference covariance C_{ref} (in Common Spatial Pattern (CSP) algorithm, for example, C_{ref}=\frac{1}{2}(S+R) and GED is performed on S and R'=S+R).

More formally: Let r \leq n be a rank of C \in S^n_+. Let Q=\operatorname{Im}(C_{ref}) be a principal subspace of C_{ref}. Let P \in M_{n \times r}(\mathbb{R}) be an isometry formed by orthonormal basis of Q. Let f:S^n_+ \to S^r_+, A \mapsto P^t A P be a "restricting" map, that restricts quadratic form q_A:\mathbb{R}^n \to \mathbb{R} to q_{A|_Q}:\mathbb{R}^n \to \mathbb{R} (in practical terms, q_A maps spatial filters to variance of the spatially filtered data X_A).

Then, the GED of S and R in the principal subspace Q of C_{ref} is performed as follows:

1. S and R are transformed to S_Q = f(S) = P^t S P and R_Q = f(R) = P^t R P, such that S_Q and R_Q are matrix representations of restricted q_{S|_Q} and q_{R|_Q}.
2. GED is performed on S_Q and R_Q: S_Q W_Q = R_Q W_Q D.
3. Eigenvectors W_Q of (S_Q, R_Q) are transformed back to \mathbb{R}^n by W = P W_Q \in \mathbb{R}^{n \times r} to obtain r spatial filters.

Note that the solution to the original optimization problem is preserved:

\frac{\mathbf{w_Q}^t S_Q \mathbf{w_Q}}{\mathbf{w_Q}^t R_Q \mathbf{w_Q}}= \frac{\mathbf{w_Q}^t (P^t S P) \mathbf{w_Q}}{\mathbf{w_Q}^t (P^t R P)
\mathbf{w_Q}} = \frac{\mathbf{w}^t S \mathbf{w}}{\mathbf{w}^t R \mathbf{w}} = \lambda

In addition to restriction, q_S and q_R can be rescaled based on the whitened C_{ref}. In this case the whitening map f_{wh}:S^n_+ \to S^r_+, A \mapsto P_{wh}^t A P_{wh} transforms A into matrix representation of q_{A|Q} rescaled according to \Lambda^{-1/2}, where \Lambda is a diagonal matrix of eigenvalues of C_{ref} and so P_{wh} = P \Lambda^{-1/2}.

In MNE-Python, the matrix P of the restricting map can be obtained using

_, ref_evecs, mask = mne.cov._smart_eigh(C_ref, ..., proj_subspace=True, ...)
restr_mat = ref_evecs[mask]

while P_{wh} using:

restr_mat = compute_whitener(C_ref, ..., pca=True, ...)