What is the objective function of PCA?

Question

Principal component analysis can use matrix decomposition, but that is just a tool to get there.

How would you find the principal components without the use of matrix algebra?

What is the objective function (goal), and what are the constraints?

Answer 1

Without trying to give a full primer on PCA, from an optimization standpoint, the primary objective function is the Rayleigh quotient. The matrix that figures in the quotient is (some multiple of) the sample covariance matrix $$\newcommand{\m}[1]{\mathbf{#1}}\newcommand{\x}{\m{x}}\newcommand{\S}{\m{S}}\newcommand{\u}{\m{u}}\newcommand{\reals}{\mathbb{R}}\newcommand{\Q}{\m{Q}}\newcommand{\L}{\boldsymbol{\Lambda}} \S = \frac{1}{n} \sum_{i=1}^n \x_i \x_i^T = \m{X}^T \m{X} / n $$ where each $\x_i$ is a vector of $p$ features and $\m{X}$ is the matrix such that the $i$th row is $\x_i^T$.

PCA seeks to solve a sequence of optimization problems. The first in the sequence is the unconstrained problem $$ \begin{array}{ll} \text{maximize} & \frac{\u^T \S \u}{\u^T\u} \;, \u \in \reals^p \> . \end{array} $$

Since $\u^T \u = \|\u\|_2^2 = \|\u\| \|\u\|$, the above unconstrained problem is equivalent to the constrained problem $$ \begin{array}{ll} \text{maximize} & \u^T \S \u \\ \text{subject to} & \u^T \u = 1 \>. \end{array} $$

Here is where the matrix algebra comes in. Since $\S$ is a symmetric positive semidefinite matrix (by construction!) it has an eigenvalue decomposition of the form $$ \S = \Q \L \Q^T \>, $$ where $\Q$ is an orthogonal matrix (so $\Q \Q^T = \m{I}$) and $\L$ is a diagonal matrix with nonnegative entries $\lambda_i$ such that $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p \geq 0$.

Hence, $\u^T \S \u = \u^T \Q \L \Q^T \u = \m{w}^T \L \m{w} = \sum_{i=1}^p \lambda_i w_i^2$. Since $\u$ is constrained in the problem to have a norm of one, then so is $\m{w}$ since $\|\m{w}\|_2 = \|\Q^T \u\|_2 = \|\u\|_2 = 1$, by virtue of $\Q$ being orthogonal.

But, if we want to maximize the quantity $\sum_{i=1}^p \lambda_i w_i^2$ under the constraints that $\sum_{i=1}^p w_i^2 = 1$, then the best we can do is to set $\m{w} = \m{e}_1$, that is, $w_1 = 1$ and $w_i = 0$ for $i > 1$.

Now, backing out the corresponding $\u$, which is what we sought in the first place, we get that $$ \u^\star = \Q \m{e}_1 = \m{q}_1 $$ where $\m{q}_1$ denotes the first column of $\Q$, i.e., the eigenvector corresponding to the largest eigenvalue of $\S$. The value of the objective function is then also easily seen to be $\lambda_1$.

---

The remaining principal component vectors are then found by solving the sequence (indexed by $i$) of optimization problems $$ \begin{array}{ll} \text{maximize} & \u_i^T \S \u_i \\ \text{subject to} & \u_i^T \u_i = 1 \\ & \u_i^T \u_j = 0 \quad \forall 1 \leq j < i\>. \end{array} $$ So, the problem is the same, except that we add the additional constraint that the solution must be orthogonal to all of the previous solutions in the sequence. It is not difficult to extend the argument above inductively to show that the solution of the $i$th problem is, indeed, $\m{q}_i$, the $i$th eigenvector of $\S$.

---

The PCA solution is also often expressed in terms of the singular value decomposition of $\m{X}$. To see why, let $\m{X} = \m{U} \m{D} \m{V}^T$. Then $n \S = \m{X}^T \m{X} = \m{V} \m{D}^2 \m{V}^T$ and so $\m{V} = \m{Q}$ (strictly speaking, up to sign flips) and $\L = \m{D}^2 / n$.

The principal components are found by projecting $\m{X}$ onto the principal component vectors. From the SVD formulation just given, it's easy to see that $$ \m{X} \m{Q} = \m{X} \m{V} = \m{U} \m{D} \m{V}^T \m{V} = \m{U} \m{D} \> . $$

The simplicity of representation of both the principal component vectors and the principal components themselves in terms of the SVD of the matrix of features is one reason the SVD features so prominently in some treatments of PCA.

Answer 2

The solution presented by cardinal focuses on the sample covariance matrix. Another starting point is the reconstruction error of the data by a q-dimensional hyperplane. If the p-dimensional data points are $x_1, \ldots, x_n$ the objective is to solve

$$\min_{\mu, \lambda_1,\ldots, \lambda_n, \mathbf{V}_q} \sum_{i=1}^n ||x_i - \mu - \mathbf{V}_q \lambda_i||^2$$

for a $p \times q$ matrix $\mathbf{V}_q$ with orthonormal columns and $\lambda_i \in \mathbb{R}^q$. This gives the best rank q-reconstruction as measured by the euclidean norm, and the columns of the $\mathbf{V}_q$ solution are the first q principal component vectors.

For fixed $\mathbf{V}_q$ the solution for $\mu$ and $\lambda_i$ (this is regression) are $$\mu = \overline{x} = \frac{1}{n}\sum_{i=1}^n x_i \qquad \lambda_i = \mathbf{V}_q^T(x_i - \overline{x})$$

For ease of notation lets assume that $x_i$ have been centered in the following computations. We then have to minimize

$$\sum_{i=1}^n ||x_i - \mathbf{V}_q\mathbf{V}_q^T x_i||^2$$

over $\mathbf{V}_q$ with orthonormal columns. Note that $P = \mathbf{V}_q\mathbf{V}_q^T$ is the projection onto the q-dimensional column space. Hence the problem is equivalent to minimizing
$$\sum_{i=1}^n ||x_i - P x_i||^2 = \sum_{i=1}^n ||x_i||^2 - \sum_{i=1}^n||Px_i||^2$$ over rank q projections $P$. That is, we need to maximize $$\sum_{i=1}^n||Px_i||^2 = \sum_{i=1}^n x_i^TPx_i = \text{tr}(P \sum_{i=1}^n x_i x_i^T) = n \text{tr}(P \mathbf{S})$$ over rank q projections $P$, where $\mathbf{S}$ is the sample covariance matrix. Now $$\text{tr}(P\mathbf{S}) = \text{tr}(\mathbf{V}_q^T\mathbf{S}\mathbf{V}_q) = \sum_{i=1}^q u_i^T \mathbf{S} u_i$$ where $u_1, \ldots, u_q$ are the $q$ (orthonormal) columns in $\mathbf{V}_q$, and the arguments presented in @cardinal's answer show that the maximum is obtained by taking the $u_i$'s to be $q$ eigenvectors for $\mathbf{S}$ with the $q$ largest eigenvalues.

The reconstruction error suggests a number of useful generalizations, for instance sparse principal components or reconstructions by low-dimensional manifolds instead of hyperplanes. For details, see Section 14.5 in The Elements of Statistical Learning.

Answer 3

See NIPALS (wiki) for one algorithm which doesn't explicitly use a matrix decomposition. I suppose that's what you mean when you say that you want to avoid matrix algebra since you really can't avoid matrix algebra here :)