In our discussion of factor analysis, we gave a way to model data $x\in {\mathbb{R}}^d$$x\in {\mathbb{R}}^d$x inR^(d)x \in \mathbb{R}^{d} as "approximately" lying in some $k$$k$kk-dimension subspace, where $k\ll d$$k\ll d$k≪dk \ll d. Specifically, we imagined that each point $x^{(i)}$$x^{(i)}$x^((i))x^{(i)} was created by first generating some $z^{(i)}$$z^{(i)}$z^((i))z^{(i)} lying in the $k$$k$kk-dimension affine space $\{\mathrm{\Lambda }z+\mu ;z\in {\mathbb{R}}^k\}$$\left\{\mathrm{\Lambda }z+\mu ;z\in {\mathbb{R}}^k\right\}${Lambda z+mu;z inR^(k)}\left\{\Lambda z+\mu ; z \in \mathbb{R}^{k}\right\}, and then adding $\mathrm{\Psi }$$\mathrm{\Psi }$Psi\Psi-covariance noise. Factor analysis is based on a probabilistic model, and parameter estimation used the iterative EM algorithm.

In this set of notes, we will develop a method, Principal Components Analysis (PCA), that also tries to identify the subspace in which the data approximately lies. However, PCA will do so more directly, and will require only an eigenvector calculation (easily done with the eig function in Matlab), and does not need to resort to EM.

Suppose we are given a dataset $\{x^{(i)};i=1,\dots ,n\}$$\left\{x^{(i)};i=1,\dots ,n\right\}${x^((i));i=1,dots,n}\left\{x^{(i)} ; i=1, \ldots, n\right\} of attributes of $n$$n$nn different types of automobiles, such as their maximum speed, turn radius, and so on. Let $x^{(i)}\in {\mathbb{R}}^d$$x^{(i)}\in {\mathbb{R}}^d$x^((i))inR^(d)x^{(i)} \in \mathbb{R}^{d} for each $i(d\ll n)$$i(d\ll n)$i(d≪n)i(d \ll n). But unknown to us, two different attributes - some $x_i$$x_i$x_(i)x_{i} and $x_j$$x_j$x_(j)x_{j}-respectively give a car's maximum speed measured in miles per hour, and the maximum speed measured in kilometers per hour. These two attributes are therefore almost linearly dependent, up to only small differences introduced by rounding off to the nearest mph or kph. Thus, the data really lies approximately on an $n-1$$n-1$n-1n-1 dimensional subspace. How can we automatically detect, and perhaps remove, this redundancy?

For a less contrived example, consider a dataset resulting from a survey of pilots for radio-controlled helicopters, where $x_1^{(i)}$$x_1^{(i)}$x_(1)^((i))x_{1}^{(i)} is a measure of the piloting skill of pilot $i$$i$ii, and $x_2^{(i)}$$x_2^{(i)}$x_(2)^((i))x_{2}^{(i)} captures how much he/she enjoys flying. Because RC helicopters are very difficult to fly, only the most committed students, ones that truly enjoy flying, become good pilots. So, the two attributes $x_1$$x_1$x_(1)x_{1} and $x_2$$x_2$x_(2)x_{2} are strongly correlated. Indeed, we might posit that that the data actually likes along some diagonal axis (the $u_1$$u_1$u_(1)u_{1} direction) capturing the intrinsic piloting "karma" of a person, with only a small amount of noise lying off this axis. (See figure.) How can we automatically compute this $u_1$$u_1$u_(1)u_{1} direction?

We will shortly develop the PCA algorithm. But prior to running PCA per se, typically we first preprocess the data by normalizing each feature to have mean 0 and variance 1. We do this by subtracting the mean and dividing by the empirical standard deviation:

where ${\mu }_j=\frac{1}{n}\sum _{i=1}^n{x}_j^{(i)}$${\mu }_j=\frac{1}{n}\sum _{i=1}^n x_j^{(i)}$mu_(j)=(1)/(n)sum_(i=1)^(n)x_(j)^((i))\mu_{j}=\frac{1}{n} \sum_{i=1}^{n} x_{j}^{(i)} and ${\sigma }_j^2=\frac{1}{n}\sum _{i=1}^n(x_j^{(i)}-{\mu }_j{)}^2$${\sigma }_j^2=\frac{1}{n}\sum _{i=1}^n (x_j^{(i)}-{\mu }_j{)}^2$sigma_(j)^(2)=(1)/(n)sum_(i=1)^(n)(x_(j)^((i))-mu_(j))^(2)\sigma_{j}^{2}=\frac{1}{n} \sum_{i=1}^{n}(x_{j}^{(i)}-\mu_{j})^{2} are the mean variance of feature $j$$j$jj, respectively.

