Principle Component Analysis

1. Notation

• Assume the data consists of $n$ datapoints of $p$ variables. Let $\mathbf{X}=(x_{ij})\in {\mathbb{R}}^{n\times p}$ with $x_{ij}$ being the $j$th variable for the $i$th datapoint.

• Denote the $i$th datapoint by ${\mathbf{x}}_i\in {\mathbb{R}}^p$ and assume ${\mathbf{x}}_1,\dots ,{\mathbf{x}}_n$ are i.i.d. samples of a random vector $X$ over ${\mathbb{R}}^p$. The $j$th dimension of $X$ will be denoted $X^{(j)}$.

2. Data Representation

• A key class of unsupervised methods are converting (input) data into alternative forms by learning some mapping $\mathbf{x}\mapsto \mathbf{z}$, where $\mathbf{x}\in {\mathbb{R}}^p$ and $\mathbf{z}\in {\mathbb{R}}^k$.

• Sometimes we map to a higher dimensional space ($k>p$) to introduce additional features (e.g., kernel method).

• More commonly, we perform dimensionality reduction by mapping to a lower dimensional space ($k<p$, typically $k\ll p$).

3. Principle Component Analysis (PCA)

PCA is a linear dimensionality reduction technique, and it finds a new basis for the data $${\mathbf{z}}_i=V^{\mathrm{\top }}{\mathbf{x}}_i.$$

• We want to choose a $V^{\mathrm{\top }}$ that captures the most information about variability in $X$.

• PCA uses variance as a proxy: it finds an orthogonal basis that maximizes variance of $Z=V^{\mathrm{\top }}X$.

• PCA is often a useful data preprocessing step (particularly with high-dimensional data) or for exploratory data analysis.

3.1. Definition

PCA is a technique to find an orthogonal basis $\{{\mathbf{v}}_1,\dots ,{\mathbf{v}}_p\}$ for the data space s.t.

• the first principal component (PC) ${\mathbf{v}}_1$ is the direction of greatest variance of data: $${\mathbf{v}}_1=\underset{\mathbf{v}:{\mathbf{v}}^{\mathrm{\top }}\mathbf{v}=1}{\arg max}\text{Var}[{\mathbf{v}}^{\mathrm{\top }}X];$$

• the $j$th PC ${\mathbf{v}}_j$ is the direction orthogonal to ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_{j-1}$ of greatest variance for $j=2,\dots ,p$: $${\mathbf{v}}_j=\underset{\mathbf{v}:{\mathbf{v}}^{\mathrm{\top }}\mathbf{v}=1\cap {\mathbf{v}}^{\mathrm{\top }}{\mathbf{v}}_i=0,\mathrm{\forall }i<j}{\arg max}\text{Var}[{\mathbf{v}}^{\mathrm{\top }}X].$$

3.2. Derivation

Recall that ${\displaystyle {\mathbf{v}}_1=\underset{\mathbf{v}}{\arg max}\text{Var}[{\mathbf{v}}^{\mathrm{\top }}X]}$ s.t. ${\mathbf{v}}^{\mathrm{\top }}\mathbf{v}=1$. Define $$\begin{array}{rl}\mathcal{L}({\mathbf{v}}_1,\lambda )& =\text{Var}[{\mathbf{v}}^{\mathrm{\top }}X]-\lambda ({\mathbf{v}}^{\mathrm{\top }}\mathbf{v}-1)\\ & ={\mathbf{v}}^{\mathrm{\top }}\text{Var}[X]\mathbf{v}-\lambda ({\mathbf{v}}^{\mathrm{\top }}\mathbf{v}-1)\\ & \approx {\mathbf{v}}^{\mathrm{\top }}S\mathbf{v}-\lambda ({\mathbf{v}}^{\mathrm{\top }}\mathbf{v}-1)\end{array}$$ where ${\displaystyle S=\frac{1}{n-1}\sum _{i=1}^n{\mathbf{x}}_i{\mathbf{x}}_i^{\mathrm{\top }}}$ is the sample covariance matrix. Let $$\frac{\partial }{\partial \mathbf{v}}\mathcal{L}({\mathbf{v}}_1,\lambda )=2S{\mathbf{v}}_1-2\lambda {\mathbf{v}}_1=0$$ and we have $S{\mathbf{v}}_1=\lambda {\mathbf{v}}_1$, i.e., ${\mathbf{v}}_1$ is the eigenvector of $S$.

