7.5 Applications to image processing and statistics

Big idea

Two applications, one factorization. An image stored as a grid of intensity values is a matrix. A data set of $N$ observations of $p$ measurements is also a matrix. In both cases the useful structure sits in a handful of directions, and the eigenvectors of a symmetric matrix say which directions those are.

For images the relevant fact is that the singular value expansion lists rank-one pieces in decreasing order of size. Keep the first few and discard the rest and you keep most of the picture while storing a fraction of the numbers. The size of what you threw away is measured by the singular values you dropped, so the error is known before you look at the result.

For data the relevant matrix is the covariance matrix, which records how the measured variables spread and how they move together. It is symmetric, so it diagonalizes orthogonally. Its eigenvectors point along the directions in which the data spreads most, its eigenvalues say how much spread each direction carries, and changing to those coordinates leaves the new variables uncorrelated. That is principal component analysis, and it is the spectral theorem applied to one specific matrix.

Decoder

The unit eigenvector of the covariance matrix with the largest eigenvalue is the first principal component of the data.

Unpacked: take the cloud of data points, move it so its center of mass sits at the origin, and look for the direction in which the cloud is longest. That direction is a unit vector, and it is the eigenvector belonging to the largest eigenvalue. The eigenvalue itself is the variance measured along that direction. “Component” here means a direction in measurement space, not a piece of the data.

Definitions and results

Data matrix. Put $N$ observations of $p$ measurements into a $p \times N$ matrix, one column per observation. The sample mean $\mathbf{M}$ is the average of the columns.

Mean-deviation form. Subtract $\mathbf{M}$ from every column to get $B$. The columns of $B$ sum to the zero vector, so the recentered cloud has its center at the origin. Every statement below is about $B$, never about the raw matrix.

Covariance matrix. The $p \times p$ matrix

$$ S = \frac{1}{N - 1}BB^{T} $$

Its $(j,j)$ entry is the sample variance of measurement $j$, and its $(i,j)$ entry is the covariance of measurements $i$ and $j$. Covariance is zero when the two measurements carry no linear relationship, positive when they tend to rise together, negative when one rises as the other falls. $S$ is symmetric and positive semidefinite, since it has the form $CC^{T}$.

Total variance. The trace of $S$, that is the sum of the individual variances. It is the single number that stands for “how much the data varies at all”.

Principal components. Orthogonally diagonalize $S$ as $PDP^{T}$ with eigenvalues ordered $\lambda_1 \geq \cdots \geq \lambda_p$. The columns of $P$ are the principal components of the data, first component first.

What the change of variable does. Set $\mathbf{y} = P^{T}\mathbf{x}$ for each recentered observation $\mathbf{x}$. The new variables $y_1, \dots, y_p$ have covariance matrix $D$: each pair is uncorrelated, and the variance of $y_i$ is $\lambda_i$. So the eigenvalues are literally the variances along the principal directions, and $y_1$ is the new variable with the largest variance of any unit-length combination of the originals.

Variance is conserved. The trace of $D$ equals the trace of $S$, because an orthogonal change of variable does not change the trace. So the total variance is redistributed, never created. The fraction of total variance explained by the $i$th component is $\lambda_i$ divided by the trace, and reporting the first few fractions is how a dimension reduction is justified.

Link to the previous section. Computing $S$ and diagonalizing it is one route; computing the singular value decomposition of $B$ directly is the other, and it finds the same directions with better numerical behavior. The left singular vectors of $B$ are the principal components, and each eigenvalue of $S$ is the square of the matching singular value of $B$ divided by $N - 1$.

Low-rank approximation. For a matrix $A$ of rank $r$ with singular value expansion $\sigma_1\mathbf{u}_1\mathbf{v}_1^{T} + \cdots + \sigma_r\mathbf{u}_r\mathbf{v}_r^{T}$, truncating after $k$ terms gives the best rank-$k$ approximation of $A$ in the least-squares sense. No other matrix of rank $k$ or less comes closer. The size of the error, measured by summing the squares of all entries of the difference and taking the square root, is $\sqrt{\sigma_{k+1}^2 + \cdots + \sigma_r^2}$.

Why that compresses an image. An $m \times n$ grid of intensities holds $mn$ numbers. A rank-$k$ truncation is described by $k$ vectors of length $m$, $k$ vectors of length $n$, and $k$ singular values, which is $k(m + n + 1)$ numbers. For a large image and a modest $k$ that is a large saving, and the singular values tell you in advance how much detail the truncation costs.

Worked examples

The best rank-one approximation of a small matrix

Take the matrix decomposed in the previous section,

$$ A = \begin{bmatrix} 3 & 0 \\ 4 & 5 \end{bmatrix} $$

