DEV Community: Ethan Davis

Linear Discriminant Analysis (LDA): Decision Boundaries

Ethan Davis — Sat, 11 Jul 2026 04:16:10 +0000

Linear Discriminant Analysis

Decision Boundary Derivation

Decision theory for classification tells us that we need to know the class posteriors $P (G ∣ X)$ for optimal classification. Suppose $f_{k} (x)$ is the conditional density of $X$ for class $G = k$ , and let $π_{k}$ be the prior probability of class $k$ , with $\sum_{ℓ = 1 K} π_{ℓ} = 1$ . Bayes theorem gives us the following [1].

P (G = k ∣ X = x) = \frac{f _{k} ( x ) π _{k}}{\sum _{ℓ = 1 K} f _{ℓ} ( x ) π _{ℓ}} .

Many machine learning techniques are based on models for the class densities. For example, linear and quadratic discriminant analysis use Gaussian densities. More flexible Gaussian mixture models (GMMs) allow for nonlinear decision boundaries. General nonparametric density estimates for each class density allow the most flexibility. Naive Bayes models are a variant of the previous case, and assume that the inputs are conditionally independent in each class [1].

Suppose that each class density is modeled as a multivariate Gaussian [1].

f_{k} (x) = (2 π)^{- p /2} ∣ Σ_{k} ∣^{- 1/2} exp (- \frac{1}{2} (x - μ_{k})^{⊤} Σ_{k - 1} (x - μ_{k})) .

Linear discriminant analysis (LDA) is the special case when we assume that the classes have a common covariance matrix $Σ_{k} = Σ \forall k$ . When we compare the log-ratio of two classes $k$ and $ℓ$ we get an equation linear in $x$ . The equal covariance matrices cause the normalization factors to cancel, as well as the quadratic part in the exponents. This linear log-odds function implies that the decision boundary between classes $k$ and $ℓ$ is linear in $x$ [1].

lo g \frac{P ( G = k ∣ X = x )}{P ( G = ℓ ∣ X = x )} = = lo g \frac{f _{k} ( x )}{f _{ℓ} ( x )} + lo g \frac{π _{k}}{π _{ℓ}} lo g \frac{π _{k}}{π _{ℓ}} - \frac{1}{2} (μ_{k} + μ_{ℓ})^{⊤} Σ^{- 1} (μ_{k} - μ_{ℓ}) + x^{⊤} Σ^{- 1} (μ_{k} - μ_{ℓ}) .

The linear discriminant functions are an equivalent description of the decision rule, with $G (x) = argmax_{k} δ_{k} (x)$ [1].

δ_{k} (x) = x^{⊤} Σ^{- 1} μ_{k} - \frac{1}{2} μ_{k ⊤} Σ^{- 1} μ_{k} + lo g π_{k} .

In practice we do not know the parameters of the Gaussian distributions, and need to estimate them using training data [1].

$π^_{k} = N_{k} / N$ , where $N_{k}$ is the number of observations from class $k$ ,
$μ^_{k} = \sum_{g_{i} = k} x_{i} / N_{k}$ ,
$Σ^= \sum_{k = 1 K} \sum_{g_{i} = k} (x_{i} - μ^_{k}) (x_{i} - μ^_{k})^{⊤} / (N - K)$ .

In the general discriminant problem, if the $Σ_{k}$ are not assumed to be equal, then the convenient cancellations in the equation do not take place. In particular the pieces quadratic in $x$ remain. We then get quadratic discriminant analysis (QDA). The decision boundary between each pair of classes $k$ and $ℓ$ is described by a quadratic equation $x : δ_{k} (x) = δ_{ℓ} (x)$ [1].

δ_{k} (x) = - \frac{1}{2} lo g ∣ Σ_{k} ∣ - \frac{1}{2} (x - μ_{k})^{⊤} Σ_{k - 1} (x - μ_{k}) + lo g π_{k} .

Both LDA and QDA perform well on an amazingly large and diverse set of classification tasks. The question is why LDA and QDA have such a good track record. The reason is likely not that the data are approximately Gaussian, and in addition for LDA that the covariances are approximately equal. More likely the reason is that the data can only support simple decision boundaries such as linear or quadratic, and the estimates provided by the Gaussian models are stable. This is a bias variance tradeoff where we put up with the bias of a linear decision boundary because it can be estimated with much lower variance than more exotic alternatives [1].

References

Trevor Hastie, Robert Tibshirani, Jerome Friedman (2009) The Elements of Statistical Learning. Springer New York.

Non-centered Parameterization for MCMC

Ethan Davis — Sat, 11 Jul 2026 04:00:31 +0000

Non-centered Parameterization

For relatively simple models, samplers mostly just work, providing us with an accurate estimate of the posterior. However, certain posterior geometries are challenging for samplers. A common example is Neal’s funnel. The shape at the top of the funnel is quite wide before narrowing into a small neck. Samplers function by taking steps from one set of parameter values to another, and a key setting is how large of a step to take when exploring the posterior surface. In complex geometries, such as with Neal’s funnel, a step size that works well in one area fails in another [1].

In hierarchical models the posterior geometry is largely defined by the correlation of hyperpriors to other parameters, which can result in funnel geometry that are difficult to sample. These inconvenient posteriors can create inefficient MCMC caused by steep curvature in the negative log-likelihood and lead to divergent transitions. Transforming the prior to be smoother helps. There is a relatively simple tweak to models, referred to as a non-centered parameterization, that helps alleviate the issue [1].

Example: Devil’s Funnel prior

v x \sim Normal (0, σ_{v}) \sim Normal (0, exp (v)) .

As $σ_{v}$ increases, a trough forms in the prior probability surface known as centered parameterization that is difficult for Hamiltonian dynamics to sample [2, 3].

To handle this, we could configure a smaller steps size to manage the steepness better, but then sampling would take longer to explore the posterior. Alternatively, we can re-parameterize (transform) to make the surface smoother known non-centered parameterization. For this purpose, we add an auxiliary variable $z$ that has a smooth probability surface. We then sample that auxiliary variable, and transform it to obtain the target variable distribution. For the case of the devil’s funnel prior we do the following [2, 3].

v z x \sim Normal (0, σ_{v}) \sim Normal (0, 1) = z exp (v) .

We can also transform sampling using Cholesky factors. The Cholesky factor $L$ provides an efficient means for encoding a correlation matrix $Ω$ since it requires fewer floating point numbers than the full correlation matrix. It is defined by $Ω = L L^{⊤}$ where $L$ is a lower-triangular matrix. We can sample data with correlation $Ω$ , and standard deviation $σ_{i}$ , using the Cholesky factorization $L$ of $Ω$ as follows where $Z$ is a matrix of $z$ -scores sampled from a standard normal distribution [2, 3].

X_{Ω} \sim diag (σ_{i}) L Z .

References

