Primer: independent component analysis

On By @thinquanautIn Primers

Independent component analysis, ICA, solves a problem that sounds impossible. Several signals were mixed together before you recorded them, you have only the mixtures, and you want the originals back. It is the cocktail party: two people talking, two microphones, each microphone hearing both. In neuroscience the “people” are an eye blink, a heartbeat, a muscle, a patch of cortex, and the “microphones” are electrodes on the scalp, all hearing everything. This primer shows how ICA does it, writes the algorithm in twenty lines, and uses it to pull the eye blinks out of a real EEG recording. It assumes lesson 4, because ICA starts where PCA stops.

The problem, in two channels

Start with a case where we know the answer. Two sources, a sine wave and a square wave. Two sensors, each recording a different blend of the two. Nobody tells the algorithm the blend.

import numpy as np
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
n = 2000
time = np.linspace(0, 8, n)
s1 = np.sin(2 * np.pi * time)                              # a smooth rhythm
s2 = np.sign(np.sin(2 * np.pi * 1.7 * time + 1))           # a square wave at a different rate
S = np.vstack([s1, s2])                                    # 2 sources x n samples
A = np.array([[1.0, 0.6],                                  # how much of each source reaches each sensor
              [0.4, 1.0]])
X = A @ S + 0.05 * rng.standard_normal(S.shape)            # 2 sensors x n samples: mixtures plus a little noise
print("sensor 1 = 1.0 * s1 + 0.6 * s2;  sensor 2 = 0.4 * s1 + 1.0 * s2")
# sensor 1 = 1.0 * s1 + 0.6 * s2;  sensor 2 = 0.4 * s1 + 1.0 * s2

Why PCA cannot do it, and what can

The obvious tool is PCA from lesson 4: find the directions the data vary along most. But PCA has a built-in assumption that gives the wrong answer here. Its components are always at right angles to each other, and the directions the sources mix along are not at right angles unless you are lucky. PCA will decorrelate the sensors, which removes the linear dependence between them, but a sine and a square wave that are uncorrelated can still be thoroughly mixed.

The two sensor readings plotted against each other. The cloud is two parallel bands: the square wave jumps between them, the sine wave slides along them. Left: PCA's axes, forced to be perpendicular, cut across both bands. Right: ICA's axes run along the bands and between them, which are the directions the two sources were actually added along.
The two sensor readings plotted against each other. The cloud is two parallel bands: the square wave jumps between them, the sine wave slides along them. Left: PCA’s axes, forced to be perpendicular, cut across both bands. Right: ICA’s axes run along the bands and between them, which are the directions the two sources were actually added along.

ICA uses a different clue: independence, which is stronger than uncorrelatedness, and a fact about sums. When you add independent signals together the result is always closer to a bell curve than the ingredients were; this is the central limit theorem doing what it does. So a mixture is more Gaussian than its sources, and the way to unmix is to search for the directions along which the data look least Gaussian. A sine wave and a square wave are both very un-bell-shaped, and in the right directions you see that; in the wrong directions you see a blur.

FastICA, the standard algorithm, does exactly this. First it whitens the data with PCA, which handles the scaling and leaves only a rotation to find. Then it rotates, step by step, to maximise a measure of non-Gaussianity, keeping the rows orthogonal so that it finds different sources rather than the same one twice. The non-Gaussianity measure is built from tanh, for reasons that are in the paper and do not matter here.

def whiten(X, k):
    """PCA-whiten: rotate to the top k components and scale each to unit variance."""
    Xc = X - X.mean(axis=1, keepdims=True)
    cov = Xc @ Xc.T / Xc.shape[1]
    variance, vectors = np.linalg.eigh(cov)
    order = np.argsort(variance)[::-1][:k]
    variance, vectors = variance[order], vectors[:, order]
    W = (vectors / np.sqrt(variance)).T                    # k x channels
    unwhiten = vectors * np.sqrt(variance)                 # channels x k: undoes W
    return W @ Xc, unwhiten

def fastica(Z, seed=0, steps=200):
    """Find a rotation of whitened data Z that makes the rows as non-Gaussian as possible."""
    k = Z.shape[0]
    W = np.random.default_rng(seed).standard_normal((k, k))
    for _ in range(steps):
        Y = W @ Z
        g, g_prime = np.tanh(Y), 1 - np.tanh(Y) ** 2       # the non-linearity and its derivative
        W_new = g @ Z.T / Z.shape[1] - g_prime.mean(axis=1)[:, None] * W
        u, _, vt = np.linalg.svd(W_new)
        W_new = u @ vt                                     # keep the rows orthogonal (decorrelated)
        if np.max(np.abs(np.abs((W_new * W).sum(axis=1)) - 1)) < 1e-6:
            break                                          # nothing changed: converged
        W = W_new
    return W