with $\sigma_1 = 3\sqrt{5}$, $\sigma_2 = \sqrt{5}$, $\mathbf{u}_1 = \tfrac{1}{\sqrt{10}}(1,3)$ and $\mathbf{v}_1 = \tfrac{1}{\sqrt{2}}(1,1)$. The best rank-one approximation is the first term of the expansion:

$$ \sigma_1\mathbf{u}_1\mathbf{v}_1^{T} = \frac{3\sqrt{5}}{\sqrt{20}}\begin{bmatrix} 1 \\ 3 \end{bmatrix}\begin{bmatrix} 1 & 1 \end{bmatrix} = \frac{3}{2}\begin{bmatrix} 1 & 1 \\ 3 & 3 \end{bmatrix} = \begin{bmatrix} 1.5 & 1.5 \\ 4.5 & 4.5 \end{bmatrix} $$

using $\sqrt{20} = 2\sqrt{5}$. Both rows are multiples of $(1,1)$, so the rank is $1$ as promised. Subtract to see the error:

$$ A - \sigma_1\mathbf{u}_1\mathbf{v}_1^{T} = \begin{bmatrix} 1.5 & -1.5 \\ -0.5 & 0.5 \end{bmatrix} $$

Sum the squares of its entries: $2.25 + 2.25 + 0.25 + 0.25 = 5$, so the error size is $\sqrt{5}$. That is exactly $\sigma_2$, matching the rule that the discarded singular values measure the loss. With only two singular values to begin with there is nothing to compress here, but the arithmetic is the same arithmetic a large image goes through.

Principal components of a small data set

Five observations of two measurements:

$$ (2, 4), \quad (4, 2), \quad (5, 5), \quad (6, 8), \quad (8, 6) $$

Each measurement sums to $25$, so both means are $5$. Subtract to get the mean-deviation matrix

$$ B = \begin{bmatrix} -3 & -1 & 0 & 1 & 3 \\ -1 & -3 & 0 & 3 & 1 \end{bmatrix} $$

Confirm the centering: each row sums to $0$. Now form $BB^{T}$, whose entries are dot products of the rows:

• row one with itself: $9 + 1 + 0 + 1 + 9 = 20$
• row two with itself: $1 + 9 + 0 + 9 + 1 = 20$
• row one with row two: $3 + 3 + 0 + 3 + 3 = 12$

With $N - 1 = 4$,

$$ S = \frac{1}{4}\begin{bmatrix} 20 & 12 \\ 12 & 20 \end{bmatrix} = \begin{bmatrix} 5 & 3 \\ 3 & 5 \end{bmatrix} $$

Both measurements have variance $5$, they have covariance $3$, and the total variance is the trace, $10$.

The eigenvalues are $5 \pm 3$, that is $\lambda_1 = 8$ and $\lambda_2 = 2$, with unit eigenvectors $\tfrac{1}{\sqrt{2}}(1,1)$ and $\tfrac{1}{\sqrt{2}}(1,-1)$. So the first principal component points along $(1,1)$, which is the diagonal the cloud is stretched along, and it explains $8/10$, or $80$ percent, of the total variance. The second component picks up the remaining $20$ percent. Trace check: $8 + 2 = 10$.

Verify the variance claim rather than trusting it. The first new variable is $y_1 = \tfrac{1}{\sqrt{2}}(x_1 + x_2)$ applied to each recentered observation, giving

$$ \frac{1}{\sqrt{2}}(-4, \; -4, \; 0, \; 4, \; 4) $$

These already average to $0$, so the sample variance is the sum of squares divided by $4$: $\tfrac{1}{2}(16 + 16 + 0 + 16 + 16)/4 = 64/8 = 8$, which is $\lambda_1$. The second new variable is $y_2 = \tfrac{1}{\sqrt{2}}(x_1 - x_2)$, giving $\tfrac{1}{\sqrt{2}}(-2, 2, 0, -2, 2)$ and variance $\tfrac{1}{2}(4 + 4 + 0 + 4 + 4)/4 = 2$, which is $\lambda_2$.

Last, check that the new variables are uncorrelated. Multiply the two lists entry by entry and add: $\tfrac{1}{2}\left(8 - 8 + 0 - 8 + 8\right) = 0$. The covariance is zero, so in the rotated coordinates the two measurements carry no linear relationship, which is the point of the change of variable. Dropping $y_2$ now costs $20$ percent of the variance and reduces two measurements to one.

Practice

Both applications start by finding the eigenvalues of a symmetric matrix: the covariance matrix in the statistics case, $A^{T}A$ in the image case. The eigenvalues are also what the percentages are computed from.

The directions themselves are eigenvectors, and a principal component is a unit eigenvector, so normalizing is part of the answer.

Quiz

Five items on eigenvalues and eigenvectors of small symmetric matrices, the kind that appear as covariance matrices here.