Performing Principal Component Analysis (PCA) with scipy.linalg

Principal Component Analysis (PCA) boils down to finding the directions in your data along which variance is maximized. These directions, called principal components, are orthogonal to each other and capture the essence of the data’s structure by reducing dimensionality without losing too much information.

Mathematically, suppose you have a data matrix X where each row is an observation and each column is a feature. The first step is usually to center X by subtracting the mean of each column, ensuring the data is zero-centered. This step very important because PCA is sensitive to the scaling and centering of data.

Once centered, PCA looks for a set of vectors v such that the variance of the projections Xv is maximized. This can be expressed as the optimization problem:

maximize   vT Σ v   subject to   vT v = 1

Here, Σ is the covariance matrix of the centered data. The constraint forces v to be a unit vector, preventing the variance from inflating infinitely.

The solution to this problem comes from the eigen-decomposition of Σ. The eigenvectors of Σ represent the directions of maximum variance, and the corresponding eigenvalues quantify the amount of variance captured along those directions. Sorting eigenvalues in descending order gives the principal components in order of importance.

Another way to see PCA is through Singular Value Decomposition (SVD). If you decompose your zero-centered data matrix X as X = U Σ VT, the columns of V are your principal directions, and the singular values in Σ relate to the variance explained by each component. This approach can be more numerically stable and efficient, especially for large datasets.

It’s worth noting that PCA is fundamentally a linear transformation that projects data onto a new coordinate system defined by these principal components. The first principal component captures the most variance, the second captures the most remaining variance orthogonal to the first, and so on.

Interpreting PCA also means understanding what the principal components represent in the context of your data, which can sometimes be tricky. The coefficients of the eigenvectors (or loadings) reveal how much each original feature contributes to a component, helping you make sense of the directions in terms of your original variables.

To summarize the math succinctly:

1. Center the data:  X_centered = X - mean(X)
2. Compute covariance matrix:  Σ = (1/(n-1)) * X_centered.T @ X_centered
3. Solve eigenvalue problem:  Σ v = λ v
4. Sort eigenvectors by eigenvalues descending
5. Project data:  Y = X_centered @ V_selected

In practice, the choice of how many components to keep depends on the fraction of variance you want to retain, often visualized by a scree plot or cumulative variance graph.

Implementing PCA with scipy.linalg

Implementing PCA using scipy.linalg involves directly working with the covariance matrix and its eigen-decomposition, giving you full control over each step. This is instructive and useful for understanding the underlying mechanics without relying on higher-level abstractions.

First, let’s import the necessary libraries and generate some synthetic data for demonstration:

import numpy as np
from scipy.linalg import eigh

# Generate synthetic data: 100 samples, 3 features
np.random.seed(0)
X = np.dot(np.random.randn(100, 3), np.array([[0.5, 1.5, 0.3],
                                             [1.0, 0.0, 1.0],
                                             [0.0, 1.0, 0.5]])) + 2

Centering the data is essential before applying PCA:

X_mean = np.mean(X, axis=0)
X_centered = X - X_mean

Next, compute the covariance matrix. Here, we use the unbiased estimator dividing by n - 1:

n_samples = X_centered.shape[0]
cov_matrix = (X_centered.T @ X_centered) / (n_samples - 1)

Since the covariance matrix is symmetric and positive semi-definite, we use scipy.linalg.eigh which is optimized for such matrices and returns sorted eigenvalues and eigenvectors:

eigenvalues, eigenvectors = eigh(cov_matrix)

Note that eigh returns eigenvalues in ascending order, so we reverse them along with eigenvectors to have descending order:

eigenvalues = eigenvalues[::-1]
eigenvectors = eigenvectors[:, ::-1]

Now, you can select the number of principal components to keep. Suppose you want to keep two components:

n_components = 2
selected_vectors = eigenvectors[:, :n_components]

Project the centered data onto these principal components to get the transformed dataset:

X_pca = X_centered @ selected_vectors

To verify the amount of variance explained by each component, compute the explained variance ratio:

explained_variance_ratio = eigenvalues / np.sum(eigenvalues)
print("Explained variance ratio:", explained_variance_ratio[:n_components])

This ratio helps decide how many components are meaningful to retain. For example, if the first two components explain over 90% of the variance, you might consider them sufficient.

Putting it all together, here’s a concise function encapsulating PCA with scipy.linalg:

def pca_scipy(X, n_components):
    X_centered = X - np.mean(X, axis=0)
    cov_matrix = np.cov(X_centered, rowvar=False)
    eigenvalues, eigenvectors = eigh(cov_matrix)
    eigenvalues = eigenvalues[::-1]
    eigenvectors = eigenvectors[:, ::-1]
    components = eigenvectors[:, :n_components]
    X_pca = X_centered @ components
    explained_variance_ratio = eigenvalues / np.sum(eigenvalues)
    return X_pca, explained_variance_ratio[:n_components], components

