---
title: "Applications: Image Processing and Statistics"
module: Symmetric Matrices, Quadratic Forms, and the SVD
moduleNumber: 7
lessonNumber: 5
order: 705
summary: >
  Principal component analysis diagonalizes the covariance matrix of a data set,
  producing uncorrelated variables ordered by variance. The leading components
  capture most of the variation, which reduces dimension, compresses images
  through low-rank SVD approximation, and connects directly to the singular
  values of the data matrix.
topics: [Symmetric Matrices, Quadratic Forms, and the SVD]
sources:
  - book: Lay
    ref: "Ch. 7 — Symmetric Matrices and Quadratic Forms; §7.5 Applications to Image Processing and Statistics"
draft: false
---

Principal component analysis applies the
[singular value decomposition](/linear-algebra/symmetric-quadratic-svd/singular-value-decomposition)
and the [Spectral Theorem](/linear-algebra/symmetric-quadratic-svd/diagonalizing-symmetric-matrices)
to multivariate data: a collection of objects, each described by several
measurements, so each datum is a vector in $\mathbb{R}^p$. PCA finds new
coordinates in which the measurements are uncorrelated and sorted by decreasing
variance, then discards the coordinates with little variance.

For example, a chemical plant might test each of $300$ plastic samples on eight
properties: melting point, density, tensile strength, and so on. Each report is
a vector in $\mathbb{R}^8$, and the set forms an $8 \times 300$ matrix of
observations. A satellite image measured in three spectral bands assigns each
pixel a vector in $\mathbb{R}^3$; a two-megapixel image is then a $3 \times
(4{,}000{,}000)$ matrix. The dimension that PCA reduces is the number of
measurements per object, not the number of objects.

## Mean, deviation, and covariance

Let $[\,X_1 \; \cdots \; X_N\,]$ be a $p \times N$ matrix of observations. The
**sample mean** is the average vector

$$
M = \frac{1}{N}(X_1 + \cdots + X_N),
$$

the center of the data cloud. Subtracting it recenters every observation:
$\hat{X}_k = X_k - M$, and the matrix $B = [\,\hat{X}_1 \; \cdots \;
\hat{X}_N\,]$ is said to be in **mean-deviation form**, with zero mean.

> **Definition (Sample covariance matrix).** For data in mean-deviation form
> $B$, the $p \times p$ matrix
> $$
> S = \frac{1}{N-1}\, B B^\top.
> $$

Because $S$ has the form $\tfrac{1}{N-1}BB^\top$, it is symmetric and positive
semidefinite: $x^\top S x = \tfrac{1}{N-1}\|B^\top x\|^2 \geq 0$. Its entries
carry the spread and coupling of the measurements:

- **Variances** on the diagonal. The entry $s_{jj}$ is the variance of the
  $j$-th measurement — how widely the $j$-th coordinate ranges over the data.
- **Covariances** off the diagonal. The entry $s_{ij}$ ($i \neq j$) is the
  covariance of measurements $i$ and $j$. When $s_{ij} = 0$ the two are
  **uncorrelated**.
- **Total variance** is the trace $\operatorname{tr}(S) = s_{11} + \cdots +
  s_{pp}$, the sum of the individual variances.

> **Worked example.** Three measurements on four individuals give observation
> vectors $X_1 = (1,2,1)$, $X_2 = (4,2,13)$, $X_3 = (7,8,1)$, $X_4 = (8,4,5)$.
> The mean is $M = (5, 4, 5)$, so the mean-deviation matrix and covariance are
> $$
> B = \begin{bmatrix} -4 & -1 & 2 & 3 \\ -2 & -2 & 4 & 0 \\ -4 & 8 & -4 & 0 \end{bmatrix},
> \qquad
> S = \frac{1}{3} B B^\top = \begin{bmatrix} 10 & 6 & 0 \\ 6 & 8 & -8 \\ 0 & -8 & 32 \end{bmatrix}.
> $$
> The variance of the third measurement is $32$, much larger than the first
> measurement's $10$, so the third coordinate spreads more widely. The $(1,3)$
> entry is $0$: the first and third measurements are uncorrelated.