Subtracting ${\mu }_j$${\mu }_j$mu_(j)\mu_{j} zeros out the mean and may be omitted for data known to have zero mean (for instance, time series corresponding to speech or other acoustic signals). Dividing by the standard deviation ${\sigma }_j$${\sigma }_j$sigma_(j)\sigma_{j} rescales each coordinate to have unit variance, which ensures that different attributes are all treated on the same "scale." For instance, if $x_1$$x_1$x_(1)x_{1} was cars' maximum speed in mph (taking values in the high tens or low hundreds) and $x_2$$x_2$x_(2)x_{2} were the number of seats (taking values around 2-4), then this renormalization rescales the different attributes to make them more comparable. This rescaling may be omitted if we had a priori knowledge that the different attributes are all on the same scale. One example of this is if each data point represented a grayscale image, and each $x_j^{(i)}$$x_j^{(i)}$x_(j)^((i))x_{j}^{(i)} took a value in $\{0,1,\dots ,255\}$$\{0,1,\dots ,255\}${0,1,dots,255}\{0,1, \ldots, 255\} corresponding to the intensity value of pixel $j$$j$jj in image $i$$i$ii.

Now, having normalized our data, how do we compute the "major axis of variation" $u$$u$uu–that is, the direction on which the data approximately lies? One way is to pose this problem as finding the unit vector $u$$u$uu so that when the data is projected onto the direction corresponding to $u$$u$uu, the variance of the projected data is maximized. Intuitively, the data starts off with some amount of variance/information in it. We would like to choose a direction $u$$u$uu so that if we were to approximate the data as lying in the direction/subspace corresponding to $u$$u$uu, as much as possible of this variance is still retained.

Consider the following dataset, on which we have already carried out the normalization steps:

Now, suppose we pick $u$$u$uu to correspond the the direction shown in the figure below. The circles denote the projections of the original data onto this line.

We see that the projected data still has a fairly large variance, and the points tend to be far from zero. In contrast, suppose had instead picked the following direction:

Here, the projections have a significantly smaller variance, and are much closer to the origin.

We would like to automatically select the direction $u$$u$uu corresponding to the first of the two figures shown above. To formalize this, note that given a unit vector $u$$u$uu and a point $x$$x$xx, the length of the projection of $x$$x$xx onto $u$$u$uu is given by $x^{T}u$$x^{T}u$x^(T)ux^{T} u. I.e., if $x^{(i)}$$x^{(i)}$x^((i))x^{(i)} is a point in our dataset (one of the crosses in the plot), then its projection onto $u$$u$uu (the corresponding circle in the figure) is distance $x^{T}u$$x^{T}u$x^(T)ux^{T} u from the origin. Hence, to maximize the variance of the projections, we would like to choose a unit-length $u$$u$uu so as to maximize:

We easily recognize that the maximizing this subject to $\parallel u{\parallel }_2=1$$\parallel u{\parallel }_2=1$||u||_(2)=1\|u\|_{2}=1 gives the principal eigenvector of $\mathrm{\Sigma }=\frac{1}{n}\sum _{i=1}^n{x}^{(i)}{x}^{(i{)}^T}$$\mathrm{\Sigma }=\frac{1}{n}\sum _{i=1}^n x^{(i)}{x}^{(i{)}^T}$Sigma=(1)/(n)sum_(i=1)^(n)x^((i))x^((i)^(T))\Sigma=\frac{1}{n} \sum_{i=1}^{n} x^{(i)} x^{(i)^{T}}, which is just the empirical covariance matrix of the data (assuming it has zero mean).^[1]

To summarize, we have found that if we wish to find a 1-dimensional subspace with with to approximate the data, we should choose $u$$u$uu to be the principal eigenvector of $\mathrm{\Sigma }$$\mathrm{\Sigma }$Sigma\Sigma. More generally, if we wish to project our data into a $k$$k$kk-dimensional subspace $(k<d)$$(k<d)$(k < d)(k<d), we should choose $u_1,\dots ,u_k$$u_1,\dots ,u_k$u_(1),dots,u_(k)u_{1}, \ldots, u_{k} to be the top $k$$k$kk eigenvectors of $\mathrm{\Sigma }$$\mathrm{\Sigma }$Sigma\Sigma. The $u_i$$u_i$u_(i)u_{i}'s now form a new, orthogonal basis for the data.^[2]

