Symmetric matrices and PCA

Here are five points in the plane: $(3,\allowbreak 3)$, $(-3,\allowbreak -3)$, $(1,\allowbreak -1)$, $(-1,\allowbreak 1)$ and $(0,\allowbreak 0)$. Plot them and one question answers itself: along which direction do they spread out most? Along the diagonal line through $(1,\allowbreak 1)$. Two of the points sit on that line, three times as far from the middle as the two on the other diagonal.

Now ask the same question of 10,000 sentence embeddings, each a list of 768 numbers. No picture helps. What you want is a calculation that takes the data and returns the direction of greatest spread, then the next one, and so on. This lesson builds it, and the answer is an eigenvector. Not of any matrix, though: of a symmetric one, and symmetric matrices have unusually clean eigenvectors.

We start with what symmetric means and what it does to eigenvectors, prove it, and then use it on data. The same facts carry the next lesson, on singular values.

#What symmetric means

Most matrices from the last two lessons had eigenvectors at odd angles to each other: $(1,\allowbreak 1)$ and $(1,\allowbreak -2)$ for $\begin{bmatrix} 4 & 1 \\ 2 & 3 \end{bmatrix}$. One family of matrices never does that, and it is the family machine learning meets most.

A matrix is symmetric if it equals its own transpose, $\mathbf{A}^\top = \mathbf{A}$: the entry in row $i$, column $j$ equals the entry in row $j$, column $i$, so the matrix is a mirror image across its diagonal. $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ is symmetric; $\begin{bmatrix} 4 & 1 \\ 2 & 3 \end{bmatrix}$ is not, because its top-right 1 differs from its bottom-left 2. Symmetric matrices are everywhere in machine learning: covariance matrices (this lesson), the Hessian of a loss (the matrix of second derivatives, in gradient descent), and $\mathbf{A}^\top\mathbf{A}$ for any matrix $\mathbf{A}$ (the next lesson) are all symmetric.

See what symmetry does before we prove anything. The instrument below is the matrix instrument again, starting at the symmetric matrix $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$, with the lime eigen-lines switched on and the eigenvalues in the readouts. The coral arrow is the second column, $(1,\allowbreak 2)$. Dragging its tip sideways changes the top-right entry and leaves the bottom-left one alone, which breaks the symmetry.

In every symmetric case the eigen-lines were perpendicular, and there were always two of them. That is a theorem.

#The spectral theorem

The picture suggests a rule. Here it is in full, followed by the two proofs we need for 2x2 matrices.

This is the diagonalization $\mathbf{A} = \mathbf{P}\mathbf{D}\mathbf{P}^{-1}$ of the previous lesson with the best possible $\mathbf{P}$, one whose inverse is just its transpose. A matrix like $\mathbf{Q}$, whose columns have length 1 and are perpendicular to each other, is called orthogonal. Three facts about it:

• Its inverse is its transpose. Entry $(i,\allowbreak j)$ of $\mathbf{Q}^\top\mathbf{Q}$ is row $i$ of $\mathbf{Q}^\top$ dotted with column $j$ of $\mathbf{Q}$, which is column $i$ dotted with column $j$. That is 1 when $i = j$ (length 1) and 0 otherwise (perpendicular), so $\mathbf{Q}^\top\mathbf{Q} = \mathbf{I}$.
• It keeps every length. Write the squared length as a dot product and use the first fact: $$|\mathbf{Q}\mathbf{x}|^2 = (\mathbf{Q}\mathbf{x})^\top(\mathbf{Q}\mathbf{x}) = \mathbf{x}^\top\mathbf{Q}^\top\mathbf{Q}\,\mathbf{x} = \mathbf{x}^\top\mathbf{x} = |\mathbf{x}|^2.$$
• It keeps every angle. The same steps show $(\mathbf{Q}\mathbf{x})\cdot(\mathbf{Q}\mathbf{y}) = \mathbf{x}^\top\mathbf{Q}^\top\mathbf{Q}\,\mathbf{y} = \mathbf{x}^\top\mathbf{y} = \mathbf{x}\cdot\mathbf{y}$, and the angle between two vectors is fixed by their dot product and their lengths. So in the plane an orthogonal matrix is a rotation or a reflection (in more dimensions, a combination of them).

Read $\mathbf{A} = \mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top$ right to left, as before. $\mathbf{Q}^\top$ rotates (or reflects) the eigen-directions onto the axes, $\boldsymbol{\Lambda}$ stretches each axis by its eigenvalue (a negative one also flips it), and $\mathbf{Q}$ rotates back. No shear anywhere: a symmetric matrix is a pure stretch along perpendicular axes.

We will prove the two halves of the theorem for the cases we need. First, why the eigenvalues are real in the 2x2 case.

