A Deep Dive into Principal Component Analysis with PyTorch

Principal component analysis, or PCA, is one of the most widely used techniques for dimensionality reduction in machine learning. It allows us to transform data into a lower dimensional space while preserving the most important information. In this article, we‘ll take an in-depth look at how PCA works, derive it from a mathematical perspective, and implement it using PyTorch. By the end, you‘ll have a solid understanding of this core machine learning algorithm.

The Curse of Dimensionality

In machine learning, we often work with data that has a very large number of features or dimensions. While more features can provide additional information, they also come with many challenges as the dimensionality grows:

• Increased computational and memory requirements for processing the data
• Algorithmic instability and overfitting due to the "curse of dimensionality"
• Difficulty visualizing the data
• Redundant and noisy features

The "curse of dimensionality" refers to various phenomena that arise when working with data in high-dimensional spaces. As the number of features grows, the amount of data needed to generalize accurately grows exponentially. With a fixed number of training samples, the predictive power reduces as dimensionality increases.

Additionally, distance-based algorithms struggle in high dimensions as the distance between any two points becomes less informative. To illustrate, let‘s look at the formulas for Euclidean distance in 1D, 2D, and the general n-dimensional case:

1D: $d(p,q) = |p_x – q_x|$

2D: $d(p,q) = \sqrt{(p_x – q_x)^2 + (p_y – q_y)^2}$

nD: $d(p,q) = \sqrt{\sum_{i=1}^n (p_i – q_i)^2}$

As the number of dimensions $n$ increases, the $n$ differences are summed, making the distance less sensitive to any individual feature. Dimensionality reduction aims to mitigate these issues by projecting the data into a lower dimensional subspace that preserves the most important structure in the data.

Principal Component Analysis

PCA is an unsupervised learning algorithm that finds a new orthogonal coordinate system for the data where the axes (called principal components) align with the directions of maximal variance. Concretely, the first principal component is the direction that maximizes the variance of the projected data. The second principal component is orthogonal to the first and maximizes variance in that direction. This process can be repeated to find any desired number of principal components.

By projecting the data onto the top $k$ principal components, we can reduce the dimensionality while preserving the maximal amount of variance in the data. This projection is optimal in the sense that it minimizes the mean squared reconstruction error between the original data and its reduced dimensional representation.

PCA Step-by-Step

Now let‘s derive PCA mathematically and understand each step of the process. We‘ll assume we have a data matrix $X \in \mathbb{R}^{n \times p}$ consisting of $n$ samples and $p$ features.

Step 1: Standardize the data

First we need to standardize the data so that each feature has zero mean and unit variance. This ensures that the principal components are not influenced by the arbitrary scales of the features. To standardize, we subtract the mean and divide by the standard deviation for each feature:

$X{ij}^{standardized} = \frac{X{ij} – \mu_j}{\sigma_j}$

where $\mu_j$ and $\sigma_j$ are the mean and standard deviation of feature $j$ respectively.

Step 2: Compute the covariance matrix

The covariance matrix $\Sigma \in \mathbb{R}^{p \times p}$ captures the pairwise covariances between all features. The covariance between features $i$ and $j$ is defined as:

$\Sigma{ij} = \frac{1}{n-1} \sum{k=1}^n (X_{ki} – \mui)(X{kj} – \mu_j)$

where $\mu_i$ and $\muj$ are the means of features $i$ and $j$. The diagonal elements $\Sigma{ii}$ are just the variances of each feature. In matrix notation, the covariance matrix can be computed as:

$\Sigma = \frac{1}{n-1} X^TX$

assuming $X$ is already mean-centered (i.e. the column means have been subtracted).

Step 3: Eigendecomposition of covariance matrix

Next we perform eigendecomposition on the covariance matrix to find its eigenvectors and eigenvalues:

$\Sigma = U \Lambda U^T$

where $U \in \mathbb{R}^{p \times p}$ is an orthogonal matrix of eigenvectors and $\Lambda \in \mathbb{R}^{p \times p}$ is a diagonal matrix of corresponding eigenvalues. The eigenvectors represent the principal component directions and the eigenvalues represent the amount of variance explained by each principal component.

Step 4: Select top eigenvectors

We select the top $k$ eigenvectors with the largest eigenvalues to form a projection matrix $W \in \mathbb{R}^{p \times k}$. The columns of $W$ are the top $k$ eigenvectors.

Step 5: Project data onto new subspace

Finally, we can project the data onto the new subspace defined by the top eigenvectors:

$Z = XW$

