Finding better axes

Data arrives in the coordinates of whatever measured it: one number per pixel, per microphone, per survey question. Those axes suit the instrument. They rarely suit the question, because the things you care about (a face, a voice, a preference) are smeared across many of them at once.

Here is the problem in two dimensions. Two hidden signals were blended and picked up by two microphones. The same cloud of points is shown with three different choices of axes.

The page has three parts, in this order:

• Part 1, the singular value decomposition (SVD). What any matrix does to space: a rotation, a stretch, and another rotation.
• Part 2, principal component analysis (PCA). The same decomposition pointed at a table of data. It finds the perpendicular axes along which the data vary most.
• Part 3, independent component analysis (ICA). The step PCA cannot take. It finds axes along which the data are statistically independent, and those need not be perpendicular.

SVD comes first even though it is the most abstract of the three, because PCA is easiest to understand as an SVD, and ICA is easiest to understand as what is left to do after PCA.

You will need vectors, matrix multiplication and transposes, and the ideas of mean and variance. Each part ends with a quiz, with single-question checks and two hands-on challenges along the way. The sidebar keeps score in this browser. The page draws on nine explainers, credited where their ideas appear and listed at the end.

1.1A matrix moves the plane

After Gregory Gundersen, Singular Value Decomposition as Simply as Possible.

A \(2\times 2\) matrix \(M\) takes each point \(\mathbf x\) of the plane to a new point \(M\mathbf x\). The move is linear: straight lines stay straight, parallel lines stay parallel, and the origin stays where it is. Within those rules a matrix can stretch, squash, rotate, flip and shear, in any combination. To know the whole matrix it is enough to know where it sends the basis vectors \(\mathbf e_1\) and \(\mathbf e_2\). Those two landing points are its columns.

In the figure the left panel is the input and the right panel is the output. Drag the two square handles on the right to change the matrix. Drag the probe \(\mathbf x\) around the circle on the left and watch where it lands.

Whatever matrix you make, three things hold.

1. The circle becomes an ellipse. Sometimes a degenerate one: the Flatten preset squashes it onto a line.
2. The ellipse has a longest and a shortest radius. Call their lengths \(\sigma_1 \ge \sigma_2\) and their directions \(\mathbf u_1\) and \(\mathbf u_2\). The two are perpendicular, as the axes of an ellipse always are.
3. The inputs that land on those radii, \(\mathbf v_1\) and \(\mathbf v_2\), are perpendicular too. A shear tilts almost every right angle, yet there is always one pair of perpendicular inputs whose right angle survives the trip.

In symbols,

The \(\sigma_i\) are the singular values of \(M\), the \(\mathbf v_i\) its right singular vectors and the \(\mathbf u_i\) its left singular vectors. Gundersen states the whole idea at the level of a fourteen-year-old: turn a square the right way first, and the matrix only stretches it into a rectangle, which may come out turned.

1.2Rotate, stretch, rotate

Those two equations already contain the factorization. Put \(\mathbf v_1, \mathbf v_2\) side by side as the columns of a matrix \(V\), put \(\mathbf u_1, \mathbf u_2\) into \(U\), and put \(\sigma_1, \sigma_2\) on the diagonal of \(\Sigma\). Both equations together then read \(MV = U\Sigma\). The columns of \(V\) are perpendicular unit vectors, which makes \(V\) an orthogonal matrix: multiplying it by its transpose gives the identity, \(VV^{\mathsf T} = I\). Multiply both sides on the right by \(V^{\mathsf T}\):

Read from right to left, this is a recipe.

1. \(V^{\mathsf T}\) rotates the plane so that \(\mathbf v_1\) and \(\mathbf v_2\) lie along the coordinate axes.
2. \(\Sigma\) stretches along those axes, by \(\sigma_1\) and \(\sigma_2\).
3. \(U\) rotates the axes onto \(\mathbf u_1\) and \(\mathbf u_2\).