Martin, Osvaldo A., Kumar, Ravin, Lao, Junpeng (2021) Bayesian Modeling and Computation in Python. Chapman and Hall/CRC.
McElreath, Richard (2020) Statistical Rethinking: A Bayesian Course with Examples in R and STAN. Chapman and Hall/CRC.
Dustin Stansbury, Richard McElreath (2024) Correlated Features. Zenodo.

Diagnosing MCMC: ESS, R-hat, MCSE, and Plots

Ethan Davis — Sat, 11 Jul 2026 03:57:57 +0000

Runtime Diagnostics

The following numerical and visual diagnostic tools are illustrated with three synthetic posteriors. The first named good_chains was generated from independent and identically distributed (IID) draws from a $Beta (2, 5)$ . Samples that are IID are ideal for approximating the posterior [1].

The second is bad_chains0 generated by sorting good_chains and then adding Gaussian noise. These values are highly autocorrelated (not independent) and not identically distributed since the flattened and sorted array is being reshaped into a 2D array. The third synthetic posterior called bad_chains1 is generated from good_chains by randomly introducing portions where consecutive samples are highly correlated to each other [1].

Effective Sample Size

When using MCMC sampling methods we may wonder if a particular sample is large enough to confidently compute quantities of interest like a mean or variance. This cannot be answered directly because samples from MCMC methods have some degree of autocorrelation. The amount of information contained in autocorrelated samples will be less than what we would get from an IID sample of the same size. The effective sample size (ESS) is an estimator that takes autocorrelation into account and provides the number of draws we would have if our sample was IID [1].

`good_chains`	`bad_chains0`	`bad_chains1`
4.389e+03	2.436	111.1

The ESS when the count of actual samples in our synthetic posteriors is 4000 [1].

Convergence of Markov chains is not uniform across parameter space. It is easier to get a good approximation from the bulk of a distribution than from the tails simply because the tails are dominated by rare events. As a general rule of thumb, an ESS greater than 400 is recommended, otherwise the estimation is unreliable [1].

Potential Scale Reduction Factor (R-hat)

We need ways to estimate convergence for finite samples. One pervasive idea is to run more than one chain starting from very different points and then check the resulting chains to see if they look similar to each other. This intuitive notion can be formalized into a numerical diagnostic known as $R^$ [1].

The overall idea is that the $R^$ for the parameter $θ$ is the standard deviation of all the samples of $θ$ , that is including all chains together, divided by the root mean square of the separated within chain standard deviations. Ideally we should get a value of 1, as the variance between chains should be the same as the variance within chains. From a practical point of view, values of $R^≲ 1.01$ are considered safe [1].

`good_chains`	`bad_chains0`	`bad_chains1`
1.000	2.408	1.033

The $R^$ from our synthetic posteriors [1].

Monte Carlo Standard Error

Using MCMC methods introduces an additional layer of uncertainty as we are approximating the posterior with a finite number of samples. We can estimate the amount of error introduced using the Monte Carlo Standard Error (MCSE), which is based on Markov chain central limit theorem. We should check the MCSE only once we are sure ESS is high enough and $R^$ is low enough, otherwise MCSE is of no use. As with the ESS the MCSE varies across the parameter space and we may want to evaluate it for different regions like specific quantiles [1].

`good_chains`	`bad_chains0`	`bad_chains1`
0.002381	0.1077	0.01781

The MCSE of our example posteriors [1].

Trace Plots

Trace plots are the most common plots in Bayesian literature. They are the first plot after inference for visual inspection. A trace plot is made by drawing the sampled values at each iteration step. From these plots were are able to see if different changes converged to the same distribution, getting a sense of the degree of autocorrelation [1].

Autocorrelation Plots

As see when discussing effective sample size, autocorrelation decreases the amount of information contained in a sample and thus is something we want to keep to a minimum. An autocorrelation plot lets us qualitatively inspect the results from ESS. The figure is the bar plot of the autocorrelation function, where the gray band represents the 95% confidence interval [1].

Rank Plots

Rank plots are histograms of the ranked samples. The ranks are computed by first combining all chains but then plotting the results separately for each chain. If all the chains are targeting the same distribution, we expect the ranks to have a uniform distribution. Additionally, if rank plots of all chains look similar, this indicates good mixing of the chains. Rank plots can be more sensitive than trace plots, are thus are recommended after refinement of sampling and modeling definitions. Trace plots are especially useful during early stages when we need to explore many different alternatives [1].

Divergences

The above diagnostics are summaries of the generated samples. Another way to perform a diagnostic is through monitoring the sampling method. One prominent example are the divergences present in Hamiltonian Monte Carlo (HMC) methods. Divergences can be visualized in trace plots as vertical bars at the bottom of the KDEs indicating where something went wrong during sampling [1].

References

Martin, Osvaldo A., Kumar, Ravin, Lao, Junpeng (2021) Bayesian Modeling and Computation in Python. Chapman and Hall/CRC.

Hamiltonian Monte Carlo (HMC)

Ethan Davis — Sat, 11 Jul 2026 03:55:25 +0000

Hamiltonian Monte Carlo

Hamiltonian Monte Carlo (HMC) is a type of MCMC method that makes use of gradients to generate new proposed states. HMC attempts to avoid the random walk behavior typical of Metropolis-Hastings by using the gradient to propose new positions far from the current one with high acceptance probability. This allows HMC to better scale to higher dimensions and in principle more complex geometries [1].

In physical systems, the Hamiltonian is a description of the total energy. The total energy can be decomposed into two terms, the kinetic and the potential energy. We can write the Hamiltonian of such a system as the following where $K (p, q)$ is the kinetic energy and $V (q)$ is the potential energy [1].

H (q, p) = K (p, q) + V (q) .

The probability of being at a particular position with a particular momentum is given as follows [1]. The exponential is a non-negative quantity and can therefore be viewed as an unnormalized probability distribution. The introduction of the minus sign in the exponential is simply a convention, meaning that higher values of energy correspond to lower values of probability [2].

p (q, p) = exp (- H (q, p)) .

To simulate such systems the Hamiltonian equations must be solve [1].

\frac{d q}{d t} \frac{d p}{d t} = \frac{\partial H}{\partial p} = \frac{\partial K}{\partial p} + \frac{\partial V}{\partial p} = - \frac{\partial H}{\partial q} = - \frac{\partial K}{\partial q} - \frac{\partial V}{\partial q} .

The potential energy $V (q)$ is given by the probability we wish to sample from $p (q)$ . For the momentum $K (p, q)$ we can invoke an auxiliary variable, and by choosing $p (p ∣ q)$ we can write $p (q, p) = p (p ∣ q) p (q)$ . This ensures we can recover the target distribution by marginalizing out the momentum. By introducing the auxiliary variable we can rewrite the equation [1].

H (q, p) = - lo g p (p ∣ q) - lo g p (q) .

