Lecture 3: Introduction to Principal Component Analysis (PCA)

UBC Master of Data Science program, 2024-25

Imports and learning outcomes

Imports

import os
import sys

import numpy as np
import pandas as pd

sys.path.append(os.path.join(os.path.abspath(".."), "code"))

import matplotlib.pyplot as plt
from plotting_functions import *
from sklearn.decomposition import PCA
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

plt.rcParams["font.size"] = 16
plt.rcParams["figure.figsize"] = (5, 4)
%matplotlib inline
pd.set_option("display.max_colwidth", 0)

import plotly.io as pio
# pio.renderers.default = "png"

%config InlineBackend.figure_formats = ['svg']

plt.rcParams.update({'font.size': 12, 'axes.labelweight': 'bold', 'figure.figsize': (6, 4)})

DATA_DIR = os.path.join(os.path.abspath(".."), "data/")

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 10
      7 sys.path.append(os.path.join(os.path.abspath(".."), "code"))
      9 import matplotlib.pyplot as plt
---> 10 from plotting_functions import *
     11 from sklearn.decomposition import PCA
     12 from sklearn.pipeline import make_pipeline

File ~/MDS/2024-25/563/DSCI_563_unsup-learn/lectures/code/plotting_functions.py:7
      5 import matplotlib.pyplot as plt
      6 from matplotlib.colors import ListedColormap, colorConverter, LinearSegmentedColormap
----> 7 from scipy.spatial import distance
      8 from sklearn.metrics import euclidean_distances
      9 from sklearn.manifold import MDS

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/spatial/__init__.py:110
      1 """
      2 =============================================================
      3 Spatial algorithms and data structures (:mod:`scipy.spatial`)
   (...)
    107    QhullError
    108 """  # noqa: E501
--> 110 from ._kdtree import *
    111 from ._ckdtree import *  # type: ignore[import-not-found]
    112 from ._qhull import *

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/spatial/_kdtree.py:4
      1 # Copyright Anne M. Archibald 2008
      2 # Released under the scipy license
      3 import numpy as np
----> 4 from ._ckdtree import cKDTree, cKDTreeNode  # type: ignore[import-not-found]
      6 __all__ = ['minkowski_distance_p', 'minkowski_distance',
      7            'distance_matrix',
      8            'Rectangle', 'KDTree']
     11 def minkowski_distance_p(x, y, p=2):

File _ckdtree.pyx:11, in init scipy.spatial._ckdtree()

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/sparse/__init__.py:315
    312 from ._sputils import get_index_dtype, safely_cast_index_arrays
    314 # For backward compatibility with v0.19.
--> 315 from . import csgraph
    317 # Deprecated namespaces, to be removed in v2.0.0
    318 from . import (
    319     base, bsr, compressed, construct, coo, csc, csr, data, dia, dok, extract,
    320     lil, sparsetools, sputils
    321 )

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/sparse/csgraph/__init__.py:187
    158 __docformat__ = "restructuredtext en"
    160 __all__ = ['connected_components',
    161            'laplacian',
    162            'shortest_path',
   (...)
    184            'csgraph_to_masked',
    185            'NegativeCycleError']
--> 187 from ._laplacian import laplacian
    188 from ._shortest_path import (
    189     shortest_path, floyd_warshall, dijkstra, bellman_ford, johnson, yen,
    190     NegativeCycleError
    191 )
    192 from ._traversal import (
    193     breadth_first_order, depth_first_order, breadth_first_tree,
    194     depth_first_tree, connected_components
    195 )

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/sparse/csgraph/_laplacian.py:7
      5 import numpy as np
      6 from scipy.sparse import issparse
----> 7 from scipy.sparse.linalg import LinearOperator
      8 from scipy.sparse._sputils import convert_pydata_sparse_to_scipy, is_pydata_spmatrix
     11 ###############################################################################
     12 # Graph laplacian

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/sparse/linalg/__init__.py:129
      1 """
      2 Sparse linear algebra (:mod:`scipy.sparse.linalg`)
      3 ==================================================
   (...)
    126 
    127 """
--> 129 from ._isolve import *
    130 from ._dsolve import *
    131 from ._interface import *

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/sparse/linalg/_isolve/__init__.py:4
      1 "Iterative Solvers for Sparse Linear Systems"
      3 #from info import __doc__
----> 4 from .iterative import *
      5 from .minres import minres
      6 from .lgmres import lgmres

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/sparse/linalg/_isolve/iterative.py:5
      3 from scipy.sparse.linalg._interface import LinearOperator
      4 from .utils import make_system