Since ${\mathbf{v}}^{\mathrm{\top }}S\mathbf{v}={\mathbf{v}}^{\mathrm{\top }}\lambda \mathbf{v}=\lambda$, then if we want to maximize, we want eigenvector corresponding to the largest eigenvalue. Similarly, ${\mathbf{v}}_j$ is the eigenvector of $S$ with $j$th largest eigenvalue.

The PCA components can be found by the eigendecomposition $$S=\frac{1}{n-1}{\mathbf{X}}^{\mathrm{\top }}\mathbf{X}=V\mathrm{\Lambda }{V}^{\mathrm{\top }}$$ where $\mathrm{\Lambda }$ is a diagonal matrix with eigenvalues (variances along each PC) ${\lambda }_1\ge \cdots \ge {\lambda }_p\ge 0$, and $V$ (loadings matrix) is a $p\times p$ orthogonal matrix whose columns are the $p$ eigenvectors of $S$, i.e., the PCs ${\mathbf{v}}_1,\dots ,{\mathbf{v}}_p$.

3.3. Optimal Linear Reconstruction

• The $k$-dimensional representation of data item ${\mathbf{x}}_i$ is the vector of projections of ${\mathbf{x}}_i$ onto first $k$ PCs: $${\mathbf{z}}_i=V_{1:k}^{\mathrm{\top }}{\mathbf{x}}_i=[{\mathbf{v}}_i^{\mathrm{\top }}{\mathbf{x}}_i,\dots ,{\mathbf{v}}_k^{\mathrm{\top }}{\mathbf{x}}_i{]}^{\mathrm{\top }}\in {\mathbb{R}}^k$$ where $V_{1:k}=\left[\begin{array}{ccc}{\mathbf{v}}_1& \cdots & {\mathbf{v}}_k\end{array}\right]$.

• Transformed data matrix, a.k.a. scores matrix: $$\mathbf{Z}=\mathbf{X}{V}_{1:k}\in {\mathbb{R}}^{n\times k}.$$

• Reconstruction of ${\mathbf{x}}_i$: $${\hat{\mathbf{x}}}_i=V_{1:k}{V}_{1:k}^{\mathrm{\top }}{\mathbf{x}}_i.$$

• PCA gives the optimal linear reconstruction of the original data based on a $k$-dimensional compression.

3.4. Property

• $S$ is a real symmetric matrix, and thus eigenvectors (PCs) are orthogonal.

• Projection to the $j$th PC, $Z^{(j)}={\mathbf{v}}_j^{\mathrm{\top }}X$, has sample variance ${\lambda }_j$, since $$\widehat{\text{Var}}[Z^{(j)}]={\mathbf{v}}_j^{\mathrm{\top }}S{\mathbf{v}}_j={\lambda }_j{\mathbf{v}}_j^{\mathrm{\top }}{\mathbf{v}}_j={\lambda }_j.$$

• Projections to PCs are uncorrelated, since for $i\ne j$, $$\widehat{\text{Cov}}(Z^{(i)},Z^{(j)})={\mathbf{v}}_i^{\mathrm{\top }}S{\mathbf{v}}_j={\lambda }_j{\mathbf{v}}_i^{\mathrm{\top }}{\mathbf{v}}_j=0.$$