Orthogonal matrices preserve lengths and angles, so they are rotations, possibly combined with a mirror flip. All of the change of shape is in \(\Sigma\).

Nothing here is special to two dimensions or to square matrices. An \(m\times n\) matrix maps the unit sphere of \(\mathbb R^n\) to an ellipsoid in \(\mathbb R^m\). \(V\) is \(n\times n\), \(U\) is \(m\times m\), and \(\Sigma\) is \(m\times n\) with the singular values down its diagonal, by convention non-negative and sorted from largest to smallest. For picturing it, Gundersen passes on Geoffrey Hinton's advice: visualize three dimensions and say "fourteen" to yourself very loudly.

1.3One direction at a time

After Callum McDougall, Six (and a half) intuitions for SVD, intuitions 4 and 5.

Multiply the three factors out and the matrix becomes a sum of simple pieces, one per singular value:

Each term works like a channel with three parts. \(\mathbf v_i\) is an input direction: the channel measures how much of \(\mathbf x\) points that way. \(\sigma_i\) is a gain. \(\mathbf u_i\) is an output direction: the channel writes its result along it. The channels do not interfere with one another, because the input directions are mutually perpendicular and so are the output directions.

Read this way, the singular values say what a matrix keeps and what it destroys.

• If \(\sigma_i = 0\) the channel is dead. Whatever \(\mathbf x\) had along \(\mathbf v_i\) is erased. The number of nonzero singular values is the rank of the matrix.
• If \(\sigma_i\) is tiny, that part of the input survives only faintly. Getting it back means dividing by a tiny number, which amplifies any noise.
• If every \(\sigma_i\) is positive (and the matrix is square), the matrix can be undone one channel at a time: \(M^{-1} = V\,\Sigma^{-1}\,U^{\mathsf T}\). Rotate back, divide by each gain, rotate back.

1.4A matrix can also be a table

After Jeremy Kun, SVD Part 1: Perspectives on Linear Algebra and Part 2, and Peter Bloem, The Singular Value Decomposition.