where $Z \in \mathbb{R}^{n \times k}$ is the reduced dimensionality representation of the data. Each column of $Z$ corresponds to a principal component score.

PCA in PyTorch

Now that we understand the mathematical derivation of PCA, let‘s see how to implement it in PyTorch. PyTorch is a popular deep learning framework that provides a NumPy-like API with GPU acceleration and automatic differentiation. Here‘s a function that performs PCA on a PyTorch tensor:

def pca(X, k):
    # Standardize the data
    X_mean = torch.mean(X, dim=0)
    X = X - X_mean.expand_as(X)

    # Compute covariance matrix
    cov_mat = torch.matmul(X.T, X) / (X.shape[0] - 1)

    # Eigendecomposition of covariance matrix
    eigvals, eigvecs = torch.linalg.eigh(cov_mat)

    # Select top k eigenvectors
    eigvecs = eigvecs[:, -k:]

    # Project data onto new subspace
    Z = torch.matmul(X, eigvecs)

    return Z

To use this function, simply pass in a PyTorch tensor X and the desired number of principal components k. The function returns the projected data Z.

Let‘s compare this to the scikit-learn implementation on a sample dataset:

from sklearn.decomposition import PCA
from sklearn.datasets import load_iris

# Load iris dataset  
iris = load_iris()
X = iris.data

# Perform PCA with scikit-learn
pca_sk = PCA(n_components=2)  
Z_sk = pca_sk.fit_transform(X)

# Perform PCA with PyTorch
X_torch = torch.from_numpy(X).float()
Z_torch = pca(X_torch, k=2).numpy()

# Compare results
print(np.allclose(Z_sk, Z_torch))  # True

As you can see, our PyTorch implementation produces identical results to scikit-learn‘s PCA. However, the PyTorch version can be easily run on a GPU for acceleration on larger datasets:

X_torch = X_torch.to(‘cuda‘)  
Z_torch = pca(X_torch, k=2).cpu().numpy() 

Choosing the Number of Components

One important hyperparameter in PCA is the number of principal components $k$ to retain. There are a few common approaches:

1. Specify the desired amount of variance to preserve (e.g. 95%) and select $k$ to exceed this threshold.
2. Use the "elbow method" by plotting the cumulative explained variance ratio against $k$ and looking for an elbow point where adding more components has diminishing returns.
3. Analyze the scree plot of eigenvalue magnitudes and look for a drop-off.
4. Use domain knowledge to determine how many components are interpretable.

In practice, it‘s often helpful to try multiple values of $k$ and examine the results.

Kernel PCA

One limitation of standard PCA is that it can only find linear subspaces. If the data lies on a nonlinear manifold, PCA may not be able to find a good low-dimensional representation. Kernel PCA is an extension that allows for nonlinear dimensionality reduction.

The idea is to first apply a nonlinear kernel function $\phi$ to map the data into a higher dimensional space. We then perform standard PCA in this high-dimensional feature space. Mathematically, kernel PCA is equivalent to standard PCA using the kernel matrix $K$ instead of the covariance matrix, where:

$K_{ij} = \phi(x_i)^T \phi(x_j)$

Common kernel functions include the polynomial kernel and the radial basis function (RBF) kernel. The advantage is that the kernel matrix can be computed without explicitly evaluating the high-dimensional $\phi$, allowing for efficient computation.

PyTorch does not have a built-in kernel PCA implementation, but it can be easily coded using the kernel trick. See this kernel PCA tutorial for a step-by-step walkthrough.

Applications of PCA

PCA has numerous applications across different fields. Some common use cases include:

• Dimensionality reduction for visualization (plotting data in 2D or 3D)
• Feature extraction and denoising
• Preprocessing for supervised learning algorithms
• Compressing data while preserving important information
• Identifying latent factors or themes in data

For example, PCA is often used in image processing to extract visual features and compress images. It can also be used to visualize and explore high-dimensional datasets such as single-cell RNA-sequencing data.

Conclusion

In this article, we took a deep dive into principal component analysis and its implementation in PyTorch. We started by discussing the curse of dimensionality and the need for dimensionality reduction. We then derived PCA mathematically and walked through each step of the process.

Next we implemented PCA in PyTorch and compared it to the scikit-learn version. We discussed how to choose the number of principal components and extend PCA to nonlinear subspaces with kernel PCA.

I hope this has given you a solid understanding of PCA and how to apply it in practice. PCA is a powerful tool in the machine learning toolbox and PyTorch makes it easy to implement even on large datasets. Try applying PCA to a dataset of your own and see what insights you can uncover!