Z, unwhiten = whiten(X, 2)
W = fastica(Z)
recovered = W @ Z
for i in range(2):
    print(f"recovered {i}: correlation with s1 {abs(np.corrcoef(recovered[i], s1)[0, 1]):.2f}, with s2 {abs(np.corrcoef(recovered[i], s2)[0, 1]):.2f}")
# recovered 0: correlation with s1 0.99, with s2 0.00
# recovered 1: correlation with s1 0.03, with s2 1.00
Top: the two sources, which the algorithm never sees. Middle: what the two sensors record. Bottom: what ICA recovers from the sensors alone. The recovered signals match the originals to a correlation of 0.99 and 1.00, up to two things ICA cannot know: their order and their sign.
Top: the two sources, which the algorithm never sees. Middle: what the two sensors record. Bottom: what ICA recovers from the sensors alone. The recovered signals match the originals to a correlation of 0.99 and 1.00, up to two things ICA cannot know: their order and their sign.

That last sentence is the first thing to remember about ICA in practice: the components come out in no particular order, scaled arbitrarily, and possibly upside down. Nothing in the data says which source is “first” or what its units are. You identify them afterwards by what they look like, which is what the EEG example is about.

A real cocktail party: 64 channels of EEG

The file holds one minute of resting EEG from one person, 64 electrodes at 160 samples per second, from the PhysioNet motor-imagery dataset. Plot four channels. Fp1 sits on the forehead above the left eye.

eeg = np.load("eeg_s001r01.npz")
E = eeg["eeg"].astype(float)                               # 64 channels x 9760 samples, microvolts
fs = float(eeg["fs"])                                      # 160 samples per second
names = list(eeg["channels"])
E = E - E.mean(axis=1, keepdims=True)
te = np.arange(E.shape[1]) / fs
print(E.shape, "channels x samples,", E.shape[1] / fs, "seconds")
# (64, 9760) channels x samples, 61.0 seconds

fig, axes = plt.subplots(4, 1, figsize=(10, 6), sharex=True)
for ax, ch in zip(axes, ["Fp1", "Fz", "Cz", "Oz"]):
    ax.plot(te, E[names.index(ch)], lw=0.7)
    ax.set(ylabel=ch, xlim=(0, 20))
axes[-1].set_xlabel("time (s)")
plt.show()
Twenty seconds on four channels. The large slow deflections on Fp1, several hundred microvolts, are eye blinks: the eyeball is a dipole, and turning it under the lid produces a voltage far bigger than anything the cortex makes. The same blinks are visible, smaller, on Fz and even Cz.
Twenty seconds on four channels. The large slow deflections on Fp1, several hundred microvolts, are eye blinks: the eyeball is a dipole, and turning it under the lid produces a voltage far bigger than anything the cortex makes. The same blinks are visible, smaller, on Fz and even Cz.

Blinks are the classic ICA target because they satisfy every assumption: one source, fixed location, mixed into every channel in fixed proportions, and about as non-Gaussian as a signal gets. Run the same two functions on 64 channels, keeping 20 components, and ask which component looks most like the forehead.

Z, unwhiten = whiten(E, 20)                                # keep the 20 strongest components of 64 channels
W = fastica(Z)
sources = W @ Z                                            # 20 independent components x samples
patterns = unwhiten @ W.T                                  # 64 x 20: how much of each component each channel sees

front = (E[names.index("Fp1")] + E[names.index("Fp2")]) / 2
similarity = [abs(np.corrcoef(s, front)[0, 1]) for s in sources]
blink = int(np.argmax(similarity))
print(f"component {blink} looks most like the forehead channels: correlation {similarity[blink]:.2f}")
strongest = np.argsort(-np.abs(patterns[:, blink]))[:6]
print("channels where it is strongest:", [names[i] for i in strongest])
# component 18 looks most like the forehead channels: correlation 0.93
# channels where it is strongest: ['Fp1', 'Fp2', 'Af8', 'Fpz', 'Af7', 'Af3']
Left: six of the twenty components, with the blink component on top: three large blinks in twenty seconds, and nothing else. Right: how much of that component each electrode sees, drawn on a schematic head. It lives on the forehead, strongest at Fp1, Fpz and Fp2, and fades to nothing by the central row.
Left: six of the twenty components, with the blink component on top: three large blinks in twenty seconds, and nothing else. Right: how much of that component each electrode sees, drawn on a schematic head. It lives on the forehead, strongest at Fp1, Fpz and Fp2, and fades to nothing by the central row.