Second, why the eigenvectors are perpendicular. The proof is three lines. It uses only the transpose rule $(\mathbf{A}\mathbf{B})^\top = \mathbf{B}^\top\mathbf{A}^\top$ from shapes, transposes and batches and the dot product written as a matrix product, $\mathbf{a}\cdot\mathbf{b} = \mathbf{a}^\top\mathbf{b}$.

Worked example. $\mathbf{A} = \begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ has trace 4 and determinant $4 - 1 = 3$, so $\lambda^2 - 4\lambda + 3 = (\lambda - 3)(\lambda - 1) = 0$ and the eigenvalues are 3 and 1. For $\lambda = 3$, the row $(-1,\allowbreak 1)$ of $\mathbf{A} - 3\mathbf{I}$ gives the eigenvector $(1,\allowbreak 1)$. For $\lambda = 1$, the row $(1,\allowbreak 1)$ of $\mathbf{A} - \mathbf{I}$ gives $(1,\allowbreak -1)$. Their dot product is $1 \cdot 1 + 1 \cdot (-1) = 0$, as promised. Both have length $\sqrt{2}$, so dividing by $\sqrt{2}$ gives unit vectors:

Multiply back to check. $\mathbf{Q}\boldsymbol{\Lambda}$ scales the columns of $\mathbf{Q}$ by 3 and 1, and this $\mathbf{Q}$ happens to equal its own transpose, so

Compare the non-symmetric $\begin{bmatrix} 4 & 1 \\ 2 & 3 \end{bmatrix}$ from the first lesson: its eigenvectors $(1,\allowbreak 1)$ and $(1,\allowbreak -2)$ have dot product $1 - 2 = -1$, so they are not perpendicular.

To recap: symmetric means $\mathbf{A}^\top = \mathbf{A}$; such a matrix has real eigenvalues and perpendicular eigenvectors, so it factors as $\mathbf{Q}\boldsymbol{\Lambda}\mathbf{Q}^\top$ with an orthogonal $\mathbf{Q}$ that keeps lengths and angles. Take one apart yourself, including the check.

#Principal components

Now the payoff promised at the start: the direction in which data varies most. It takes three steps. Measure how the features vary together (the covariance matrix), turn "the spread along a direction" into a formula, and maximize that formula.

The data and its covariance. You have $n$ examples, each with $d$ features, stored as the rows of a data matrix $\mathbf{X}$ of shape $n \times d$ (the row convention used in code). First center the data: subtract each column's mean, so every feature averages to zero. Call the result $\tilde{\mathbf{X}}$. The covariance matrix is

It is $d \times d$. Its diagonal entry $C_{jj}$ is the variance of feature $j$, the average squared distance from the mean. The off-diagonal entry $C_{jk}$ is the covariance of features $j$ and $k$, positive when they tend to rise together. It is symmetric because $C_{jk}$ and $C_{kj}$ are the same sum. (Dividing by $n - 1$ is the usual sample convention, and it is what np.cov does. Module 0 divided by $n$, as np.var does; the two choices differ by the same factor for every variance and change no direction.)

For the five points from the opening, the mean is $(0,\allowbreak 0)$, so they are already centered. Add up the three kinds of product over the points:

Divide by $n - 1 = 4$:

Each feature has variance 5, and the covariance 4 is positive: when $x$ is large, $y$ tends to be large too.

The spread along a direction. Take a unit vector $\mathbf{u}$ and project every centered example onto it: $z_i = \mathbf{u}\cdot\tilde{\mathbf{x}}_i$, where $\tilde{\mathbf{x}}_i$ is row $i$ of $\tilde{\mathbf{X}}$ written as a column. Each $z_i$ is how far example $i$ lies along $\mathbf{u}$. The $z_i$ average to zero, and their variance turns out to be a single number built from $\mathbf{u}$, $\mathbf{C}$ and $\mathbf{u}$ again (an expression of this kind is called a quadratic form):

Check it on the five points with $\mathbf{u} = (1,\allowbreak 1)/\sqrt{2}$. The projections are $(3 + 3)/\sqrt{2} = 3\sqrt{2}$, then $-3\sqrt{2}$, then $(1 - 1)/\sqrt{2} = 0$, $0$ and $0$. Their squares add to $18 + 18 = 36$, and $36/4 = 9$. The formula gives the same: $\mathbf{C}\mathbf{u} = (5 + 4,\allowbreak \ 4 + 5)/\sqrt{2} = (9,\allowbreak 9)/\sqrt{2}$, and $\mathbf{u}^\top\mathbf{C}\mathbf{u} = (9 + 9)/2 = 9$. Along the x-axis, $\mathbf{u} = (1,\allowbreak 0)$, it gives $C_{11} = 5$, and along $(1,\allowbreak -1)/\sqrt{2}$ it gives $(5 - 4 - 4 + 5)/2 = 1$. The diagonal really is the long direction.