Analysis is easiest when the covariance matrix is diagonal, meaning all
measurements are mutually uncorrelated; PCA produces exactly this.

## Principal component analysis

PCA looks for an orthogonal change of variable that decorrelates the data. With
the observations in mean-deviation form, find an orthogonal $P = [\,u_1 \;
\cdots \; u_p\,]$ giving new coordinates $X = PY$, equivalently $Y = P^\top X$,
so that the new variables $y_1, \dots, y_p$ are uncorrelated and ordered by
decreasing variance. The covariance matrix of the transformed data is $P^\top S
P$, and making it diagonal is orthogonal diagonalization of $S$.

> **Definition (Principal components).** For a covariance matrix $S = PDP^\top$
> with eigenvalues $\lambda_1 \geq \cdots \geq \lambda_p \geq 0$ and orthonormal
> eigenvectors $u_1, \dots, u_p$, the eigenvectors are the **principal
> components** of the data. The first principal component $u_1$ points along the
> direction of greatest variance.

Since $S$ is symmetric, the Spectral Theorem supplies the diagonalizing $P$.
The eigenvalue $\lambda_j$ is the variance of the new variable $y_j$, and the
$j$-th principal component defines that variable as a weighted combination of the
originals: if $u_1$ has entries $c_1, \dots, c_p$, then

$$
y_1 = u_1^\top X = c_1 x_1 + c_2 x_2 + \cdots + c_p x_p.
$$

The weights are the entries of the eigenvector, and $y_1$ has the largest
variance any unit-weighted combination can achieve. This last claim is the extreme-values-on-the-unit-sphere theorem
of the [constrained optimization](/linear-algebra/symmetric-quadratic-svd/constrained-optimization)
lesson: for a unit vector $u$, the variance of $y = u^\top X$ is $u^\top S u$,
whose maximum over unit $u$ is the top eigenvalue $\lambda_1$, attained at $u_1$.
Each later component maximizes variance subject to being uncorrelated with the
ones before.

$$
% caption: The first principal component u₁ lies along the long axis of the data
% cloud — the direction of greatest variance — and u₂ is perpendicular; the mean
% M is the center from which deviations are measured.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  % scatter of points along a tilted line
  \foreach \p in {(-2.4,-1.1),(-1.7,-0.6),(-1.9,-1.0),(-1.0,-0.6),(-0.6,-0.1),
                  (-0.2,-0.4),(0.2,0.3),(0.7,0.2),(1.1,0.7),(1.6,0.6),
                  (1.9,1.1),(2.4,1.0),(0.4,-0.1),(-1.3,-0.2),(1.3,0.3)}
    \fill[acc!70] \p circle (1.7pt);
  % mean
  \fill[black] (0,0) circle (2.2pt);
  \node[black] at (0.25,-0.45) {$M$};
  % PC1 (long axis, slope ~0.44)
  \draw[->, acc, very thick] (0,0) -- (2.7,1.2) node[right, black] {$u_1$};
  \draw[acc, thick, dashed] (0,0) -- (-2.7,-1.2);
  % PC2 (perpendicular, short)
  \draw[->, acc, very thick] (0,0) -- (-0.55,1.24) node[above, black] {$u_2$};
\end{tikzpicture}
$$

### Multispectral imaging