The component the algorithm found correlates 0.93 with the forehead channels, and its scalp pattern is a textbook blink: largest at the three electrodes nearest the eyes, falling off smoothly backwards. Nothing told the algorithm where the eyes are. The pattern comes from the mixing proportions it discovered, which are the physics of a dipole under the forehead.

Removing it

Because ICA gives you both the sources and the mixing pattern, you can set one source to zero and rebuild the channels without it. That is artifact removal, and it is what every EEG pipeline does with blinks, heartbeats and muscle.

sources_clean = sources.copy()
sources_clean[blink] = 0                                   # silence the blink component
E_clean = patterns @ sources_clean                         # and rebuild the 64 channels without it

for ch in ["Fp1", "Fz", "Cz"]:
    i = names.index(ch)
    print(f"{ch}: standard deviation {E[i].std():5.1f} uV before, {E_clean[i].std():5.1f} uV after")
# Fp1: standard deviation 110.3 uV before,  46.1 uV after
# Fz: standard deviation  60.8 uV before,  41.6 uV after
# Cz: standard deviation  54.1 uV before,  47.7 uV after

fig, axes = plt.subplots(3, 1, figsize=(10, 5.5), sharex=True)
axes[0].plot(te, sources[blink], lw=0.8)
axes[0].set(ylabel="blink component")
axes[1].plot(te, E[names.index("Fp1")], lw=0.7, label="Fp1 before")
axes[1].plot(te, E_clean[names.index("Fp1")], lw=0.7, label="Fp1 after")
axes[1].legend()
axes[2].plot(te, E[names.index("Cz")], lw=0.7, label="Cz before")
axes[2].plot(te, E_clean[names.index("Cz")], lw=0.7, label="Cz after")
axes[2].legend()
axes[2].set(xlabel="time (s)", xlim=(0, 20))
plt.show()
Top: the blink component. Middle: Fp1 before and after removing it; the blinks are gone and the brain signal underneath is untouched. Bottom: the same for Cz, where the blinks were small to begin with and the change is slight.
Top: the blink component. Middle: Fp1 before and after removing it; the blinks are gone and the brain signal underneath is untouched. Bottom: the same for Cz, where the blinks were small to begin with and the change is slight.

Fp1’s standard deviation drops from 110 to 46 microvolts, which is the blinks leaving. Cz barely changes, because little of the blink reached it. What remains on Fp1 is the EEG that was always there under the blinks, and that is what you analyse.

What ICA assumes, and when it fails

  • Sources are independent. Not just uncorrelated. Two brain regions that co-activate are not independent and ICA may split them strangely.
  • At most one source is Gaussian. The method finds non-Gaussianity; two Gaussian sources cannot be told apart by any rotation.
  • Mixing is linear and fixed. True for the volume conduction of EEG and MEG, which is why ICA works so well there. Not true if the source moves, which is why a subject who shifts in the chair halfway through gives two blink components.
  • You need enough data. The usual rule of thumb is at least 20 to 30 times the square of the number of components. One minute at 160 Hz, 9,760 samples, is enough for 20 components (20 x 400 = 8,000) and nowhere near enough for 64 (which would want 80,000 or more).
  • It does not tell you what the components are. Order, sign and scale are arbitrary, and identifying the blink was our job. Automated labelling exists (ICLabel is the common one) but it is a classifier trained on human labels, not part of ICA.

Exercises

  1. Run fastica with a different seed. Which components change order? Which change sign? Do any change shape?
  2. Keep 10 components instead of 20, then 40. At what point does the blink component stop being clean?
  3. Find the component with the highest kurtosis (((s - s.mean())**4).mean() / s.var()**2 - 3) instead of the highest correlation with Fp1. Is it the blink? If not, plot it and its scalp pattern and work out what it is.
  4. Harder: the recording also contains a heartbeat somewhere, a sharp spike about once a second with a pattern that is strongest at the back and sides of the head. Find it.

Next primer: Bayesian inference. Data: EEG Motor Movement/Imagery Dataset (Schalk et al. 2004) via PhysioNet (Goldberger et al. 2000), subject 1, run 1, Open Data Commons Attribution licence. FastICA: Hyvärinen and Oja, 2000.

Keep reading