The best direction is the top eigenvector. $\mathbf{C}$ is symmetric, so by the spectral theorem it has perpendicular unit eigenvectors $\mathbf{q}_1,\allowbreak \ldots,\allowbreak \mathbf{q}_d$. Number them so the eigenvalues come largest first, $\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_d \ge 0$. Any unit $\mathbf{u}$ can be written as $\mathbf{u} = c_1\mathbf{q}_1 + \cdots + c_d\mathbf{q}_d$, and then $c_1^2 + \cdots + c_d^2 = 1$. The reason: $|\mathbf{u}|^2 = \mathbf{u}\cdot\mathbf{u}$, and multiplying out that dot product gives a term $c_i c_j\,(\mathbf{q}_i\cdot\mathbf{q}_j)$ for every pair; the pairs with $i \ne j$ are 0 (perpendicular) and each $\mathbf{q}_i\cdot\mathbf{q}_i$ is 1, so $|\mathbf{u}|^2 = c_1^2 + \cdots + c_d^2$, which is 1 for a unit vector. The box below shows that then

with equality when $\mathbf{u} = \mathbf{q}_1$. The variance along $\mathbf{u}$ is an average of the eigenvalues, weighted by $c_i^2$, and an average is never above its largest value.

Check it on the five points. Their covariance has the unit eigenvectors $\mathbf{q}_1 = (1,\allowbreak 1)/\sqrt{2}$ and $\mathbf{q}_2 = (1,\allowbreak -1)/\sqrt{2}$, with eigenvalues 9 and 1 (the worked example below confirms this; they are the variances along the two diagonals found above). The x-axis $\mathbf{u} = (1,\allowbreak 0)$ has $c_1 = \mathbf{u}\cdot\mathbf{q}_1 = 1/\sqrt{2}$ and $c_2 = \mathbf{u}\cdot\mathbf{q}_2 = 1/\sqrt{2}$ (dotting $\mathbf{u} = c_1\mathbf{q}_1 + c_2\mathbf{q}_2$ with $\mathbf{q}_1$ leaves only $c_1$, because $\mathbf{q}_1\cdot\mathbf{q}_1 = 1$ and $\mathbf{q}_1\cdot\mathbf{q}_2 = 0$), so $c_1^2 = c_2^2 = \tfrac12$, and the formula gives $9 \cdot \tfrac12 + 1 \cdot \tfrac12 = 5$: the variance along the x-axis, $C_{11}$, halfway between the two eigenvalues.

So the direction of greatest variance is the top eigenvector, and the variance along it is the top eigenvalue. The same argument, restricted to directions perpendicular to $\mathbf{q}_1$, picks out $\mathbf{q}_2$, and so on down the list.

Worked example. For the five points, $\mathbf{C} = \begin{bmatrix} 5 & 4 \\ 4 & 5 \end{bmatrix}$ has the same pattern as $\begin{bmatrix} 2 & 1 \\ 1 & 2 \end{bmatrix}$ (equal diagonal entries, equal off-diagonal entries), so its eigenvectors are again $(1,\allowbreak 1)$ and $(1,\allowbreak -1)$. The eigenvalues are $5 + 4 = 9$ and $5 - 4 = 1$: check that $\mathbf{C}(1,\allowbreak 1) = (9,\allowbreak 9)$ and $\mathbf{C}(1,\allowbreak -1) = (1,\allowbreak -1)$, and that $9 + 1 = 10$ is the trace and $9 \times 1 = 25 - 16$ the determinant. The first principal component is $(1,\allowbreak 1)/\sqrt{2}$, with variance 9, the number we computed from the projections above. It carries $9/(9 + 1) = 90\%$ of the total variance.

Every other direction has less spread than the first component. Measure one yourself.

Now a larger cloud: 200 random points stretched by 2 along one axis and 0.5 along the other, then tilted by 30° and moved away from the origin, so you can watch PCA recover the tilt.

The code prints an angle of 27.5°, close to the 30° tilt but not equal to it, because 200 random points only approximate the shape they were drawn from. For the same reason the two variances come out as 3.352 and 0.199 rather than $2^2 = 4$ and $0.5^2 = 0.25$, the spreads the cloud was drawn with. The first component still holds 94.4% of the variance. Use np.linalg.eigh, not eig, for symmetric matrices: it is faster, guarantees real results and perpendicular eigenvectors, and sorts the eigenvalues (in increasing order, which is why the code reverses them).

To recap: center the data, form the symmetric covariance matrix $\mathbf{C}$, and its eigenvectors, largest eigenvalue first, are the principal components; each eigenvalue is the variance along its direction. Now implement it, and test it on the five points and on a 3D cloud with a known shape.