Principal Component Analysis and Factor Extraction › Mathematical Foundation of PCA

PCA turns a set of $p$ correlated variables into $p$ uncorrelated principal components, ordered by the amount of variance they explain. The step is linear and orthogonal, meaning it preserves distances and angles while rotating the coordinate system to align with directions of maximum variance. This section develops the mathematical framework that makes this step possible, starting from the raw data and building toward the eigenvalue problem that lies at the heart of PCA.

Suppose we have $n$ observations of $p$ variables, organized in an $n \times p$ data matrix $\mathbf{X}$. Each row represents one observation (e.g., one day's yield curve), and each column represents one variable (e.g., the yield at a specific maturity). This matrix structure provides a natural way to organize our data: reading across a row shows us the complete cross-section of variables at one point in time, while reading down a column reveals the time series of a single variable.

The first step is to center the data by subtracting the column means. Centering is essential because PCA seeks directions of maximum variance, and variance is measured relative to the mean. Without centering, the first principal component would simply point toward the overall mean of the data rather than capturing the directions of greatest spread. Let $\bar{\mathbf{x}}$ denote the $p$-dimensional vector of column means. The centered data matrix is:

$$\tilde{\mathbf{X}} = \mathbf{X} - \mathbf{1}_n \bar{\mathbf{x}}^T$$

where:

• $\tilde{\mathbf{X}}$: centered data matrix
• $\mathbf{X}$: original data matrix
• $\mathbf{1}_n$: vector of ones of size $n$
• $\bar{\mathbf{x}}$: vector of column means

In this expression, the outer product $\mathbf{1}_n \bar{\mathbf{x}}^T$ creates an $n \times p$ matrix where every row is identical to the transpose of the mean vector. Subtracting this from the original data ensures that each column of the resulting matrix has mean zero.

The sample covariance matrix is then:

$$\mathbf{\Sigma} = \frac{1}{n-1} \tilde{\mathbf{X}}^T \tilde{\mathbf{X}}$$

where:

• $\mathbf{\Sigma}$: sample covariance matrix
• $n$: number of observations
• $\tilde{\mathbf{X}}$: centered data matrix

To understand why this formula produces the covariance matrix, consider what happens when we multiply $\tilde{\mathbf{X}}^T$ by $\tilde{\mathbf{X}}$. The result is a $p \times p$ matrix where the $(i,j)$ entry equals the sum of products of the centered values for variables $i$ and $j$ across all observations. Dividing by $n-1$ (rather than $n$) gives us the unbiased sample covariance, correcting for the fact that we have estimated the means from the same data. This $p \times p$ symmetric positive semi-definite matrix captures how the variables co-move. The diagonal elements are variances; the off-diagonal elements are covariances.

As we discussed in Part I's treatment of linear algebra, any symmetric matrix can be decomposed into its eigenvalues and eigenvectors. This spectral theorem is the mathematical foundation of PCA, and gives the key insight that allows us to find orthogonal directions of maximum variance. For the covariance matrix:

$$\mathbf{\Sigma} = \mathbf{V} \mathbf{\Lambda} \mathbf{V}^T$$

where:

• $\mathbf{\Sigma}$: covariance matrix
• $\mathbf{\Lambda}$: diagonal matrix of eigenvalues $\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p \geq 0$
• $\mathbf{V}$: orthogonal matrix of eigenvectors
• $\mathbf{v}_k$: $k$-th eigenvector (column of $\mathbf{V}$)
• $p$: total number of variables

The eigenvalues are ordered from largest to smallest, which reflects the need for each direction for explaining variance. Because the covariance matrix is positive semi-definite (a property inherited from its construction as $\tilde{\mathbf{X}}^T \tilde{\mathbf{X}}$ scaled by a positive constant), all eigenvalues are non-negative. The eigenvectors form an orthonormal basis, meaning they are mutually perpendicular and each has unit length. This orthogonality matters because it ensures that the principal components are uncorrelated.

The $k$-th principal component is defined as the projection of the centered data onto the $k$-th eigenvector:

$$\mathbf{z}_k = \tilde{\mathbf{X}} \mathbf{v}_k$$

where:

• $\mathbf{z}_k$: vector of scores for the $k$-th principal component
• $\tilde{\mathbf{X}}$: centered data matrix
• $\mathbf{v}_k$: $k$-th eigenvector (loading vector)

This projection operation computes a weighted sum of the original variables for each observation, where the weights come from the eigenvector. The result is an $n$-dimensional vector of principal component scores, with one score for each observation. The eigenvector $\mathbf{v}_k$ is called the $k$-th loading vector because it specifies how each original variable contributes to the component. A large positive loading means that variable contributes strongly and positively to the component; a large negative loading means it contributes strongly but in the opposite direction.

The first principal component is the linear combination of the original variables with maximum variance. This is more than a convenient property; it is the defining characteristic of PCA. To see why finding this direction of maximum variance corresponds to finding the first eigenvector, consider finding the unit vector $\mathbf{w}$ that maximizes the variance of $\tilde{\mathbf{X}} \mathbf{w}$.

The variance of the projection is:

$$\begin{aligned} \text{Var}(\tilde{\mathbf{X}} \mathbf{w}) &= \frac{1}{n-1} (\tilde{\mathbf{X}} \mathbf{w})^T (\tilde{\mathbf{X}} \mathbf{w}) && \text{(sample variance definition)} \\ &= \frac{1}{n-1} (\mathbf{w}^T \tilde{\mathbf{X}}^T) (\tilde{\mathbf{X}} \mathbf{w}) && \text{(transpose property)} \\ &= \mathbf{w}^T \left( \frac{1}{n-1} \tilde{\mathbf{X}}^T \tilde{\mathbf{X}} \right) \mathbf{w} && \text{(rearrange terms)} \\ &= \mathbf{w}^T \mathbf{\Sigma} \mathbf{w} && \text{(substitute covariance matrix)} \end{aligned}$$

where:

• $\mathbf{w}$: projection vector (weights)
• $\mathbf{\Sigma}$: covariance matrix
• $\tilde{\mathbf{X}}$: centered data matrix
• $n$: number of observations

This derivation reveals a key connection: the variance of any linear combination of our variables equals a quadratic form in the covariance matrix. The quadratic form $\mathbf{w}^T \mathbf{\Sigma} \mathbf{w}$ measures how the covariance structure "stretches" the unit vector $\mathbf{w}$. Our goal is to find the direction that experiences maximum stretching.

We want to maximize this subject to $\mathbf{w}^T \mathbf{w} = 1$. The constraint is necessary because without it, we could increase variance indefinitely by scaling $\mathbf{w}$. By restricting ourselves to unit vectors, we focus on finding the direction of maximum variance rather than an arbitrarily long vector. Using Lagrange multipliers (recall our optimization discussion from Part I), we set:

$$\mathcal{L} = \mathbf{w}^T \mathbf{\Sigma} \mathbf{w} - \lambda (\mathbf{w}^T \mathbf{w} - 1)$$

where:

• $\mathcal{L}$: Lagrangian function
• $\lambda$: Lagrange multiplier (eigenvalue)
• $\mathbf{w}$: projection vector
• $\mathbf{\Sigma}$: covariance matrix

The Lagrangian incorporates our constraint through the multiplier $\lambda$, which will turn out to have a deep interpretation as the eigenvalue. Taking the derivative with respect to $\mathbf{w}$ and setting it to zero:

$$\begin{aligned} \frac{\partial \mathcal{L}}{\partial \mathbf{w}} &= \frac{\partial}{\partial \mathbf{w}} (\mathbf{w}^T \mathbf{\Sigma} \mathbf{w}) - \frac{\partial}{\partial \mathbf{w}} [\lambda (\mathbf{w}^T \mathbf{w} - 1)] && \text{(differentiate terms)} \\ &= 2 \mathbf{\Sigma} \mathbf{w} - 2\lambda \mathbf{w} && \text{(compute gradients)} \\ &= 0 && \text{(first-order condition)} \end{aligned}$$

where:

• $\partial \mathcal{L}/\partial \mathbf{w}$: gradient of the Lagrangian with respect to $\mathbf{w}$
• $\mathbf{\Sigma}$: covariance matrix
• $\mathbf{w}$: projection vector
• $\lambda$: Lagrange multiplier (eigenvalue)

The gradient of the quadratic form $\mathbf{w}^T \mathbf{\Sigma} \mathbf{w}$ equals $2\mathbf{\Sigma} \mathbf{w}$ because of the symmetry of $\mathbf{\Sigma}$, and the gradient of $\mathbf{w}^T \mathbf{w}$ equals $2\mathbf{w}$. This yields the eigenvalue equation:

$$\mathbf{\Sigma} \mathbf{w} = \lambda \mathbf{w}$$

where:

• $\mathbf{\Sigma}$: covariance matrix
• $\mathbf{w}$: eigenvector
• $\lambda$: eigenvalue

This equation has a geometric interpretation. The covariance matrix $\mathbf{\Sigma}$ acts as a linear transformation that stretches and rotates vectors. An eigenvector is special because when $\mathbf{\Sigma}$ acts on it, the result points in the same direction as the original vector, merely scaled by the eigenvalue $\lambda$. These are precisely the directions along which the covariance structure exhibits pure stretching without rotation.

Substituting this into the variance definition:

$$\begin{aligned} \text{Var}(\tilde{\mathbf{X}} \mathbf{w}) &= \mathbf{w}^T \mathbf{\Sigma} \mathbf{w} && \text{(variance definition)} \\ &= \mathbf{w}^T (\lambda \mathbf{w}) && \text{(substitute eigenvalue equation)} \\ &= \lambda (\mathbf{w}^T \mathbf{w}) && \text{(factor out scalar)} \\ &= \lambda && \text{(constraint: $\mathbf{w}^T \mathbf{w} = 1$)} \end{aligned}$$

where:

• $\text{Var}(\tilde{\mathbf{X}} \mathbf{w})$: variance of the projected data
• $\mathbf{w}$: projection vector
• $\mathbf{\Sigma}$: covariance matrix
• $\lambda$: eigenvalue

This result shows that the variance of the projection onto an eigenvector exactly equals the corresponding eigenvalue. To maximize variance, we choose the eigenvector with the largest eigenvalue. This eigenvector points in the direction of greatest spread in our data.

For subsequent components, we impose the additional constraint of orthogonality to previously selected directions. This sequential optimization produces the full set of eigenvectors, ordered by decreasing eigenvalue. The second principal component is the direction of maximum variance among all directions perpendicular to the first, and so on. Because the eigenvectors of a symmetric matrix are automatically orthogonal, the eigendecomposition provides exactly what we need.

The eigenvalues have a direct interpretation: $\lambda_k$ equals the variance of the $k$-th principal component. This is not a coincidence but follows directly from the derivation above. When we project our data onto the $k$-th eigenvector, the resulting variance equals $\lambda_k$. The total variance in the data equals the sum of all eigenvalues (which also equals the trace of the covariance matrix):

$$\text{Total variance} = \sum_{k=1}^{p} \lambda_k = \text{tr}(\mathbf{\Sigma})$$

where:

• $\lambda_k$: $k$-th eigenvalue (variance of $k$-th component)
• $\text{tr}(\mathbf{\Sigma})$: trace of the covariance matrix
• $p$: total number of variables

This identity reflects a fundamental property of the eigendecomposition: the trace of a matrix (the sum of its diagonal elements) equals the sum of its eigenvalues. Since the diagonal elements of the covariance matrix are the variances of individual variables, the total variance is conserved when we rotate to the principal component basis. PCA does not create or destroy variance; it merely redistributes it across orthogonal directions.

The proportion of variance explained by the first $m$ components is:

$$R^2_m = \frac{\sum_{k=1}^{m} \lambda_k}{\sum_{k=1}^{p} \lambda_k}$$

where:

• $R^2_m$: cumulative proportion of variance explained
• $m$: number of retained components
• $\lambda_k$: variance of the $k$-th component
• $p$: total number of variables

This metric guides the choice of how many components to retain. In financial applications, it's common to retain enough components to explain 90-95% of variance, though the appropriate threshold depends on the application. For yield curve modeling, where the first three components typically explain over 99% of variance, the choice is clear. For equity returns, where variance is more dispersed, the decision requires more judgment about the tradeoff between parsimony and fidelity.