> **Worked example.** For a three-band satellite image of Railroad Valley,
> Nevada, each pixel is a vector of three signal intensities, and the covariance
> matrix over four million pixels is
> $$
> S = \begin{bmatrix} 2382.78 & 2611.84 & 2136.20 \\ 2611.84 & 3106.47 & 2553.90 \\ 2136.20 & 2553.90 & 2650.71 \end{bmatrix}.
> $$
> Its eigenvalues and unit eigenvectors are
> $$
> \lambda_1 = 7614.23, \; u_1 = \begin{bmatrix} .54 \\ .63 \\ .56 \end{bmatrix};
> \quad
> \lambda_2 = 427.63, \; u_2 = \begin{bmatrix} .49 \\ .30 \\ -.82 \end{bmatrix};
> \quad
> \lambda_3 = 98.10, \; u_3 = \begin{bmatrix} -.68 \\ .72 \\ -.14 \end{bmatrix}.
> $$
> The first principal component defines the new variable
> $y_1 = .54 x_1 + .63 x_2 + .56 x_3$, a near-equal blend of the three bands. An
> image whose pixels are shaded by $y_1$ carries most of the scene's information
> in a single grayscale channel.

## Reducing dimension

An orthogonal change of variable preserves total variance, because $P$ does not
change lengths: $\operatorname{tr}(D) = \operatorname{tr}(P^\top S P) =
\operatorname{tr}(S)$. So the variance is redistributed across the new
variables,

$$
\underbrace{\lambda_1 + \cdots + \lambda_p}_{\text{new variables}}
= \operatorname{tr}(S) = \underbrace{s_{11} + \cdots + s_{pp}}_{\text{original variables}},
$$

and the fraction of total variance captured by $y_j$ is $\lambda_j /
\operatorname{tr}(S)$. When a few eigenvalues dominate, the corresponding
components describe most of the data.

For the Railroad Valley image, the total variance is $\operatorname{tr}(S) =
7614.23 + 427.63 + 98.10 = 8139.96$, split as

| Component | Eigenvalue | Fraction of variance |
| --- | --- | --- |
| First | $7614.23$ | $93.5\%$ |
| Second | $427.63$ | $5.3\%$ |
| Third | $98.10$ | $1.2\%$ |

The first component alone holds $93.5\%$ of the variation, so a single grayscale
image displays almost all the information the three bands carry. The third
component is nearly constant across the scene; the data effectively lie in a
plane, and largely along one line within it. Three-dimensional data has been
compressed to essentially one dimension with a $1.2\%$ loss.

$$
% caption: The variance spectrum of the covariance matrix falls off sharply;
% the cumulative curve shows the first component already explains most of the
% total, so truncating after it discards little.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[black, ->] (-0.2,0) -- (6.4,0) node[right, black] {component};
  \draw[black, ->] (-0.2,0) -- (-0.2,3.2) node[above, black] {share};
  % bars
  \fill[acc!30, draw=acc] (0.4,0) rectangle (1.4,2.8);
  \fill[acc!20, draw=acc] (2.0,0) rectangle (3.0,0.48);
  \fill[acc!15, draw=acc] (3.6,0) rectangle (4.6,0.12);
  \node[black] at (0.9,3.0) {$93.5\%$};
  \node[black] at (2.5,0.75) {$5.3\%$};
  \node[black] at (4.1,0.4) {$1.2\%$};
  % cumulative curve
  \draw[black, thick] (0.9,2.8) -- (2.5,3.05) -- (4.1,3.1);
  \fill[black] (0.9,2.8) circle (1.6pt) (2.5,3.05) circle (1.6pt) (4.1,3.1) circle (1.6pt);
  \node[black] at (5.4,3.1) {cumulative};
\end{tikzpicture}
$$

## The SVD route

Diagonalizing $S$ directly is not the preferred computation. Forming $S =
\tfrac{1}{N-1}BB^\top$ squares the data, magnifying errors, the same hazard
noted for $A^\top A$ in the SVD lesson. Instead, run the SVD on the data matrix
itself. Set

$$
A = \frac{1}{\sqrt{N-1}}\, B^\top,
$$

an $N \times p$ matrix. Then $A^\top A = \tfrac{1}{N-1} B B^\top = S$, so:

- The **squares of the singular values** of $A$ are the eigenvalues of $S$,
  hence the variances $\lambda_j$ of the principal components.
- The **right singular vectors** of $A$ are the eigenvectors of $A^\top A = S$,
  hence the principal components themselves.