Then, to represent $x^{(i)}$$x^{(i)}$x^((i))x^{(i)} in this basis, we need only compute the corresponding vector

Thus, whereas $x^{(i)}\in {\mathbb{R}}^d$$x^{(i)}\in {\mathbb{R}}^d$x^((i))inR^(d)x^{(i)} \in \mathbb{R}^{d}, the vector $y^{(i)}$$y^{(i)}$y^((i))y^{(i)} now gives a lower, $k$$k$kk-dimensional, approximation/representation for $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}. PCA is therefore also referred to as a dimensionality reduction algorithm. The vectors $u_1,\dots ,u_k$$u_1,\dots ,u_k$u_(1),dots,u_(k)u_{1}, \ldots, u_{k} are called the first $k$$k$kk principal components of the data.

Remark. Although we have shown it formally only for the case of $k=1$$k=1$k=1k=1, using well-known properties of eigenvectors it is straightforward to show that of all possible orthogonal bases $u_1,\dots ,u_k$$u_1,\dots ,u_k$u_(1),dots,u_(k)u_{1}, \ldots, u_{k}, the one that we have chosen maximizes $\sum _i{\parallel y^{(i)}\parallel }_2^2$$\sum _i {∥y^{(i)}∥}_2^2$sum_(i)||y^((i))||_(2)^(2)\sum_{i}\left\|y^{(i)}\right\|_{2}^{2}. Thus, our choice of a basis preserves as much variability as possible in the original data.

In problem set $4$$4$44, you will see that PCA can also be derived by picking the basis that minimizes the approximation error arising from projecting the data onto the $k$$k$kk-dimensional subspace spanned by them.

PCA has many applications; we will close our discussion with a few examples. First, compression-representing $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s with lower dimension $y^{(i)}$$y^{(i)}$y^((i))y^{(i)}'s–is an obvious application. If we reduce high dimensional data to $k=2$$k=2$k=2k=2 or $3$$3$33 dimensions, then we can also plot the $y^{(i)}$$y^{(i)}$y^((i))y^{(i)}'s to visualize the data. For instance, if we were to reduce our automobiles data to $2$$2$22 dimensions, then we can plot it (one point in our plot would correspond to one car type, say) to see what cars are similar to each other and what groups of cars may cluster together.

Another standard application is to preprocess a dataset to reduce its dimension before running a supervised learning learning algorithm with the $x^{(i)}$$x^{(i)}$x^((i))x^{(i)}'s as inputs. Apart from computational benefits, reducing the data's dimension can also reduce the complexity of the hypothesis class considered and help avoid overfitting (e.g., linear classifiers over lower dimensional input spaces will have smaller VC dimension).

Lastly, as in our RC pilot example, we can also view PCA as a noise reduction algorithm. In our example it, estimates the intrinsic "piloting karma" from the noisy measures of piloting skill and enjoyment. In class, we also saw the application of this idea to face images, resulting in eigenfaces method. Here, each point $x^{(i)}\in {\mathbb{R}}^{100\times 100}$$x^{(i)}\in {\mathbb{R}}^{100\times 100}$x^((i))inR^(100 xx100)x^{(i)} \in \mathbb{R}^{100 \times 100} was a $10000$$10000$1000010000 dimensional vector, with each coordinate corresponding to a pixel intensity value in a $100\times 100$$100\times 100$100 xx100100 \times 100 image of a face. Using PCA, we represent each image $x^{(i)}$$x^{(i)}$x^((i))x^{(i)} with a much lower-dimensional $y^{(i)}$$y^{(i)}$y^((i))y^{(i)}. In doing so, we hope that the principal components we found retain the interesting, systematic variations between faces that capture what a person really looks like, but not the "noise" in the images introduced by minor lighting variations, slightly different imaging conditions, and so on. We then measure distances between faces $i$$i$ii and $j$$j$jj by working in the reduced dimension, and computing ${\parallel y^{(i)}-y^{(j)}\parallel }_2$${∥y^{(i)}-y^{(j)}∥}_2$||y^((i))-y^((j))||_(2)\left\|y^{(i)}-y^{(j)}\right\|_{2}. This resulted in a surprisingly good face-matching and retrieval algorithm.