So far \(M\) has been a transformation. More often a matrix is a table: rows are movies and columns are viewers, or rows are days and columns are sensors, or the entries are the pixels of a picture. Kun points out that these are one object seen twice (the table's columns are where the map sends the basis vectors), and that one SVD describes both.

For a table, the sum of pieces is the useful form. Each piece \(\sigma_i\mathbf u_i\mathbf v_i^{\mathsf T}\) is a table of the simplest possible kind, in which every row is a multiple of one pattern \(\mathbf v_i^{\mathsf T}\). The pieces come sorted by size, so you can stop early:

The Eckart–Young theorem says that this shortcut is the best one available. Among all matrices of rank \(k\), \(A_k\) is the closest to \(A\) in total squared error, and that error is exactly \(\sigma_{k+1}^2 + \sigma_{k+2}^2 + \cdots\) Below, a 128 × 128 grayscale image is treated as a matrix of brightness values.

Try the four images in turn.

• Landscape is smooth and its rows resemble one another, so the singular values fall away quickly and a few dozen layers are plenty.
• Lettering has sharp diagonal and curved edges, which no small set of row patterns can describe. It needs more layers before the letters are clean.
• Plaid was made by adding three patterns together. Three layers reproduce it exactly, and every later singular value is zero.
• Static has no shared structure for the layers to exploit. The first layer captures its average brightness. After that the singular values barely fall, and the approximation stays poor until nearly all of them are in.

Real data fall between Plaid and Static, as Kun points out. Measurements are almost never exactly low rank, because noise spills into every direction. But the structure is often concentrated in a few layers while the noise is spread thinly over all of them, so cutting the sum short keeps most of the structure and discards mostly noise.

2.1Two ways to summarize a cloud

After amoeba's answer on Cross Validated, Making sense of principal component analysis, eigenvectors & eigenvalues (CC BY-SA 4.0).

Suppose you have measured a cellar of wines: how dark each one is, how strong, how old, how acidic. Many of these properties go together, so the table says less than its number of columns suggests. PCA summarizes it with a few new features. It does not pick some columns and throw the rest away. Each new feature is a weighted combination of all the old ones (darkness plus strength minus acidity, say), chosen to summarize the wines as well as possible.

amoeba gives two meanings of "summarize well", and they sound different.

• Spread. A feature on which every wine scores about the same tells you nothing. A good feature has high variance.
• Reconstruction. A good feature lets you guess the original measurements back from it with little error.

Take just two measurements, so each wine is a point in the plane, and move the origin to the center of the cloud. A candidate feature is a line through the origin, and a wine's score on that feature is where its perpendicular projection lands on the line. Turn the line and judge it both ways.

The two criteria never disagree, and the picture shows why. Each blue point, its orange projection and the origin form a right triangle. For a unit vector \(\mathbf w\) along the line, Pythagoras gives

Average over the wines: total variance = variance along the line + variance across it. The first is what the feature captures and the second is the reconstruction error. The total does not depend on the line, so whatever makes one term large makes the other small. The line that wins is the first principal component.

The springs button runs amoeba's physical version. Picture the line as a rod and each thin orange segment as a spring. A spring stores energy in proportion to the square of its length, so the rod comes to rest where the sum of squared distances is smallest.

2.2Eigenvectors of the covariance matrix

The small matrix in Figure 2.1 is the covariance matrix of the cloud, measured in the frame made by the line and its perpendicular. In the original axes it is

where \(X\) holds the centered data, one wine per row. The diagonal entries are the spread along each axis. The off-diagonal entry says how strongly the two coordinates move together. The variance along any unit direction \(\mathbf w\) is \(\mathbf w^{\mathsf T}C\,\mathbf w\).

Watch the matrix as the line nears its best angle. Two things happen together: the top-left entry reaches its maximum, and the off-diagonal entries pass through zero. In that frame the matrix is diagonal, which means the two new coordinates are uncorrelated.

An eigenvector of a matrix is a direction that the matrix stretches without turning: \(C\mathbf v = \lambda\mathbf v\), where the stretch factor \(\lambda\) is its eigenvalue. A frame in which a symmetric matrix becomes diagonal is a frame made of its eigenvectors. The spectral theorem guarantees that every symmetric matrix has one, with perpendicular axes:

So the principal components are the eigenvectors of the covariance matrix, and each eigenvalue \(\lambda_i\) is the variance of the data along its eigenvector. With more than two features the story repeats: the second component is the best direction perpendicular to the first, the third is the best perpendicular to both, and so on. (Bloem proves the spectral theorem from scratch in part 3 of his series, for readers who want it.)

2.3PCA is the SVD of the centered data

After Gundersen (final section), Jonathon Shlens, A Tutorial on Principal Component Analysis §VI, and McDougall, intuition 6A.

Take the SVD of the centered data matrix itself, \(X = U\Sigma V^{\mathsf T}\), and substitute it into the covariance:

\(U\) drops out because \(U^{\mathsf T}U = I\), and \(\Sigma^{\mathsf T}\Sigma\) is a diagonal matrix holding the squares \(\sigma_i^2\). Compare with \(C = V\Lambda V^{\mathsf T}\): the \(V\) of the SVD is the matrix of eigenvectors, and there is nothing further to compute.

In Part 1, \(\mathbf v_1\) was the unit input that a matrix stretches most. Here the matrix is the data, and \(\lVert X\mathbf w\rVert^2\) is the sum of squared projections of every sample onto \(\mathbf w\), which is \(n-1\) times the variance along \(\mathbf w\). The most-stretched direction is the direction of greatest variance. And Eckart–Young is amoeba's second criterion: the first \(k\) components give the smallest reconstruction error that any \(k\) directions can.

In practice, center the data first: without that the SVD describes how far the data sit from the origin, not how they spread around their mean. Software also computes PCA from the SVD of \(X\) directly, because forming \(X^{\mathsf T}X\) squares every singular value and throws away numerical precision.

2.4Throwing dimensions away

After Shlens, A Tutorial on Principal Component Analysis §II–IV.

The point of new axes is that you may not need all of them. Shlens opens with an experimenter who does not know what she is looking at. A ball on a spring moves back and forth along a line. Not knowing which line, she films it with three cameras placed wherever was convenient. Each camera reports the ball's position in its own picture as two numbers, so every frame gives six. The motion has one degree of freedom, but the data have six columns, all of them noisy and none aligned with the spring.

After PCA, one component carries nearly all the variance, and its score over time is the motion of the spring. The other five sit at a common low level: they are the noise. That bar chart is called a scree plot, and the usual way to decide how many components to keep is to look for the point where it flattens. Raise the wobble and a second bar grows out of the floor, because the ball now has two independent motions. Raise the noise and the floor rises until the first bar no longer stands out.

Shlens frames the goal as two wishes about the covariance matrix: large entries on the diagonal (signal) and zeros off it (no redundancy). The six camera channels are highly redundant. The principal components are uncorrelated with one another.

Components as building blocks

After Bloem, A friendly introduction to PCA.

When each sample is an image, each principal direction is an image too, and every sample is the average image plus a weighted sum of them:

Bloem shows this with photographs of real faces. The figure below uses 160 cartoon faces, each drawn with a random head width, head height, eye spacing, eye size and smile. Each one is a 32 × 32 image, so 1,024 numbers.

Six weights stand in for 1,024 pixels. You cannot beat the starting weights by hand, because they are the projections of the face onto the components, and those are the least-squares optimum.

Look also at what each slider does. The faces were made by turning five separate knobs. The first two components happen to line up with two of them, head width and head height, because those change the most pixels. The later components match no single knob: each one changes the mouth, the eyes and the outline together. PCA finds directions of high variance, and nothing obliges those to coincide with the separate causes that produced the data. Part 3 takes that problem up.

2.5What PCA assumes

After Shlens, A Tutorial on Principal Component Analysis §VII.

PCA has no settings to tune. That makes it dependable, and it also means its assumptions are fixed in advance. Shlens lists them.

• Linearity. The new features are weighted sums of the old ones.
• Large variance means important. Directions of big spread are taken to be signal and directions of small spread to be noise.
• Perpendicular axes. Each component is orthogonal to all the others.

A fourth is practical: variance has units, so the answer depends on how each feature happens to be scaled.

The crossed streaks have two preferred directions, but PCA returns neither, because all it is asked to do is remove correlation. Its components are uncorrelated: there is no straight-line trend between them. That is weaker than independent, which means that knowing one tells you nothing at all about the other. For some problems independence is what you need.

3.1The cocktail party problem

After Aapo Hyvärinen & Erkki Oja, Independent Component Analysis: Algorithms and Applications §1–2.

Two people talk at once in a room with two microphones. Each microphone records a weighted sum of the two voices, with weights that depend on where it stands:

You are given the recordings \(\mathbf x\) and nothing else. You want the voices \(\mathbf s\). If you knew the mixing matrix \(A\) you would simply invert it. You do not, so you are looking for an unmixing matrix \(W\) with \(W\mathbf x = \mathbf s\), knowing neither \(A\) nor \(\mathbf s\).

Put that way it looks hopeless: one equation, with both things on the right unknown. ICA gets out of this by assuming that the sources are statistically independent: what one voice is doing at a given instant tells you nothing about what the other is doing. The task then becomes: find the \(W\) whose outputs are independent of each other. The rest of this part shows that this is very nearly enough to pin \(W\) down, and why PCA's weaker goal of uncorrelated outputs is not.

The model is not limited to voices. In Arnaud Delorme's introduction for EEG researchers, the sensors are electrodes on the scalp and the sources are patches of cortex, muscles and blinking eyes.

3.2Whitening: how far PCA takes you

After Shlens, A Tutorial on Independent Component Analysis §V–VII, and Delorme.

Shlens organizes the whole problem around the SVD of the unknown mixing matrix, \(A = U\Sigma V^{\mathsf T}\). Unmixing means undoing its three factors in reverse order,

which is the channel-by-channel inverse from Part 1. (These are not the \(U\), \(\Sigma\) and \(V\) of Part 2. There the SVD was of the data table \(X\). Here it is of the mixing matrix \(A\).) Two of the three factors can be read off the covariance of the recordings. Assume the sources have variance 1 (their true loudness is something we could never determine anyway). Independent sources are uncorrelated, so their covariance matrix is the identity, and

\(V\) has cancelled. Diagonalizing the covariance of the recordings, which is what PCA does, gives us \(U\) and \(\Sigma\). Apply \(\Sigma^{-1}U^{\mathsf T}\) to the data:

This step is called whitening or sphering. The whitened data \(\mathbf z\) have identity covariance, and they are the sources themselves, turned by an unknown rotation \(V^{\mathsf T}\).

3.3The rotation covariance cannot see

Why not carry on with the covariance and find \(V\) as well? Because whitened data look the same to the covariance from every angle. Turn \(\mathbf z\) by any orthogonal matrix \(R\):

Every rotation of the whitened cloud is uncorrelated with unit variances. PCA cannot tell them apart and stops here.

To go further we need the distinction from the end of Part 2. Two quantities are uncorrelated when there is no linear trend between them. They are independent when knowing one tells you nothing at all about the other, which is to say their joint distribution is the product of their separate distributions. Independence implies zero correlation. The reverse fails. In the last panel of Figure 3.1 the square sits tilted, and its two coordinates are uncorrelated, yet if the first is near its largest possible value you are at a corner, and the second is confined to a narrow range. Knowing one tells you about the other.

So the plan is to find, among all rotations of the whitened data, the one whose coordinates are independent. For that we need a way to tell when we have arrived.

3.4Blends look more Gaussian

After Hyvärinen & Oja §4, and Delorme.

Measuring independence head-on is awkward. Hyvärinen and Oja give a shortcut that rests on the central limit theorem: add independent random variables together, and the sum tends to be closer to a Gaussian than the variables were.

Now take any single coordinate of the whitened data. It is some weighted sum of the sources, \(y = q_1s_1 + q_2s_2\). If both weights are nonzero it is a blend, and a blend tends to be more Gaussian than its ingredients. It is least Gaussian when one of the weights is zero, that is, when \(y\) is a single source. That gives ICA its working principle:

The simplest gauge of non-Gaussianity is excess kurtosis. For a variable with mean 0 and variance 1 it is \(\operatorname{kurt}(y) = \mathbb E[y^4] - 3\). It is zero for a Gaussian, negative for flat-topped shapes such as the uniform distribution (−1.2), and positive for spiky, heavy-tailed ones such as speech. Either sign counts as non-Gaussian, so the quantity to maximize is its absolute value.

3.5Finding the last rotation

The figure starts from whitened recordings of two independent sources. Turn the axes. The two histograms show the data's coordinates along the two axes.

Three things to try.

• Flat sources. The cloud is a tilted square. When the axes line up with its sides, both histograms turn into flat-topped boxes and the kurtoses reach −1.2.
• Spiky sources. The cloud has arms along the two source directions. At the right angle the histograms are sharply peaked and the kurtoses are large and positive.
• Gaussian sources. The histograms and the kurtoses look the same at every angle.

That last case is why Gaussian sources are forbidden. The density of whitened Gaussian data depends only on distance from the origin, so the cloud is perfectly round. Every rotation of it gives coordinates that are independent and Gaussian. All of them are valid answers and the true sources cannot be told apart from the rest. ICA needs non-Gaussian sources and can tolerate at most one Gaussian among them, since that one is whatever remains when the others have been found.

The answer also repeats every 90°. Turning the axes a quarter turn swaps the two components and flips a sign, which is the same solution with different labels. Section 3.6 returns to this.

Notice also that the scatter plot has no time axis. Basic ICA sees only this cloud of points, so the order in which the samples arrived plays no part in the answer.

FastICA and Infomax

Nobody turns a dial. FastICA, from Hyvärinen and Oja, finds each direction with a fixed-point iteration that typically settles in a handful of steps. In place of kurtosis, which a few outliers can throw off badly, it uses a smoother approximation to negentropy: how far a distribution's entropy falls short of that of a Gaussian with the same variance. (The Gaussian has the largest entropy of any distribution with a given variance, so negentropy is zero for a Gaussian and positive for everything else.) Infomax, the algorithm Delorme describes and the usual choice for EEG, takes another route: it maximizes the entropy of the outputs after passing them through a squashing nonlinearity, which, when the nonlinearity suits the sources, amounts to minimizing the mutual information between them. Shlens minimizes the sum of the outputs' entropies directly. These objectives are closely related but not identical. On well-behaved data they give much the same answer, and on messy data they can disagree. On this page, Figure 3.3 searches the kurtosis curve it shows, while the opening figure and the sound demo below use FastICA with a tanh nonlinearity.

3.6Unmixing the party

After Hyvärinen & Oja §2–3, and Delorme.

Here is the whole pipeline on sound. Two sources, a melody and a drum pattern, play in a room. Each microphone hears both, more loudly the nearer it stands. Press ▶ on any row to listen, and drag the microphones to change the mix.

Listen to the microphones, then the whitened rows, then the ICA rows. The whitened outputs are uncorrelated and still audibly blends. ICA's are the sources, apart from three things it has no way to know.

• Order. Nothing marks one source as first. The melody may come out as IC 1 or IC 2.
• Scale. Double a source and halve its column of \(A\), and the recordings are unchanged. By convention each component is given variance 1.
• Sign. A special case of scale: multiply by −1. An upside-down waveform sounds the same.

None of these matters to the ear, but they do when you compare two runs or try to rank components.

Now put both microphones in the same spot. The two recordings become the same blend, differing only in each microphone's own hiss. The mixing matrix has no inverse, and all ICA can hand back is that blend and the hiss. Move one microphone a little and the separation returns, though while the two are close the result is fragile: this is the tiny singular value of section 1.3, where undoing the mix means dividing by a number near zero. ICA needs at least as many genuinely different recordings as there are sources. The full list of what the basic method assumes:

• The mixing is linear and instantaneous. Real rooms have echoes and delays, which break this. Separating real speech needs extensions of the method.
• The sources are independent.
• At most one source is Gaussian.
• There are at least as many sensors as sources, and each source reaches the sensors in its own pattern: the columns of the mixing matrix are linearly independent.

Delorme describes what this buys in EEG. Each row of \(W\) is a spatial filter that pulls one component out of the electrodes. Each column of \(W^{-1}\) is that component's scalp map, the fixed pattern with which it reaches every electrode. To remove an eye blink, find the component that carries it, set that row of the output to zero, and multiply back by \(W^{-1}\). The recording comes back with the blink largely gone. Anything else that ended up in that component goes with it, so the cleanup is only as good as the separation.

The SVD produces that factorization for any matrix: input directions, gains, output directions. PCA is the SVD of a centered data table, read statistically: \(V\) holds the axes of greatest variance and \(\Sigma\) the spread along them. For a mixing problem \(\mathbf x = A\mathbf s\), PCA recovers \(U\) and \(\Sigma\) from the covariance and stops, because the covariance contains nothing else. ICA recovers the remaining rotation \(V\) from the shape of the distribution.

Sources

This page is a synthesis of nine explainers, kept as local reading copies in the explainers folder next to this one. Their authors hold the rights to them. The words and figures here are new, but the ideas, the examples and often the order of the argument are theirs.