Computing the SVD of $A$ iteratively is faster and more accurate than an
eigenvalue decomposition of $S$, and the gap widens with dimension. For a
hyperspectral image with $p = 224$ bands, PCA runs in seconds this way. The
[numerical SVD](/linear-algebra/numerical-linear-algebra/numerical-eigenvalues-and-svd)
lesson develops the algorithm.

## Low-rank approximation and compression

The [rank-one expansion](/linear-algebra/symmetric-quadratic-svd/singular-value-decomposition)
$A = \sum_{i=1}^r \sigma_i u_i v_i^\top$ orders the pieces of $A$ by singular
value, so truncating after $k$ terms,

$$
A_k = \sigma_1 u_1 v_1^\top + \cdots + \sigma_k u_k v_k^\top,
$$

gives the best rank-$k$ approximation of $A$: no other rank-$k$ matrix is closer
in either the operator or Frobenius norm. Image compression uses this truncation
directly. A grayscale image is an $m \times n$ matrix of pixel intensities;
storing it in full takes $mn$ numbers. Storing $A_k$ takes only $k(m + n + 1)$:
the $k$ singular values, plus $k$ vectors of length $m$ and $k$ of length $n$.
When the singular values decay quickly, a small $k$ reconstructs the image with
little visible loss at a fraction of the storage.

$$
% caption: A rank-k SVD truncation stores k(m+n+1) numbers instead of the full
% mn; when the singular values decay fast, small k reconstructs the image with
% little loss.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize,
  bar/.style={draw}]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[black, ->] (-0.2,0) -- (-0.2,2.9) node[above, black] {storage};
  % full
  \fill[black, draw=black] (0.5,0) rectangle (2.0,2.5);
  \node[black] at (1.25,2.75) {$mn$};
  \node[black] at (1.25,-0.4) {full image};
  % rank-k
  \fill[acc!20, draw=acc] (3.5,0) rectangle (5.0,0.7);
  \node[black] at (4.25,0.95) {$k(m{+}n{+}1)$};
  \node[black] at (4.25,-0.4) {rank-$k$};
  % arrow of savings
  \draw[->, acc, thick] (2.3,1.6) -- (3.3,0.9) node[midway, above, sloped, black] {compress};
\end{tikzpicture}
$$

The choice of $k$ is read off the singular value spectrum: keep enough terms to
capture a target fraction of the total energy $\sum \sigma_i^2$, and drop the
tail where the singular values are negligible. The same truncation denoises,
since noise usually spreads across the small singular values that get discarded.

The ratio $\sigma_1 / \sigma_r$ of largest to smallest retained singular value
is the [condition number](/linear-algebra/numerical-linear-algebra/conditioning-and-floating-point),
a measure of how close the matrix sits to rank-deficiency and how sensitive its
inverse problems are. PCA itself is equivalent to orthogonal regression: the
line through the mean along $u_1$ minimizes the sum of squared perpendicular
distances to the data.

$$
% caption: In the u₁ coordinate the first principal component is the horizontal
% axis; PCA places it to minimize the total squared perpendicular distance —
% the vertical drops — from the data to the line.
\begin{tikzpicture}[scale=1.0, >=stealth, font=\footnotesize]
  \definecolor{acc}{HTML}{4A6FA5}
  \draw[acc, very thick] (-3,0) -- (3.3,0) node[right, black] {$u_1$ line};
  \foreach \x/\y in {-2.6/0.42,-2.0/-0.32,-1.4/0.26,-0.8/-0.36,-0.2/0.22,
                     0.4/-0.26,1.0/0.32,1.6/-0.22,2.2/0.38,2.7/-0.3}{
    \draw[black, dashed] (\x,\y) -- (\x,0);
    \fill[acc!75] (\x,\y) circle (1.7pt);
  }
  \node[black] at (0,-1.05) {vertical drops = perpendicular distances};
\end{tikzpicture}
$$