• The total sample variance is given by $$\text{Tr}(S)=\sum _{i=1}^p{S}_{ii}={\lambda }_1+\cdots +{\lambda }_p$$ so the proportion of total variance explained by the $j$th PC is ${\displaystyle \frac{{\lambda }_j}{{\lambda }_1+\cdots +{\lambda }_p}}$.

4. Singular Value Decomposition (SVD)

• Recall that any real-valued $n\times p$ matrix $\mathbf{X}$ can be written as $\mathbf{X}=UD{V}^{\mathrm{\top }}$, where $U$ is an $n\times n$ orthogonal matrix (i.e., $U{U}^{\mathrm{\top }}=U^{\mathrm{\top }}U=I_n$), $D$ is an $n\times p$ matrix with decreasing non-negative elements on the diagonal (the singular values) and zero off-diagonal elements, and $V$ is a $p\times p$ orthogonal matrix (i.e., $V{V}^{\mathrm{\top }}=V^{\mathrm{\top }}V=I_p$).

• SVD always exists even for non-square matrices.

• Fast and numerically stable algorithms for SVD are available in most languages.

Let $\mathbf{X}=UD{V}^{\mathrm{\top }}$ be the SVD of the $n\times p$ data matrix $\mathbf{X}$. We have $$(n-1)S={\mathbf{X}}^{\mathrm{\top }}\mathbf{X}=(UD{V}^{\mathrm{\top }}{)}^{\mathrm{\top }}(UD{V}^{\mathrm{\top }})=V{D}^{\mathrm{\top }}{U}^{\mathrm{\top }}UD{V}^{\mathrm{\top }}=V{D}^{\mathrm{\top }}D{V}^{\mathrm{\top }}.$$

By analogy to the eigendecomposition $S=V\mathrm{\Lambda }{V}^{\mathrm{\top }}$:

• $V$ from SVD of $\mathbf{X}$ is the same as from eigendecomposition of $S$, and thus the columns in $V$ are the PCs;

• ${\displaystyle \mathrm{\Lambda }=\frac{1}{n-1}{D}^{\mathrm{\top }}D}$ where $D_{ii}=\sqrt{(n-1){\lambda }_i}$ for all $i\le min(n,p)$, and if $p>n$, ${\lambda }_i=0$ for all $i>n$.

5. PCA Projection from Gram Matrix

• Recall that $\mathbf{K}=\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}$, where ${\mathbf{K}}_{ij}={\mathbf{x}}_i^{\mathrm{\top }}{\mathbf{x}}_j$ is called the Gram matrix of dataset $\mathbf{X}$.

• $\mathbf{K}$ and $(n-1)S={\mathbf{X}}^{\mathrm{\top }}\mathbf{X}$ have the same nonzero eigenvalues, equal to the nonzero squared singular values of $\mathbf{X}$.

We have $$\mathbf{X}{\mathbf{X}}^{\mathrm{\top }}=(UD{V}^{\mathrm{\top }})(UD{V}^{\mathrm{\top }}{)}^{\mathrm{\top }}=UD{V}^{\mathrm{\top }}V{D}^{\mathrm{\top }}{U}^{\mathrm{\top }}=UD{D}^{\mathrm{\top }}{U}^{\mathrm{\top }}.$$ Note that the columns of $U$ represent the eigenvectors of $\mathbf{K}$, but they are different from the eigenvectors of ${\mathbf{X}}^{\mathrm{\top }}\mathbf{X}$.

If we consider projection to all PCs, the transformed data matrix is $$\mathbf{Z}=\mathbf{X}V=UD{V}^{\mathrm{\top }}V=UD.$$ Hence, $\mathbf{Z}$ can be obtained from the eigendecomposition of Gram matrix $\mathbf{K}$, where eigendecomposition gives us $U$ directly, and $D_{ii}=\sqrt{(n-1){\lambda }_i}$. When $p\gg n$, eigendecomposition of $\mathbf{K}$ requires much less computation ($\mathcal{O}(n^3)$) than the eigendecomposition of the covariance matrix ($\mathcal{O}(p^3)$).