In the first post in this series, we outlined the motivation and theory behind principal component analysis (PCA), which takes points $x_1, \ldots, x_N$ in a high dimensional space to points in a lower dimensional space while preserving as much of the original variance as possible.

In this follow-up post, we apply principal components regression (PCR), an algorithm which includes PCA as a subroutine, to a small dataset to demonstrate the ideas in practice.


To understand this post, you will need to be familiar with the following concepts:

Ordinary least squares

In ordinary least squares (OLS), we want to find a line of best fit between the points $x_1, \ldots, x_N$ and the labels $y_1, \ldots, y_N$.

Denoting by $X$ the matrix whose rows are the points and $y$ the vector whose entries are the labels, the intercept $\alpha$ and slope (a.k.a. gradient) $\beta$ are obtained by minimizing $\Vert \alpha + X \beta - y \Vert$. Some matrix calculus reveals that the minimum is obtained at the values of $\alpha$ and $\beta$ for which \begin{equation} N \alpha = y^\intercal e - \beta^\intercal X^\intercal e \end{equation} and \begin{equation} X^\intercal X \beta = X^\intercal y - \alpha X^\intercal e \end{equation} where $e$ is the vector of all ones.

Principal components regression

The idea behind PCR is simple: instead of doing OLS on the high dimensional space, we first map the points to a lower dimensional space obtained by PCA and then do OLS. In more detail, we

  1. pick a positive integer $k < p$,
  2. construct the matrix $V_k$ whose columns are the first $k$ principal components of $X$,
  3. compute $Z_k = X V_k$, a matrix whose rows are the original points transformed to a lower dimensional “PCA space”, and
  4. perform OLS to find a line of best fit between the transformed points and $y$.

By the previous section, we know that the minimum is obtained at the values of the intercept $\alpha_k$ and gradient $\beta_k$ for which \begin{equation} N \alpha_k = y^\intercal e - \beta_k^\intercal Z_k^\intercal e \end{equation} and \begin{equation} Z_k^\intercal Z_k \beta_k = Z_k^\intercal y - \alpha_k Z_k^\intercal e \end{equation}

Once we have solved these equations for $\alpha_k$ and $\beta_k$, we can predict the label $\hat{y}$ corresponding to a new sample $x$ as $\hat{y} = \alpha_k + x^\intercal V_k \beta_k$.

Computational considerations

Due to the result below, the linear system involving $\alpha_k$ and $\beta_k$ is a (permuted) arrowhead matrix. As such, the system can be solved efficiently.

Lemma. $Z_k^\intercal Z_k = \Sigma_k^2$ where $\Sigma_k$ is the $k \times k$ diagonal matrix whose entries are the first $k$ principal components of $X$ in descending order.

Proof. Let $v_j$ denote the $j$-th column of $V_k$. Since $v_j$ is a principal component of $X$, it is also an eigenvector of $X^\intercal X$ with eigenvalue $\sigma_j^2$, the square of the $j$-th singular value. Therefore, the $(i, j)$-th entry of $Z_k^\intercal Z_k$ is \begin{equation} (X v_i)^\intercal (X v_j) = v_i^\intercal X^\intercal X v_j = \sigma_j^2 v_i^\intercal v_j = \sigma_j^2 \delta_{ij}. \end{equation} where $\delta_{ij}$ is the Kronecker delta.

The California housing dataset

The California housing dataset from [1] has 20,640 samples and 8 predictors. For each $k \leq p = 8$, we fit using PCR on 80% of the samples (the training set) and report the root mean squared error (RMSE) on both the training set and the remaining 20% samples (the test set).

from sklearn.datasets import fetch_california_housing
from sklearn.decomposition import PCA
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
import numpy as np

X_all, y_all = fetch_california_housing(return_X_y=True)
X_train, X_test, y_train, y_test = train_test_split(X_all, y_all, test_size=0.2)

def fit(degree: int = 1):
    """Fits a PCR model to the training data.

    degree: Generates all polynomial combinations of the features up to the specified degree.

    An `sklearn.pipeline.Pipeline` object corresponding to a learned model.
    pipeline = make_pipeline(
        PolynomialFeatures(degree, include_bias=False),
    ), y_train)
    return pipeline

def rmses(model, X, y):
    """Calculates the RMSE of the model using only the first k PCs for each possible value of k.

    model: The `sklearn.pipeline.Pipeline` object returned by the `fit` method.
    X: An (N, p) matrix of features.
    y: The corresponding (N) vector of targets.

    A list of RMSEs (the first entry corresponds to k=1 and the last to k=p).
    pca = model[:-1]
    Z = pca.transform(X)
    regression = model[-1]
    results = []
    for k in range(1, Z.shape[1] + 1):
        pred = regression.intercept_ + Z[:, :k] @ regression.coef_[:k]
        rmse = mean_squared_error(y, pred, squared=False)
    return results

model = fit()
rmses_train = rmses(model, X_train, y_train)
rmses_test = rmses(model, X_test, y_test)
Rank (k) Training set RMSE (in $100,000s) Test set RMSE (in $100,000s)
1 1.15492 1.14608
2 1.14128 1.12491
3 1.13947 1.12457
4 0.853976 0.844071
5 0.853972 0.844123
6 0.810381 0.807466
7 0.729905 0.731179
8 0.723369 0.72742

Both training and test set RMSEs are (roughly) decreasing functions of the rank. This suggests that using all 8 predictors does not cause overfitting.

Deriving predictors

One way to reduce the test set RMSE is to introduce more predictors into the model. Consider, as a toy example, a dataset where each sample $x_i$ has only three predictors: $x_i \equiv (a_i, b_i, c_i)$. We can replace each sample $x_i$ by a new sample $x_i^\prime \equiv (a_i, b_i, c_i, a_i^2, a_i b_i, a_i c_i, b_i^2, b_i c_i, c_i^2)$. In particular, we have added all possible quadratic monomials in $a_i, b_i, c_i$. These new entries are referred to as “derived” predictors. Note that derived predictors need not be quadratic, or even monomials; any function of the original predictors is referred to as a derived predictor.

Returning to the dataset, we add all cubic monomials. It is reasonable to expect that unlike OLS applied to $X$, OLS applied to the derived matrix $X^\prime$ will almost certainly overfit. We plot the results (obtained by setting degree=3 in the previous fitting code) of PCR below, observing the effects of overfitting at large values of the rank.

model = fit(degree=3)
rmses_train = rmses(model, X_train, y_train)
rmses_test = rmses(model, X_test, y_test)



[1] Pace, R. Kelley, and Ronald Barry. “Sparse spatial autoregressions.” Statistics & Probability Letters 33.3 (1997): 291-297.