We are free to choose the kinetic energy, and if we choose it to be Gaussian and drop the normalization constant, then we have the following where $M$ is the precision matrix [1].

K (p, q) = \frac{1}{2} p^{⊤} M^{- 1} p + lo g ∣ M ∣ .

If we choose $M$ to be the identity matrix we have $K (p, q) = \frac{1}{2} p^{⊤} p$ , making $\frac{\partial K}{\partial p} = p$ and $\frac{\partial K}{\partial q} = 0$ , and then we can simplify the Hamiltonian equations [1].

\frac{d q}{d t} \frac{d p}{d t} = p = - \frac{\partial V}{\partial q} .

Therefore, the HMC algorithm can be summarized as follows [1].

Sample $p \sim N (0, I)$
Simulate $q_{t}$ and $p_{t}$ for some amount of time $T$
$q_{T}$ is the new proposed state
Use the Metropolis acceptance criterion to accept or reject $q_{T}$ .

In an idealized model, energy is perfectly conserved. However, we use the Metropolis acceptance criterion to correct for errors introduced by the numerical simulation of the Hamiltonian equations. To compute the Hamiltonian equations we have to compute a trajectory. That is, all the intermediate points between one state and the next. In practice this involves computing a series of small integration steps. The most popular method is the leapfrog integrator. Leapfrog integration updates positions $q_{t}$ and momentums $p_{t}$ at interleaved time points, staggered in such that they leapfrog over each other [1].

An efficient HMC run requires proper tuning of its hyperparameters [1]:

The time discretization (step size of the leapfrog)
The integration time (number of leapfrog steps)
The precision matrix $M$ that parametrizes the kinetic energy.

For example, if the step size is too large then the leapfrog integrator will be inaccurate and too many proposals will be rejected. However, if it is too small we waste computation resources. If the number of steps is too small, the simulated trajectory at each iteration will be too short and sampling will fall back to random walk. But if too large the trajectory might run in circles and we again waste computation resources. If the estimated covariance (inverse of the precision matrix) is too different from the posterior covariance, the proposal momentum will be suboptimal and movement in the position space will be too large or too small in some dimension [1].

References

Martin, Osvaldo A., Kumar, Ravin, Lao, Junpeng (2021) Bayesian Modeling and Computation in Python. Chapman and Hall/CRC.
Bishop, Christopher M., Bishop, Hugh (2023) Deep Learning: Foundations and Concepts. Springer Cham.

Markov Chain Monte Carlo: Theoretical Foundations

Ethan Davis — Sat, 11 Jul 2026 03:52:53 +0000

Markov Chain Monte Carlo

Almost as soon as computers were invented, they were used for simulation. Markov chain Monte Carlo (MCMC) was invested as Los Alamos, Metropolis et al (1953) simulated a liquid in equilibrium with its gas phase. Their tour de force was the realization that they did not need to simulate the exact dynamics, they only needed to simulate some Markov chain with the same equilibrium distribution. The Metropolis algorithm was widely used by chemists and physicists, but was not widely known among statisticians until after 1990. Hastings (1970) generalized the Metropolis algorithm, and simulations following his scheme are said to use the Metropolis-Hastings (MH) algorithm [1].

A special case of the MH algorithm was introduced by Geman et al (1984) discussing optimization to find the posterior mode rather than simulation. Algorithms following their scheme are said to use the Gibbs sampler. It took some time for the spatial statistics community to understand that the Gibbs sampler simulated the posterior distribution, thus enabling full Bayesian inference of all kinds. Gelfand et al (1990) made the wider Bayesian community aware of the Gibbs sampler, and then it was rapidly realized that most Bayesian inference could be done using MCMC, whereas very little could be done without MCMC. Green (1995) generalized the MH algorithm as much as it could be generalized [1].

Theoretical Foundations

A sequence $X_{1}, X_{2}, \dots$ of random elements of some set is a Markov chain if the conditional distribution of $X_{n + 1}$ given $X_{1}, \dots, X_{n}$ depends on $X_{n}$ only. The set in which the $X_{i}$ take values is called the state space of the Markov chain. A Markov chain has stationary transition probabilities if the conditional distribution of $X_{n + 1}$ given $X_{n}$ does not depend on $n$ . This is the main kind of Markov chain of interest in MCMC. The joint distribution of a Markov chain is determined by the following [1].

The marginal distribution of $X_{1}$ called the initial distribution, and
The conditional distribution of $X_{n + 1}$ given $X_{n}$ called the transition probability distribution.

If the state space is finite or countable, written $x_{1}, \dots, x_{n}$ , then the initial distribution can be associated with a vector $λ = (λ_{1}, \dots, λ_{n})$ defined by $P (X_{1} = x_{i}) = λ_{i}$ for $i = 1, \dots, n$ , and the transition probabilities can be associated with a matrix $P$ having elements $p_{ij}$ defined by $P (X_{n + 1} = x_{j} ∣ X_{n} = x_{i})$ where $i = 1, \dots, n$ and $j = 1, \dots, n$ . When the state space is uncountable, we must think of the initial distribution and transition probability distribution as unconditional and conditional probability distributions [1].

A stochastic process is stationary if for every positive integer $k$ the distribution of the $k$ -tuple $(X_{n + 1}, \dots, X_{n + k})$ does not depend on $n$ . An initial distribution is said to be stationary, invariant, or equilibrium for some transition probability distribution if the Markov chain specified by this initial distribution and transition probability distribution is stationary. Stationarity implies stationary transition probabilities, but not vice versa. The Metropolis-Hastings-Green (MHG) algorithm constructs a transition probability distribution that preserves a specified equilibrium distribution [1].

A transition probability distribution is reversible with respect to an initial distribution if for its Markov chain $X_{1}, X_{2}, \dots$ , the distribution of pairs $(X_{i}, X_{i + 1})$ is exchangeable. Reversibility implies stationarity, but not vice versa. A reversible Markov chain has the same laws running forward and backward. That is, for any $i$ and $k$ , the distributions $(X_{i + 1}, \dots, X_{i + k})$ and $(X_{i + k}, \dots, X_{i + 1})$ are the same. All known methods for constructing transition probabilities that preserve a specified equilibrium are special cases of the MHG algorithm, and all elementary updates from the MHG algorithm are reversible [1].

A bit of compute code that makes a pseudorandom change to its state is an update mechanism. An update mechanism is elementary if it is not made up of parts that are themselves update mechanisms preserving the specified distribution. Suppose the specified distribution (the desired stationary distribution of the MCMC sampler) has unnormalied density $h$ . The Metropolis-Hastings update does the following [1]:

When the current state is $x$ , propose a move to $y$ , having conditional probability density given $x$ denoted $q (x, y)$ .
Calculate the Hastings ratio

$r (x, y) = \frac{h ( y ) q ( y , x )}{h ( x ) q ( x , y )} .$
Accept the proposed move $y$ with probability

