Principal component analysis

Table of contents

Resources

• Principal component analysis by Hervé Abdi and Lynne J. Williams is excellent at explaining PCA interpretation. It also covers some extensions to PCA.
• A Tutorial on Principal Component Analysis by Jonathon Shlens goes into more detail on the intuition behind PCA, while also discussing its applicability and limits.
• I learnt PCA from these lecture notes from Xavier Gendre. He provides a comprehensive walkthrough using a dataset of Skyrim bows.

Data

PCA assumes you have a dataframe consisting of numerical variables. This includes booleans and integers.

As an example, let’s use a dataset describing the energy mix of each country in the world, in 2019. The dataset is normalized into percentages, which makes countries comparable with one another.

import prince

dataset = prince.datasets.load_energy_mix(year=2019, normalize=True)
dataset.head().style.format('{:.0%}')  

Fitting

The PCA estimator implements the fit/transform API from scikit-learn.

pca = prince.PCA(
    n_components=3,
    n_iter=3,
    rescale_with_mean=True,
    rescale_with_std=True,
    copy=True,
    check_input=True,
    engine='sklearn',
    random_state=42
)
pca = pca.fit(dataset)

The available parameters are:

• n_components — the number of components that are computed. You only need two if your intention is to visualize the two major components.
• n_iter — the number of iterations used for computing the SVD.
• rescale_with_mean — whether to substract each column’s mean
• rescale_with_std — whether to divide each column by it’s standard deviation
• copy — if False then the computations will be done inplace which can have possible side-effects on the input data
• engine — what SVD engine to use (should be one of ['fbpca', 'sklearn'])
• random_state — controls the randomness of the SVD results.

These methods are common to other Prince estimators.

Eigenvalues

The importance of a principal component is indicated by the proportion of dataset inertia it explains. This is called the explained inertia, and is obtained by dividing the eigenvalues obtained with SVD by the total inertia.

The ideal situation is when a large share of inertia is explained by a low number of principal components.

pca.eigenvalues_summary

In this dataset, the first three components explain ~62% of the dataset inertia. This isn’t great, but isn’t bad either.

The columns in the above table are also available as separate attributes:

pca.eigenvalues_

array([1.96301214, 1.56069986, 1.40327823])

pca.percentage_of_variance_

array([24.5376518 , 19.50874822, 17.54097785])

pca.cumulative_percentage_of_variance_

array([24.5376518 , 44.04640003, 61.58737788])

Eigenvalues can also be visualized with a scree plot.

pca.scree_plot()

Coordinates

The row principal coordinates can be obtained once the PCA has been fitted to the data.

pca.transform(dataset).head()

The transform method is in fact an alias for row_coordinates:

pca.transform(dataset).equals(pca.row_coordinates(dataset))

True

This is transforming the original dataset into factor scores.

The column coordinates are obtained during the PCA training process.

pca.column_coordinates_

Visualization

The row and column coordinates be visualized together with a scatter chart.

pca.plot(
    dataset,
    x_component=0,
    y_component=1,
    color_rows_by='continent',
    show_row_markers=True,
    show_column_markers=True,
    show_row_labels=False,
    row_labels_column=None,  # for DataFrames with a MultiIndex
    show_column_labels=False
)

Contributions

pca.row_contributions_.head().style.format('{:.0%}')  

Observations with high contributions and different signs can be opposed to help interpret the component, because these observations represent the two endpoints of this component.

Column contributions are also available.

pca.column_contributions_.style.format('{:.0%}')

Cosine similarities

pca.row_cosine_similarities(dataset).head()

pca.column_cosine_similarities_

Column correlations

pca.column_correlations

(pca.column_correlations ** 2).equals(pca.column_cosine_similarities_)

True

Inverse transformation

You can transform row projections back into their original space by using the inverse_transform method.

reconstructed = pca.inverse_transform(pca.transform(dataset))
reconstructed.head()

This allows measuring the reconstruction error of the PCA. This is usually defined as the $L^2$ norm between the reconstructed dataset and the original dataset.

import numpy as np

np.linalg.norm(reconstructed.values - dataset.values, ord=2)

1.1927159589872411

Supplementary data

Active rows and columns make up the dataset you fit the PCA with. Anything you provide afterwards is considered supplementary data.

For example, we can fit the PCA on all countries that are not part of North America.

active = dataset.query("continent != 'North America'")
pca = prince.PCA().fit(active)

The data for North America can still be projected onto the principal components.

pca.transform(dataset).query("continent == 'North America'")

As for supplementary, they must be provided during the fit call.

pca = prince.PCA().fit(dataset, supplementary_columns=['hydro', 'wind', 'solar'])
pca.column_correlations

There can be supplementary rows and columns at the same time.

pca = prince.PCA().fit(active, supplementary_columns=['hydro', 'wind', 'solar'])

Performance

Under the hood, Prince uses a randomised version of SVD. This is much faster than traditional SVD. However, the results may have a small inherent randomness. This can be controlled with the random_state parameter. The accurary of the results can be increased by providing a higher n_iter parameter.

By default prince uses scikit-learn’s randomized SVD implementation – the one used in TruncatedSVD.

Prince supports different SVD backends. For the while, the only other supported randomized backend is Facebook’s randomized SVD implementation, called fbpca:

prince.PCA(engine='fbpca')

PCA(engine='fbpca')

PCA(engine='fbpca')

You can also use a non-randomized SVD implementation, using the scipy.linalg.svd function:

prince.PCA(engine='scipy')

PCA(engine='scipy')

PCA(engine='scipy')

Here is a very crude benchmark that compares each engine on different dataset sizes.

import itertools
import numpy as np
import pandas as pd

N = [100, 10_000, 100_000]
P = [3, 10, 30]
ENGINES = ['sklearn', 'fbpca', 'scipy']

perf = []

for n, p, engine in itertools.product(N, P, ENGINES):

    # Too slow
    if engine == 'scipy' and n > 10_000:
        continue

    X = pd.DataFrame(np.random.random(size=(n, p)))
    duration = %timeit -q -n 1 -r 3 -o prince.PCA(engine=engine, n_iter=3).fit(X)
    perf.append({
        'n': n,
        'p': p,
        'engine': engine,
        'seconds': duration.average
    })

(
    pd.DataFrame(perf)
    .set_index(['n', 'p'])
    .groupby(['n', 'p'], group_keys=False)
    .apply(lambda x: x[['engine', 'seconds']]
    .sort_values('seconds'))
)