----> 5 from scipy.linalg import get_lapack_funcs
      7 __all__ = ['bicg', 'bicgstab', 'cg', 'cgs', 'gmres', 'qmr']
     10 def _get_atol_rtol(name, b_norm, atol=0., rtol=1e-5):

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/linalg/__init__.py:203
      1 """
      2 ====================================
      3 Linear algebra (:mod:`scipy.linalg`)
   (...)
    200 
    201 """  # noqa: E501
--> 203 from ._misc import *
    204 from ._cythonized_array_utils import *
    205 from ._basic import *

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/linalg/_misc.py:3
      1 import numpy as np
      2 from numpy.linalg import LinAlgError
----> 3 from .blas import get_blas_funcs
      4 from .lapack import get_lapack_funcs
      6 __all__ = ['LinAlgError', 'LinAlgWarning', 'norm']

File ~/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/linalg/blas.py:213
    210 import numpy as np
    211 import functools
--> 213 from scipy.linalg import _fblas
    214 try:
    215     from scipy.linalg import _cblas

ImportError: dlopen(/Users/kvarada/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/linalg/_fblas.cpython-312-darwin.so, 0x0002): Library not loaded: @rpath/libgfortran.5.dylib
  Referenced from: <0B9C315B-A1DD-3527-88DB-4B90531D343F> /Users/kvarada/miniforge3/envs/jbook/lib/libopenblas.0.dylib
  Reason: tried: '/Users/kvarada/miniforge3/envs/jbook/lib/libgfortran.5.dylib' (duplicate LC_RPATH '@loader_path'), '/Users/kvarada/miniforge3/envs/jbook/lib/libgfortran.5.dylib' (duplicate LC_RPATH '@loader_path'), '/Users/kvarada/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/linalg/../../../../libgfortran.5.dylib' (duplicate LC_RPATH '@loader_path'), '/Users/kvarada/miniforge3/envs/jbook/lib/python3.12/site-packages/scipy/linalg/../../../../libgfortran.5.dylib' (duplicate LC_RPATH '@loader_path'), '/Users/kvarada/miniforge3/envs/jbook/bin/../lib/libgfortran.5.dylib' (duplicate LC_RPATH '@loader_path'), '/Users/kvarada/miniforge3/envs/jbook/bin/../lib/libgfortran.5.dylib' (duplicate LC_RPATH '@loader_path'), '/usr/local/lib/libgfortran.5.dylib' (no such file), '/usr/lib/libgfortran.5.dylib' (no such file, not in dyld cache)

Learning outcomes

From this lecture, students are expected to be able to:

• Explain some issues caused by high-dimensional data and the need for dimensionality reduction.

• Explain the intuition behind Principal Component Analysis (PCA).

• Describe the role and shapes of four matrices \(X\), \(W\), \(Z\), and \(\hat{X}\) in the context of dimensionality reduction techniques;

• Explain how to get \(Z\) from \(X\) and \(W\).

• Explain how to get \(\hat{X}\) from \(Z\) and \(W\).

• State the loss function of PCA.

• Explain the difference between PCA and linear regression.

• Broadly explain how PCA is learned using SVD.

• Explain how PCA can be used in data compression, better representation, and visualization.

• Use sklearn.decomposition.PCA to perform Principal Component Analysis.

• Use sklearn’s inverse_transform to get reconstructions.

Dimensionality reduction: Motivation and introduction [video]

Motivation

• Suppose you’re shown the picture below and you are told that this is Eva.

• Do you have to remember every pixel in the image to recognize other pictures of Eva?

• For example, if you are asked which one is Eva in the following pictures, it’ll be fairly easy for you to identify her just based on some high-level features.

• Just remembering important features such as shape of eyes, nose, mouth, shape and colour of hair etc. suffice to tell her apart from other people.

• Can we learn such high-level features or the summary of the given raw features with machine learning models?

• Yes! With dimensionality reduction techniques!

• As data scientists, given a dataset we either want to understand some phenomenon or build predictive models.

• Very often the data we work with is clouded, complex, unclear, or even redundant.

• But in reality the underlying phenomenon we are trying to understand or the relationship between variables in the data is much simpler.

• Dimensionality reduction is useful in such scenarios.

Toy example: nutritional value of pizzas

• Suppose we want to analyze nutritional value of pizzas of different brands.

• Here is a toy dataset for this problem.

pizza_df = pd.read_csv(DATA_DIR + "Pizza.csv")
pizza_df.head()

	brand	id	mois	prot	fat	ash	sodium	carb	cal
0	A	14069	27.82	21.43	44.87	5.11	1.77	0.77	4.93
1	A	14053	28.49	21.26	43.89	5.34	1.79	1.02	4.84
2	A	14025	28.35	19.99	45.78	5.08	1.63	0.80	4.95
3	A	14016	30.55	20.15	43.13	4.79	1.61	1.38	4.74
4	A	14005	30.49	21.28	41.65	4.82	1.64	1.76	4.67