$a (x, y) = min (1, r (x, y)),$ that is, the state after the update is $y$ with probability $a (x, y)$ , and the state after the update is $x$ with probability $1 - a (x, y)$ .

For example, consider the probability density function $h (x)$ is given by $h (x) = \frac{g ( x )}{\int _{- \infty \infty} g ( u ) d u}$ where we cannot or do not desire to solve the integral in the denominator analytically. Therefore, the distribution is only known up to some unknown constant: $h (x) \propto g (x)$ . Notice that the ratio $\frac{h ( y )}{h ( x )}$ from $r (x, y)$ does not depend on the normalizing constant. The latter term $\frac{q ( y , x )}{q ( x , y )}$ corrects for biases from the proposal distribution. The special case of the Metropolis-Hastings algorithm where the proposal distribution is symmetric meaning $q (y, x) = q (x, y)$ is referred to as the Metropolis algorithm [1].

The Metropolis-Hastings update is reversible with respect to $h$ , meaning that the transition probability that describes the update is an exact sampler of the specified distribution. If $X_{n}$ is the current state and $Y_{n}$ is the proposal, we have $X_{n} = X_{n + 1}$ whenever the proposal is rejected. The distribution of $(X_{n}, X_{n + 1})$ given rejection is exchangeable, and we must show that $(X_{n}, Y_{n})$ is exchangeable given acceptance. That is, we must show that for any function $f$ that has expectation, we can interchange the arguments of $f$ [1].

E [f (X_{n}, Y_{n}) a (X_{n}, Y_{n})] = \iint f (x, y) h (x) a (x, y) q (x, y) d x d y = E [f (Y_{n}, X_{n}) a (X_{n}, Y_{n})] .

This follows if we can interchange $x$ and $y$ in $h (x) a (x, y) q (x, y)$ . Only the set of $x$ and $y$ such that $h (x) > 0$ , $q (x, y) > 0$ , and $a (x, y) > 0$ contribute to the integral or sum in the discrete case. These inequalities further imply that $h (y) > 0$ and $q (y, x) > 0$ . Thus $r (y, x) = \frac{1}{r ( x , y )}$ for all $x$ and $y$ . Suppose that $r (x, y) \leq 1$ , so $r (x, y) = a (x, y)$ and $a (y, x) = 1$ . Then we have the following [1].

h (x) a (x, y) q (x, y) = h (x) r (x, y) q (x, y) = h (y) q (y, x) = h (y) a (y, x) q (y, x) .

Conversely, suppose that $r (x, y) > 1$ , so $a (x, y) = 1$ and $a (y, x) = r (y, x)$ , and we have the following [1].

h (x) a (x, y) q (x, y) = h (x) q (x, y) = h (y) r (y, x) q (y, x) = h (y) a (y, x) q (y, x) .

References

Brooks, Steve, Gelman, Andrew, Jones, Galin, Meng, Xiao-Li (Eds.) (2011) Handbook of Markov Chain Monte Carlo. Chapman and Hall/CRC.

EEG Feature Engineering and Selection

Ethan Davis — Fri, 10 Jul 2026 23:55:11 +0000

Feature Engineering & Selection

The event related features consist of chunks of time series concatenated from all the channels, resulting from a low pass or band pass filtering and a down sampling step. This category of features is relevant when considering evoked activity after the presentation of a stimulus (visual, auditory, or sensory). They are therefore of interest when one is expecting significant changes in several amplitudes occurring at a given moment [1].

The spectral features are used in the case of detection for an oscillatory activity, when changes in EEG rhythms amplitudes are expected. The features are associated with the power spectra estimated in a given channel and in a given frequency band for a specific time window. Power spectra can be computed through many methods like, most notably, spectrogram, Morlet wavelet scalogram, and autoregressive models [1].

Spatial filtering can be a valuable tool both for event related and spectral features. It relies on the combination of signals, recorded from different sensors, to obtain a new one, associated with an improved signal-to-noise ratio. Spatial filtering methods can be divided into three categories [1].

First, and not data driven, this category relies on physical considerations regarding the way the signals propagate through different brain tissues. The most famous illustration of this category is the Laplacian filter. In its simplest version, the small Laplacian consists of a derivation of the EEG waveform from the average signal computed with the four nearest neighbors for each electrode location [1].

Second, and data driven and unsupervised, this spatial filtering category can rely on a principal component analysis (PCA) approach. Third, and data driven and supervised, has its most famous examples in common spatial patterns (CSP) for spectral patterns and xDawn for event related features. CSP consists of a linear combination of EEG signals to maximize the difference between two classes in terms of variance. The xDawn approach aims at improving the signal-to-noise ratio through a projection of the raw EEG signals onto an estimated subspace of the evoked potentials [1].

Even though spectral and event related features are the most used in EEG literature, alternative features have been considered in past years. Firstly, features relying on covariance matrices have recently been used extensively, in particular for Riemannian geometry based classification. Despite an unclear neurophysiological interpretation, they have reached state of the art performance and won a large number of competitions. Secondly, new features, which take into account the interconnected nature of brain functioning and intensity of interactions, have emerged [1].

Estimating brain interactions in real time is not trivial: it consists of finding a compromise between ensuring the quasi-stationary nature of the signals and the statistical reliability of the functional connectivity estimation. Recent studies have considered the use of brain network metrics as potential features. These include local and global network metrics [1].

At the local scale, the node degree counts the number of connections linking one node to another. In weighted networks, it is referred to as node strength and consists of summing the weights of the connections of the considered node. Another local scale property of interest is the betweenness centrality defined as the extent to which a node lies between other pairs of nodes. This metric enables the identification of the nodes that are crucial for information transfer between distant regions [1].

At the global scale, there is the characteristic path length and clustering coefficient. The characteristic path length indicates the global tendency of the nodes in the network to integrate and exchange information. The clustering coefficient measures the tendency of having nodes whose neighbors are mutually interconnected [1].

The feature selection is a crucial step as it prevents redundancy, ensures the reliability of the features, reduces dimensionality tuned, and helps to provide interpretable results. It can be divided into three categories: embedded, filtered, and wrapper methods [1].

In filter methods, the feature selection is performed independently and before evaluation. For example, the $R^{2}$ score criterion can be used to assess to which extent a given feature is influenced by a task performed by the subject. In wrapper methods, the feature selection utilizes the classification. In other words, in an iterative process, the relevance of each subset of features is assessed through classification performance until a given criterion is met. The embedded method consists of integrating both the feature selection and classification in the same process, for example, with a decision tree [1].

References

Corsi, Marie-Constance (2023) Electroencephalography and Magnetoencephalography. Springer US.

Cleaning EEG: Artifacts, Preprocessing, and Source Reconstruction

Ethan Davis — Fri, 10 Jul 2026 23:48:28 +0000

Data Analysis