Using this function, you can quickly perform PCA and retrieve the transformed data, variance ratios, and principal directions:

X_pca, var_ratio, components = pca_scipy(X, 2)
print("Projected data shape:", X_pca.shape)
print("Variance explained:", var_ratio)

One subtlety to watch out for is the scale of your features. PCA is sensitive to the units and variances of the original features, so if they differ significantly, it is common to standardize the data (zero mean and unit variance) before applying PCA.

If you want to standardize, you can use:

X_mean = np.mean(X, axis=0)
X_centered = X - X_mean

Then apply the PCA function on X_scaled rather than X. This ensures each feature contributes equally to the analysis, preventing dominance by features with larger scales.

Finally, while this manual approach provides transparency, for larger datasets or more convenience, consider using sklearn.decomposition.PCA, which leverages optimized routines internally and offers additional functionality like automatic dimensionality selection and inverse transforms.

However, mastering the underlying implementation with scipy.linalg is invaluable for grasping what PCA truly does under the hood, which is essential for debugging, customization, or educational purposes. Next, we will explore how to interpret and visualize these results to make them actionable in real-world scenarios, focusing on loading plots and variance graphs that reveal the structure of your data and the meaning of each principal component…

Interpreting and visualizing PCA results

Once you have computed the principal components, interpreting and visualizing these results becomes crucial for understanding the underlying structure of your data. One effective way to visualize PCA results is through loading plots, which show how much each original feature contributes to the principal components.

To create a loading plot, you can use the eigenvectors, which represent the directions of the principal components. Each component’s eigenvector coefficients indicate the contribution of each feature to that component. A higher absolute value means a stronger contribution. Here’s how to generate a loading plot using matplotlib:

import matplotlib.pyplot as plt

def plot_loadings(components, feature_names):
    plt.figure(figsize=(10, 6))
    for i in range(components.shape[1]):
        plt.scatter(np.arange(len(feature_names)), components[:, i], label=f'PC {i + 1}')
    plt.xticks(np.arange(len(feature_names)), feature_names, rotation=45)
    plt.axhline(0, color='grey', lw=1)
    plt.title('PCA Loadings')
    plt.xlabel('Features')
    plt.ylabel('Loading values')
    plt.legend()
    plt.tight_layout()
    plt.show()

To use this function, you would pass the principal component loadings and the feature names:

feature_names = ['Feature 1', 'Feature 2', 'Feature 3']
plot_loadings(components, feature_names)

Another important visualization is the explained variance ratio plot. This plot allows you to see how much variance each principal component captures, helping you decide how many components to retain. Here’s a simple implementation:

def plot_explained_variance(explained_variance_ratio):
    plt.figure(figsize=(10, 6))
    plt.plot(range(1, len(explained_variance_ratio) + 1), explained_variance_ratio, marker='o', linestyle='--')
    plt.title('Explained Variance by Principal Components')
    plt.xlabel('Principal Component')
    plt.ylabel('Explained Variance Ratio')
    plt.axhline(y=0.95, color='r', linestyle='--')
    plt.xticks(range(1, len(explained_variance_ratio) + 1))
    plt.grid()
    plt.tight_layout()
    plt.show()

You can call this function with the explained variance ratio from your PCA results:

plot_explained_variance(var_ratio)

These visualizations not only enhance your understanding of the PCA results but also provide an intuitive way to communicate findings to stakeholders. They help elucidate how the original features influence the new principal components, making the analysis actionable and relatable.

Additionally, you can visualize the transformed data in a lower-dimensional space. For instance, if you reduce your data to two dimensions, a scatter plot can reveal clustering or separation patterns among different classes or groups. Here’s how you can plot the PCA-transformed data:

def plot_pca_scatter(X_pca, labels=None):
    plt.figure(figsize=(10, 6))
    if labels is not None:
        plt.scatter(X_pca[:, 0], X_pca[:, 1], c=labels, cmap='viridis', edgecolor='k', s=50)
    else:
        plt.scatter(X_pca[:, 0], X_pca[:, 1], edgecolor='k', s=50)
    plt.title('PCA Transformed Data')
    plt.xlabel('Principal Component 1')
    plt.ylabel('Principal Component 2')
    plt.grid()
    plt.tight_layout()
    plt.show()

To visualize your PCA results with labels, you would call:

plot_pca_scatter(X_pca, labels)

In this way, PCA not only simplifies the data but also reveals insights that are often obscured in higher dimensions. By interpreting and visualizing your PCA results effectively, you can uncover the latent structures in your data, enhancing both analysis and decision-making processes.