X_pizza = pizza_df.drop(columns=["id", "brand"])
y_pizza = pizza_df["brand"]
X_pizza.head()

	mois	prot	fat	ash	sodium	carb	cal
0	27.82	21.43	44.87	5.11	1.77	0.77	4.93
1	28.49	21.26	43.89	5.34	1.79	1.02	4.84
2	28.35	19.99	45.78	5.08	1.63	0.80	4.95
3	30.55	20.15	43.13	4.79	1.61	1.38	4.74
4	30.49	21.28	41.65	4.82	1.64	1.76	4.67

X_pizza.shape

(300, 7)

We have features such as amount of moisture, amount of protein, amount of fat, amount of ash, amount of sodium, and amount of carbohydrates, and amount of calories per 100 grams in the dataset.

Let’s examine correlations between different variables.

corr_heatmat(X_pizza.corr(), w=6, h=4)
plt.show();

• There is redundancy in the data; many features are correlated.

• Can we summarize these features in some meaningful way so that the data is cleaner and less redundant?

• Can we just discard some redundant features?

• We have seen some (not very satisfactory) feature selection methods to identify least important features in a greedy way and throw away such features.

• This week we are going to look at a class of more sophisticated approaches for this, which are typically referred to as dimensionality reduction.

What is dimensionality reduction?

Dimensionality reduction is the task of summarizing data or reducing a dataset in high dimension (e.g., 1000) to low dimension (e.g., 10) while retaining the most “important” characteristics of the data.

Dimensionality reduction is also used to reduce the dimensionality similar to feature selection. But

• We will not be just dropping columns as we did in feature selection.

• The idea of (linear) dimensionality reduction is to project high dimensional data to low dimensional space while retaining the most “important” characteristics of the data.

• We can also reconstruct the original data (with some error) from this transformed data.

How do we reduce the dimensions?

• The techniques we are going to look at this week summarize the data by creating new features which are linear combinations of the original features.

• Example: $\(\text{new\_feature} = 0.44 \times fat + 0.47 \times ash - 0.42 \times carb \dots \)$

Dimensionality reduction toy example

• Let’s apply a popular dimensionality reduction technique called Principal Component Analysis (PCA) using sklearn’s PCA on our nutritional value of pizzas toy data.

• Learning a PCA model and transforming data is similar to applying preprocessing transformations in sklearn.

• You can learn a PCA model and transform the data using fit and transform methods, respectively.

n_components = (
    2  # summarize the data with only two features (typically called components)
)
pipe_pca = make_pipeline(
    StandardScaler(), PCA(n_components=n_components)
)  # scaling before PCA is crucial. We'll see the reason later.
Z = pipe_pca.fit_transform(X_pizza)  # transform the data

How does the data look like after dimensionality reduction?

Z_labels = ["Z" + str(i + 1) for i in range(n_components)]
pd.DataFrame(Z, columns=Z_labels, index=X_pizza.index).head()

	Z1	Z2
0	5.010343	-2.679215
1	5.023755	-2.529295
2	4.805439	-2.673700
3	4.469543	-2.285029
4	4.471893	-2.159152

• We have reduced dimensionality from original 7 features to 2 features.

• The two new features can be thought of as the summary of the original features.

• It has learned the “most informative” linear combinations of the features.

• Each new feature (principal component) has a coefficient associated with each of the original features and the value of the new feature is a linear combination of the original features.

W_labels = ["PC" + str(i + 1) for i in range(n_components)]
W = pipe_pca.named_steps["pca"].components_
pd.DataFrame(W, columns=X_pizza.columns, index=W_labels)

	mois	prot	fat	ash	sodium	carb	cal
PC1	0.064709	0.378761	0.446666	0.47189	0.435703	-0.424914	0.244487
PC2	0.628276	0.269707	-0.234379	0.11099	-0.201662	-0.320312	-0.567458

\[\text{Z1} = 0.064709 \times \text{mois} + 0.378761 \times \text{prot} + \dots + -0.424914 \times \text{carb} + 0.244487 \times \text{cal}\]

\[\text{Z2} = -0.628276 \times \text{mois} + -0.628276 \times \text{prot} + \dots + 0.320312 \times \text{carb} + 0.567458 \times \text{cal}\]

np.round(Z[0, :], 4)  # transformed values for the 0th example

array([ 5.0103, -2.6792])

x0_scaled = pipe_pca.named_steps["standardscaler"].transform(X_pizza)[0, :]
np.round((np.dot(x0_scaled, W[0, :]), np.dot(x0_scaled, W[1, :])), 4)

array([ 5.0103, -2.6792])

pca_components = plot_pca_w_vectors(W, W_labels, X_pizza.columns, width=800, height=800)
pca_components.show()