Artifacts are signals that make recording more difficult and may hamper the analysis of the brain activity recorded with EEG. Such artifacts can be divided into two categories: neurophysiological artifacts and environmental noise [1].

Neurophysiological artifacts correspond to noise generated by the subjects themselves, whether it is voluntary or not. The brain is far from being the only organ that generates electromagnetic activity. In particular, the eyes and heart produce electromagnetic activity. As a result, the main neurophysiological artifacts are related to cardiac activity and ocular activity (blinks and saccades), and can be visually seen during an EEG recording. A possible way to reduce ocular artifacts is to instruct the subject to avoid moving their eyes and, for short recordings only, to avoid eye blinking [1].

Another neurophysiological artifact, can be induced by the subject’s voluntary motion. Motion engenders muscular activity that can distort the recorded brain signals. Typical examples are jaw clenching and swallowing. A possible way to reduce the muscular artifacts is to instruct the subjects to remain as quiet as possible and to avoid moving their jaws [1].

Environmental noise refers to artifacts generated by the environment that surrounds the experimental setup. They can be linked with mechanical vibrations like the presence of a tramway nearby, or associated with power lines occurring at 50 or 60 Hz. Use of a Faraday cage (shielded room) may help to achieve the absolute best signal-to-noise ratio. Modern EEG equipment, especially systems with active electrodes, can also reduce noise. Additionally, system noise refers to artifacts generated by the sensors themselves that could be broken [1].

Preprocessing for artifact removal is probably the most crucial part when analyzing EEG data. The point here is to remove noise without eliminating information of interest of distorting the signal. The first thing to do is to extensively study the dataset, in particular, inspecting the EEG signals but also the associated broadband power spectra. This preliminary step enables identification of most artifacts and, most importantly, if they have a specific temporal and frequency signature (presence of periodic artifacts) [1].

From this point, it is possible to choose a specific strategy to remove the observed noise. In the case of cardiac and ocular artifacts, given their clear pattern, an efficient way to isolate and reduce them consists of applying independent component analysis (ICA). In the case of artifacts at a specific frequency, like power like noise at 50 or 60 Hz, one can consider applying notch filters. In the case of muscular activity, applying a low pass filter with a cutoff frequency at 40 Hz can be of interest [1].

It is possible to directly analyze the signals recorded by the sensors. In such a case, one will say that the analysis is performed in the space of the sensors. However, it is possible to go one step further and estimate the activity within the brain. This processing step is called source reconstruction and consists of estimating neural correlations of signals locations. It can be performed when one wants to have access to a higher spatial resolution and to provide a more accurate description, and interpretation, of the neurophysiological activity occurring. For that purpose, both direct and inverse problems need to be solved [1].

For the direct problem, we aim at modeling the electromagnetic field produced. For that purpose, it is necessary to consider both a physical model and a model that will predict how the electromagnetic field will be generated at the scalp. The simplest model is the spherical model, which considers the head as an ensemble of spheres (brain, cerebrospinal fluid, skull, scalp) each characterized by a given conductivity. More realistic models rely on geometrical reconstruction of the different layers that form the head tissues, directly extracted from magnetic resonance imaging (MRI). They consist of building meshes of the interfaces between different tissues [1].

The purpose of the inverse problem is to reconstruct brain activity in the physical model. One of the main challenges is the nonuniqueness of its solution. In other words, a large number of brain activity patterns could generate the same signature detected at the sensor level. Therefore, some constraints or assumptions are essential to lead to a unique solution that best reflects the acquired data. Dipole modeling methods rely on modeling with a reduced number of equivalent dipoles and, therefore, are based on an a priori hypothesis of the required number of sources. Scanning methods consist of estimating the probability of presence of a current dipole inside each voxel [1].

References

Corsi, Marie-Constance (2023) Electroencephalography and Magnetoencephalography. Springer US.

Recording EEG: Electrodes, Montages, and Protocol

Ethan Davis — Fri, 10 Jul 2026 23:45:37 +0000

Experiment Protocol

EEG signals are recorded through the use of electrodes placed over the scalp. The EEG relies on different types of electrodes: wet/dry and active/passive electrodes. Wet electrodes need an electrolytic gel to enable conduction between the skin and electrode. Dry electrodes behave as a conductor between the skin and electrode. Active electrodes contain an electronic module that performs a pre-amplification of the signal to ensure stability of the system from changes in impedance and noise. Passive electrodes do not use a pre-amplification module [1].

Even though some differences can be found from one EEG device to another, there are standardized ways to name and localize EEG sensors (also called channels). Each channel is referred to by a letter and number. Odd channels are located on the left hemisphere and even ones on the right hemisphere. The letters correspond to the area: frontal, temporal, parietal, central, and occipital. In addition to the sensors, one can also find landmarks: nasion, inion, and pre-auricular points [1].

Depending on the scientific question to be addressed and, therefore, the brain areas of interest, different montages can be found. One can build an EEG montage from less than 5 electrodes and up to 256 electrodes. EEG measurements rely on a difference of electrical potentials. For this purpose, two montages can be considered: the referential and bipolar montage [1].

In the referential montage each difference of electrical potentials considers an electrode placed over the scalp and a reference. The choice of the reference is crucial. The most commonly chosen locations for the reference electrodes are the mastoids (temporal bone behind the ears), though several studies prefer placing the reference at the vertex (Cz, midline of the scalp). The bipolar montage consists of performing the difference between two electrodes after the experiment. Another electrode referred to as the ground electrode is used. Among the chosen locations is the scapula (shoulder blade) [1].

There has been increased interest in developing wearable EEG, to remove wires and reduce its dimension but also to enable long lasting recordings in a less constrained environment. Three bottlenecks need to be overcome: the EEG electrodes, hard to put on and keep in place on the head; the EEG hardware, to make it less power consuming and miniaturized; and the EEG software, to propose the most intelligible and reliable information regarding the captured brain activity [1].

For data acquisition, it will consist of cleaning the locations where electrodes will be in contact with the skin, like forehead and mastoids. The electrodes and EEG cap are then placed. Several key distances can be measured to verify the cap is well placed. Then, the experimenter needs to ensure that the communication between electrodes and scalp is established. For that purpose, an assessment of the impedance is made for each electrode. The lower it is, the better it is. Once impedance are lower than a certain threshold, typically a few kOhms, then the experiment can start [1].

Once the subject is correctly in place, the experimenter can start some pre-recordings to check quality of the signal and give specific instructions to the subjects accordingly, like loosening the jaw to avoid muscular artifacts. Finally, the experimenter can give further instructions regarding the task to perform before starting the recordings. After the session ends, the data is stored in specific servers to be processed [1].

References

Corsi, Marie-Constance (2023) Electroencephalography and Magnetoencephalography. Springer US.

EEG Basics: How Brain Activity Becomes a Signal

Ethan Davis — Fri, 10 Jul 2026 23:43:59 +0000

Basic Principles

The origin of electroencephalogram (EEG) signals comes from neurons creating electrical signals that are transmitted to other cells through synapses. An action potential (AP) arrives at a synaptic cleft (step 1 in the figure) where it transmits chemical information through neurotransmitters (step 2) that generate postsynaptic potentials (PSPs) and local current (step 3). A PSP will create a current sink and will propogate to the cell body to generate a current source (step 4) [1].

As a result, the postsynaptic potential creates an electrical dipole consisting of a negative pole (the sink) and a positive pole (the source). The dipole generates primary (intracellular) currents and secondary (extracellular) currents. EEG signals result from postsynaptic potentials. More specifically, EEG signals result from the spatial and temporal summation of activity from a large population of synchronous neurons [1].

EEG signals correspond to a difference between electrical potentials, mainly due to extracellular currents. It is sensitive to both radial currents (actively located at the gyrus level) and to tangential currents (generated within sulci) even though it has strong sensitivity to radial currents. EEG, though, is also strongly attenuated and deformed by crossing through the skull [1].

Characteristics of EEG [1]

Property	Value
Measurement	Difference of potentials, extracellular currents
Spatial resolution	2–3 cm
Temporal resolution	1 ms or less
Amplitudes	≈ 100 μVolts
Advantages	Portable, low cost
Drawbacks	Need of a reference, affected by bone, diffuse

There are two main types of electrophysiological activity of interest that are exploited for EEG: the evoked and oscillatory activity. Evoked responses are weak variations of electromagnetic activity resulting from a stimulation, such as in response to a task performed by the participant. There is a specific way to identity and describe evoked responses according to their latency, amplitude, shape, and polarity [1].

The figure first shows a positive deflection occurring 300 ms after the presentation of a stimulation. These waves of evoked activity can reflect different mechanisms: the early components are mostly exogenous and are related to stimulus characteristics, and late components are endogenous and are related to the performed task and to the subject’s state. Alternatively, oscillatory activity, or induced activity, results from the summation of activity in a given brain region. These rhythms are mainly defined by their frequency, amplitude, shape, location, and duration [1].

Each frequency band is referred to by a Greek letter and corresponds to a subject’s state. Delta (0.5–3 Hz) and Theta (3–7 Hz) rhythms are detected in deep and slight sleeps respectively. Alpha (8–12 Hz in posterior areas) and Mu (7–13 Hz in central areas) rhythms are observed in quiet watch and resting state (with eyes closed for Alpha). Beta (13–30 Hz) rhythm is detected during active watch and cognitive tasks such as motor imagery. Gamma rhythm (divided into two sub-rhythms: slow in 30–70 Hz and fast beyond 70 Hz) is observed during specific cognitive processing [1].

References

Corsi, Marie-Constance (2023) Electroencephalography and Magnetoencephalography. Springer US.

The Tangent Space of the SPD Manifold for EEG Classification

Ethan Davis — Fri, 10 Jul 2026 22:58:44 +0000

Tangent Space

Current state-of-the-art (SOTA) machine learning (ML) for electroencephalogram (EEG) uses across channel covariance matrices and Riemannian geometric statistics on the symmetric positive definite (SPD) manifold for classification [1, 2]. These methods either discriminate directly on the SPD manifold or in the vector space that is tangent to the manifold. In this section we analyze the latter method of discrimination that performs classification in the tangent space (TS) of the SPD manifold.

Definitions

A manifold is a collection of points that locally, but not globally, resembles Euclidean space. A Riemannian metric is defined by a smoothly varying collection of scalar products $⟨ \cdot, \cdot ⟩_{x}$ in each tangent space $T_{x} M$ at points $x$ on a manifold $M$ . For example, the inner product gives a norm $∣ \cdot ∣_{x} : T_{x} M \to R$ by $∣ v ∣_{x 2} = ⟨ v, v ⟩_{x}$ . The shortest distance between two points on a manifold is a geodesic $γ (t)$ , and it is computed by integrating the norm along its curve [3].

We can map vectors in the tangent space to the manifold using geodesics. The vector $v \in T_{x} M$ can be mapped to the point of the manifold that is reached after a unit time $t = 1$ by the geodesic $γ (t)$ starting at $x$ with initial velocity $γ^{'} (0) = v$ . This mapping $Exp_{x} : T_{x} M \to M$ is called the exponential map where $Exp_{x} (v) = γ (1)$ . The inverse of the exponential map is the logarithmic map $Log_{x} (y)$ and is the smallest vector $v$ as measured by the Riemannian metric so that $Exp_{x} (v) = y$ [3].

In general, finding geodesics involves solving second-order ordinary differential equations (ODEs) [3]. However, SPD matrices have additional Lie group structure that can be used to simplify algorithms and speed up computations of Riemannian metrics and geodesics [4]. A group is a set $G$ with an associative binary operation satisfying $(x y) z = x (yz)$ for all $x, y, z \in G$ , contains an identity element $e \in G$ so that $e x = x e = x$ for all $x \in G$ , and contains an inverse element $x^{- 1} \in G$ for each $x \in G$ where $x^{- 1} x = x x^{- 1} = e$ . A Lie group is a smooth manifold and also a group [3].

For each $x, y$ in a Lie group $G$ , left-translations by $y$ are denoted $L_{y} : x \mapsto y x$ . The differential $(L_{y})_{*}$ of left-translation maps the tangent space $T_{x} G$ to the tangent space $T_{y x} G$ . In particular, $(L_{y})_{*}$ maps any vector $u \in T_{e} G$ to the vector $(L_{y})_{*} u \in T_{y} G$ . The vector field $u ~ (y) = (L_{y})_{*} u$ is called left-invariant since it is invariant under left-translations $u ~ \circ L_{y} = (L_{y})_{*} u ~$ for all $y \in G$ [3]. In other words, pullback to the identity and pushforward to another object are commutative under left-translation.

(L_{y})_{*} u ~ (x) = (L_{y})_{*} (L_{x})_{*} u = (L_{y x})_{*} u = u ~ (y x) = u ~ (L_{y} x) = (u ~ \circ L_{y}) x .

The left-translation maps give a useful way of defining Riemannian metrics on Lie groups. By definition the tangent space $T_{e} G$ of a Lie group $G$ at the identity element, typically denoted $g$ , is called a Lie algebra. Given an inner product $⟨ \cdot, \cdot ⟩_{g}$ on the Lie algebra, we can extend it to an inner product on tangent spaces at all elements of the group by setting $⟨ u, v ⟩_{g} = ⟨(L_{g^{- 1}})_{*} u, (L_{g^{- 1}})_{*} v ⟩_{g}$ . This defines a left-invariant Riemannian metric on $G$ since $⟨(L_{h})_{*} u, (L_{h})_{*} v ⟩_{h g} = ⟨ u, v ⟩_{g}$ for any $u, v \in T_{g} G$ [3].

⟨(L_{h})_{*} u, (L_{h})_{*} v ⟩_{h g} = ⟨(L_{(h g)^{- 1}})_{*} (L_{h})_{*} u, (L_{(h g)^{- 1}})_{*} (L_{h})_{*} v ⟩_{g} = ⟨(L_{(h g)^{- 1} h})_{*} u, (L_{(h g)^{- 1} h})_{*} v ⟩_{g} = ⟨(L_{g^{- 1} h^{- 1} h})_{*} u, (L_{g^{- 1} h^{- 1} h})_{*} v ⟩_{g} = ⟨(L_{g^{- 1}})_{*} u, (L_{g^{- 1}})_{*} v ⟩_{g} = ⟨ u, v ⟩_{g} .

For a matrix $P \in R^{n \times n}$ and any non-zero vector $x \in R^{n}$ , if $P = P^{⊤}$ and $x^{⊤} P x > 0$ , then $P$ is an SPD matrix [5]. The space of SPD matrices is a smooth manifold, and not a vector space since it lacks an additive identity, additive inverses, and zero and negative real scalar multiplication [4]. The SPD manifold $Sym_{n +}$ of all SPD matrices is a Lie group with matrix multiplication as its group operation [3]. Therefore, we can derive invariant metrics and geodesics on the SPD manifold.

Let $A$ be a matrix from the general linear group $GL (n)$ of non-singular matrices. Then $A P A^{⊤} \in Sym_{n +}$ since $(A P A^{⊤})^{⊤} = (A^{⊤})^{⊤} P A^{⊤} = A P A^{⊤}$ and for any non-zero vector $x \in R^{n}$ then $x^{⊤} (A P A^{⊤}) x = (x^{⊤} A) P (A^{⊤} x) = (A^{⊤} x)^{⊤} P (A^{⊤} x) > 0$ . Given two matrices $P, Q \in Sym_{n +}$ we can derive the distance between them as a norm from the identity by choosing $A = P^{- 1/2}$ [4]. That is, with left-invariance we can pullback to the identity before pushforward.

dist (P, Q) = dist (Id, P^{- 1/2} Q P^{- 1/2}) = N (P^{- 1/2} Q P^{- 1/2}) .

By left-invariance at the identity of a Lie group we have the vector field $u ~ (x) = (L_{x})_{*} u$ . The geodesic starting at the identity $e$ with initial velocity $u$ satisfies $x (0) = e$ and $x^{'} (t) = u ~ (x (t))$ . For matrix groups, this ODE becomes $x^{'} (t) = x (t) u$ and is uniquely solved by the matrix exponential $x (t) = exp (t u)$ . In other words, the exponential map is the matrix exponential and its inverse the logarithmic map is the matrix logarithm [3].

From the figure take $R = P^{- 1/2} Q P^{- 1/2}$ . The matrix logarithm $lo g (R)$ produces the vector for the norm $N (R)$ . This norm is given by the Frobenius norm $∣ \cdot ∣_{F}$ . Furthermore, since $R$ is SPD, it has eigendecomposition $R = V Λ V^{⊤}$ where $Λ = diag (λ_{1}, \dots, λ_{n})$ and $λ_{i} > 0$ . By the matrix logarithm $lo g (R) = V diag (lo g λ_{1}, \dots, lo g λ_{n}) V^{⊤}$ . Therefore, by the orthogonal invariance of the Frobenius norm, the norm $N (R)$ can be written as the square root of the sum of squared eigenvalue logarithms [4].

dist (P, Q) = N (P^{- 1/2} Q P^{- 1/2}) = ∣ lo g (R) ∣_{F} = (i = 1 \sum n lo g^{2} λ_{i})^{1/2} .

the figure gives a closed form solution for the distance between two SPD matrices on $Sym_{n +}$ . Not only is it exact and so does not require optimization, but it can be computed in the vector space of the tangent space and so does not require integration on the curve. Furthermore, by the nature of logarithms we see that matrices with zero or negative eigenvalue are in fact infinite distance from SPD matrices on $Sym_{n +}$ , contrary to Euclidean space and metrics [4].

We can also find closed form solutions for the logarithmic and exponential maps on $Sym_{n +}$ . Rather than solve for the norm of the logarithmic map, we can solve for the vector itself. This is done by pushing forward to $P$ after we pulled back to the identity from $P$ . Furthermore, as the inverse of the logarithmic map, the closed form solution for the exponential map takes a vector as input and outputs the point on $Sym_{n +}$ reached after unit time elapsed along the geodesic.

Log_{P} (Q) = P^{1/2} lo g (P^{- 1/2} Q P^{- 1/2}) P^{1/2} .

Exp_{P} (W) = P^{1/2} exp (P^{- 1/2} W P^{- 1/2}) P^{1/2} .

At the foundation of statistics is the notion of a distance. The closed form solution of the mean in Euclidean space $\frac{1}{N} \sum_{i = 1 N} X_{i}$ relies on the additive structure of vector space and does not generalize to Riemannian space. However, there are defining properties of the mean that do generalize: The geometric mean is a least-squares centroid that minimizes the sum-of-squared distances. A natural strategy for computing the geometric mean is gradient descent optimization [6].

y ˉ = y \in M arg min i = 1 \sum N dist (y, y_{i})^{2} .

The exponential and logarithmic map, distance, and mean on $Sym_{n +}$ are common Riemannian geometric statistics used by SOTA ML for EEG that takes across channel sample covariance matrices for classification [1, 2]. The minimum distance to the mean (MDM) classifier is trained by memorizing the mean of each class and is tested by predicting classes based on the minimum distance to the memorized means [1]. All operations occur on the SPD manifold. In the next section we describe in more detail an alternative algorithm that performs classification in the tangent space of $Sym_{n +}$ .

Example

Consider two three-dimensional tensors of EEG recordings $X_{1}, X_{2} \in R^{p \times q \times r}$ from separate classes where the dimensions $p, q, r$ are the trials, channels, and time, respectively. For each trial, let sample covariance matrices (SCMs) be computed between channels across time so that $Σ_{1 (k)}, Σ_{2 (k)} \in R^{q \times q}$ for $k = 1, \dots, p$ . Assume each SCM is conditioned so that it is SPD. Furthermore, let us define $Σ ˉ_{1}, Σ ˉ_{2} \in R^{q \times q}$ as geometric mean SCMs averaged across all trials $k$ of each class.

With the geometric mean of each class $Σ ˉ_{1}, Σ ˉ_{2} \in Sym_{n +}$ we define $Σ ˉ \in Sym_{n +}$ as the geometric mean of these means. Then, using the logarithmic map we find the vectors between all points in all clusters to $Σ ˉ$ . We represent these vectors in the tangent space $T_{Σ ˉ} Sym_{n +}$ of the mean point $Σ ˉ$ . This linearization of SPD matrices is done in a way that respects the metric space of $Sym_{n +}$ . Once in the tangent space $T_{Σ ˉ} Sym_{n +}$ we can train standard ML algorithms that are optimized for vector space.

Once in the tangent space $T_{Σ ˉ} Sym_{n +}$ we can train standard ML algorithms that are optimized for vector space. This method of linearization from the space of SPD matrices to vector space is the standard used by software libraries like PyRiemann [7]. As of now, SOTA EEG classifiers are those that learn and predict in the tangent space of the SPD manifold rather than on the manifold itself [1, 2]. An explanation for this is that machine learning is traditionally optimized for vector space.

References

Barachant, Alexandre, Bonnet, Stéphane, Congedo, Marco, Jutten, Christian (2012) Multiclass Brain–Computer Interface Classification by Riemannian Geometry. IEEE Transactions on Biomedical Engineering.
Sylvain Chevallier, Emmanuel K. Kalunga, Quentin Barthélemy, Florian Yger (2018) Riemannian Classification for SSVEP-Based BCI: Offline versus Online Implementations. CRC Press.
Stefan Sommer, Tom Fletcher, Xavier Pennec (2020) Introduction to differential and Riemannian geometry. Academic Press.
Xavier Pennec (2020) Manifold-valued image processing with SPD matrices. Academic Press.
Kevin P. Murphy (2022) Probabilistic Machine Learning: An Introduction. MIT Press.
Tom Fletcher (2020) Statistics on manifolds. Academic Press.
Alexandre Barachant, Quentin Barthélemy, Jean-Rémi King, Alexandre Gramfort, Sylvain Chevallier, Pedro L. C. Rodrigues, Emanuele Olivetti, Vladislav Goncharenko, Gabriel Wagner vom Berg, Ghiles Reguig, Arthur Lebeurrier, Erik Bjäreholt, Maria Sayu Yamamoto, Pierre Clisson, Marie-Constance Corsi, Igor Carrara, Apolline Mellot, Bruna Junqueira Lopes, Brent Gaisford, Ammar Mian, Anton Andreev, Gregoire Cattan, Arthur Lebeurrier (2025) pyRiemann. Zenodo.

Common Spatial Pattern (CSP): A Correctness Proof

Ethan Davis — Fri, 10 Jul 2026 22:58:37 +0000

Common Spatial Pattern

Correctness Proof

The common spatial pattern (CSP) allows one to maximize the variance of signals from one condition and at the same time minimize the variance of signals from another condition [1]. For example, consider two three-dimensional tensors ( $trials \times channels \times time$ ) of electroencephalogram (EEG) from separate classes $X_{1}, X_{2} \in R^{p \times q \times r}$ . For each trial let us assume zero mean and compute between channel across time sample covariance matrices (SCMs) $Σ_{i (k)} = \frac{1}{r - 1} X_{i (k)} X_{i (k) ⊤} \in R^{q \times q}$ for $i = 1, 2$ and $k = 1, \dots, p$ . Furthermore, let us take $Σ ˉ_{i} \in R^{q \times q}$ as the arithmetic mean SCM for each class across all trials $k$ .

Our example is a standard data and preprocessing pipeline for CSP used by software libraries like PyRiemann [2]. We assume the SCMs are conditioned so that they are symmetric positive definite (SPD). By definition, for a matrix $A \in R^{n \times n}$ and any non-zero vector $x \in R^{n}$ , if $A = A^{⊤}$ and $x^{⊤} A x > 0$ , then $A$ is SPD. As a result, all eigenvalues are positive $λ_{i} > 0$ from the eigendecomposition $A = Q D Q^{⊤}$ where $D = diag (λ_{1}, \dots, λ_{n})$ [3]. Furthermore, the sum of two SPD matrices $A, B \in R^{n \times n}$ is also SPD since $x^{⊤} (A + B) x = x^{⊤} A x + x^{⊤} B x > 0$ . Let us diagonalize the mean SCMs from our example.

Σ ˉ_{1} + Σ ˉ_{2} = V Λ V^{⊤} .

The whitening transformation of an SCM results in zero mean, unit variance, and zero covariance [1]. In other words, the whitening transformation maps an SCM to the identity matrix $Σ \mapsto I$ . We define the whitening matrix as $W = V Λ^{- 1/2} V^{⊤} \in R^{n \times n}$ . The result is no longer an ellipse when plotted but rather the unit circle.

W (Σ ˉ_{1} + Σ ˉ_{2}) W^{⊤} = V Λ^{- 1/2} V^{⊤} (V Λ V^{⊤}) V Λ^{- 1/2} V^{⊤} = I .

We can rewrite our whitening transformation to show that the variance of the signals in our first class $Σ ˉ_{1}$ is maximized, while the variance of the signals in our second class $Σ ˉ_{2}$ is also minimized. We do this by distributing the whitening matrix throughout the sum of our SCMs. The result is two transformed SCMs we denote by $Σ ˉ_{1'}, Σ ˉ_{2'} \in R^{n \times n}$ whose sum is the identity matrix.

W (Σ ˉ_{1} + Σ ˉ_{2}) W^{⊤} = W Σ ˉ_{1} W^{⊤} + W Σ ˉ_{2} W^{⊤} = Σ ˉ_{1'} + Σ ˉ_{2'} = I .

Let us use $V_{i'} Λ_{i'} V_{i'⊤}$ where $i = 1, 2$ to denote the diagonalization of our transformed SCMs. Since $Σ ˉ_{1'} + Σ ˉ_{2'} = I$ the sum of the diagonals $d_{ij} = diag (Λ_{i'})_{j}$ is $d_{1 j} + d_{2 j} = 1$ . Furthermore, $d_{ij} > 0$ since $Σ ˉ_{i'}$ is SPD. Therefore, the diagonals are bounded between 0 and 1. When the variance $d_{aj} ≫ d_{bj}$ for $a \neq = b$ , the classification problem is more or less solved for unseen EEG trials. However, when $d_{aj}$ is close to $d_{bj}$ the discrimination between classes is more ambiguous [1].

References

Benjamin Blankertz (2018) Gentle Introduction to Signal Processing and Classification for Single-Trial EEG Analysis. CRC Press.
Alexandre Barachant, Quentin Barthélemy, Jean-Rémi King, Alexandre Gramfort, Sylvain Chevallier, Pedro L. C. Rodrigues, Emanuele Olivetti, Vladislav Goncharenko, Gabriel Wagner vom Berg, Ghiles Reguig, Arthur Lebeurrier, Erik Bjäreholt, Maria Sayu Yamamoto, Pierre Clisson, Marie-Constance Corsi, Igor Carrara, Apolline Mellot, Bruna Junqueira Lopes, Brent Gaisford, Ammar Mian, Anton Andreev, Gregoire Cattan, Arthur Lebeurrier (2025) pyRiemann. Zenodo.
Kevin P. Murphy (2022) Probabilistic Machine Learning: An Introduction. MIT Press.