9.2 Introduction to Principal Component Analysis (PCA)#
Principal Component Analysis (PCA) is a technique used to emphasize variation and bring out strong patterns in a dataset. It’s often used to make data easy to explore and visualize or to reduce the number of variables in a dataset while preserving the variation in the dataset as much as possible. It does this by creating new variables called principal components that are linear combinations of the original variables. The first principal component accounts for as much of the variation in the data as possible, and each succeeding component accounts for as much of the remaining variation as possible. PCA is closely related to singular value decomposition (SVD).
PCA in a Nutshell#
PCA is one of the central applications of the SVD, where we can transform high-dimensional correlated data. In PCA we pre-processes the data by mean subtraction and setting the variance to unity before performing the SVD. The geometry of the resulting coordinate system is determined by principal components (PCs) that are uncorrelated (orthogonal) to each other, but have maximal correlation with the measurements.
We again consider a dataset of many (real valued) Measurements stored in a matrix \(X \in \mathbb{R}^{m \times n}\), where \(n\) is the number of measurements and \(m\) is the number of features (please note that we changed the dimensions compared to the SVD section). We can write this as:
where each measurement is \(x_i \in \mathbb{R}^m\) and the whole dataset is \(X \in \mathbb{R}^{m \times n}\).
First we compute the mean of \(X\) for each measurement:
and construct the mean Matrix \(X_M\):
Then we subtract the mean from each measurement resulting in the mean substracted data \(B\):
Using the mean subtracted data \(B\) we compute the covariance matrix of \(B\):
The principal components of \(X\) are also the eigenvectors of \(C_X\). Hence, if we calculate the SVD of \(X\), the columns of matrix \(V\) contain the eigenvectors of \(B^TB = C_X\) (see the description of \(V\) in the SVD section). Therefore, the columns of \(V\) are the principal components of \(X\).
In summary we can write the following pseudo code for PCA:
# PCA - Perform PCA using SVD.
# Input: data - MxN matrix of input data (M dimensions, N trials)
[M,N] = size(data)
mean_x = mean(X, axis=2)
data = data - repmat(mean_x,1,N)
B = data* / sqrt(N-1)
[u,S,PC] = svd(B)
# calculate the variances
S = diag(S)
V = S .* S
# project the original data
signals = PC’ * data
Example: Measuring a spring-mass system using cameras#
We will use an example from “A Tutorial on Principal Component Analyses” by Jonathon Shlens (https://arxiv.org/pdf/1404.1100.pdf) to illustrate the use of PCA.
Consider the following example. We have an ideal spring-mass system of a ball with mass m attached to a massless, frictionless spring. The spring is attached to a wall and the ball is free to move in the horizontal direction. The ball is initially at rest at the equilibrium position of the spring. We then pull the ball to the right and release it. The ball will then oscillate back and forth. We are interested in measuring the position of the mass over time. We have three video cameras which we can use to track the position of the mass (see figure below).
Since we do not know anything about the system, we did not position the cameras in a position which would allow us to directly measure the position of the mass. Instead, we positioned the cameras in arbitrary angles. By doing this we can only measure the position of the mass in the x- and y-directions of each camera frame. Hence we want to use the measured camera positions to calculate the position of the mass in its x-direction. How do we get from this data set to a simple equation of x?
First let start with the problem and the data.
import Pkg
Pkg.instantiate()
import Pkg
# # Pkg.generate("PCA")
Pkg.activate("PCA")
Pkg.add("Revise")
Pkg.add("Plots")
Pkg.add("LinearAlgebra")
Pkg.add("Statistics")
Pkg.add("Random")
Pkg.add("Revise")
Activating project at `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA`
Resolving package versions...
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Project.toml`
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Manifest.toml`
Resolving package versions...
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Project.toml`
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Manifest.toml`
Resolving package versions...
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Project.toml`
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Manifest.toml`
Resolving package versions...
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Project.toml`
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Manifest.toml`
Resolving package versions...
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Project.toml`
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Manifest.toml`
Resolving package versions...
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Project.toml`
No Changes to `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA/Manifest.toml`
import Pkg
Pkg.activate("PCA")
using Revise
include("./SpringMass.jl") # this file implements the generation of our spring-mass toy problem
using .SpringMass
using Statistics, LinearAlgebra, Random, Plots
Activating project at `~/Library/Mobile Documents/com~apple~CloudDocs/Projects/numerical_methods/09-PrincipalComponentAnalysis/PCA`
Let’s generate our dataset and take a first look at it:
amplitude = 2.0 # amplitude of the motion
X, x_true, y_true = generate_dataset(amplitude)
([2.0340980695392608 1.0401178217442442 … 0.5406565701027691 2.1916883467806425; 1.79238689925546 0.8712672617758772 … 0.882720529981148 2.238250204375316; … ; 1.7256267194994948 0.6334691274186726 … 0.37045375973481676 1.8678621037055916; 1.642483629682659 0.9927447197856961 … 0.36319428701059026 1.964475863750956], [2.0, 1.9960534568565431, 1.9842294026289558, 1.9645745014573774, 1.9371663222572622, 1.902113032590307, 1.859552971776503, 1.809654104932039, 1.7526133600877272, 1.6886558510040302 … 1.6886558510040273, 1.7526133600877276, 1.809654104932036, 1.8595529717764974, 1.9021130325903046, 1.9371663222572653, 1.964574501457378, 1.9842294026289553, 1.9960534568565433, 2.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
So what do we have here?
The dataset X is a matrix where each row corresponds to a different time step and each column corresponds to a different measurement. The measurements are the x and y coordinates of the ball’s position as recorded by three different cameras. The cameras are not aligned with the true x and y axes, but instead have their own coordinate systems that are rotated at arbitrary angles with respect to the true axes.
Before we take a look at the data we plot the true motion of the ball we want to measure. This will give us an idea what we expect to see in the data.
# Plot the "true" motion
p = plot(x_true, label="True Motion", title="True Motion of the Ball", xlabel="Time", ylabel="x")
Now let’s take a look at the data. We have three cameras which measure the position of the ball in the x- and y-directions. We can plot the position of the ball in the x- and y-directions for each camera:
# Scatter Plot the 2-dim data
p1 = scatter(X[:, 1], X[:, 2], label="Camera 1")
p2 = scatter(X[:, 3], X[:, 4], label="Camera 2")
p3 = scatter(X[:, 5], X[:, 6], label="Camera 3")
plot(p1, p2, p3, title="Scatter Plot of Data", xlabel="x", ylabel="y")
Here’s a breakdown of the columns in X:
Columns 1 and 2: The x and y coordinates of the ball’s position as recorded by camera 1 (Plot 1).
Columns 3 and 4: The x and y coordinates of the ball’s position as recorded by camera 2 (Plot 2).
Columns 5 and 6: The x and y coordinates of the ball’s position as recorded by camera 3 (Plot 3).
So these 6 columns are the measurements we have. Where each measurement \(x_{i}\) at time step \(i\) can be viewed as:
Each row in X represents a different time step \(i\). The time steps are evenly spaced and cover the entire duration of the ball’s motion. We simulated for 10 seconds with a time step of 0.01 seconds, so there are 1000 rows in X.
The values in X are real numbers that represent the position of the ball in the coordinate system of each camera. These values have been corrupted by Gaussian noise to simulate measurement error.
In order to analyze this dataset, we will perform Principal Component Analysis (PCA) to find the directions of greatest variance in the data. These directions, known as the principal components, can reveal the true motion of the ball, as well as the relationships between the different measurements.
PCA on the spring-mass system#
Now let us compute the principal components of our data set. We will preprocess the data and use the eigen() function od the LinearAlgebra package to compute the principal components as explained above.
# Perform PCA
μ = mean(X, dims=1) # compute the mean
X_centered = X .- μ # center the data
Σ = X_centered' * X_centered / (size(X, 1) - 1) # compute the covariance matrix
λ, V = eigen(Σ) # compute the eigenvalues and eigenvectors
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
6-element Vector{Float64}:
0.057666530374654114
0.0592954555159329
0.060266344778793626
0.06322771740789107
0.06828922302185023
6.100006964226926
vectors:
6×6 Matrix{Float64}:
-0.374712 0.622082 -0.39322 -0.216319 -0.0249224 -0.520161
-0.770163 -0.522829 0.229146 0.0163289 -0.146459 -0.243464
-0.305434 0.581104 0.628052 0.235513 -0.0308123 0.343748
-0.093659 0.0278164 -0.210895 -0.403731 -0.753349 0.464159
0.0831243 0.0279124 -0.219372 0.817322 -0.494869 -0.176851
0.396823 0.0206969 0.55301 -0.257797 -0.405646 -0.552517
# The principal components are given by the eigenvectors
PCs = V
6×6 Matrix{Float64}:
-0.374712 0.622082 -0.39322 -0.216319 -0.0249224 -0.520161
-0.770163 -0.522829 0.229146 0.0163289 -0.146459 -0.243464
-0.305434 0.581104 0.628052 0.235513 -0.0308123 0.343748
-0.093659 0.0278164 -0.210895 -0.403731 -0.753349 0.464159
0.0831243 0.0279124 -0.219372 0.817322 -0.494869 -0.176851
0.396823 0.0206969 0.55301 -0.257797 -0.405646 -0.552517
# Plot the original data and the principal components
p1 = plot(X, label=["Camera 1 x" "Camera 1 y" "Camera 2 x" "Camera 2 y" "Camera 3 x" "Camera 3 y"], title="Original Data")
p2 = plot(X * PCs, label=["PC 1" "PC 2" "PC 3" "PC 4" "PC 5" "PC 6"], title="Principal Components")
plot(p1, p2, layout=(1, 2), size=(1200, 400), legend=:topleft)
Things we can do with PCA#
Correlation: If your data variables are correlated, PCA can help identify these correlations. Highly correlated variables will contribute to the same principal component.
Dimensionality Reduction: If your data has many dimensions, PCA can be used to reduce its dimensionality while preserving as much of the data’s variation as possible. The reduced-dimension data can be easier to analyze and visualize.
Variance Explained: The eigenvalues of the covariance matrix in PCA represent the amount of variance explained by each principal component. The first principal component explains the most variance, the second principal component explains the second most, and so on. By examining these eigenvalues, you can get a sense of how much of the total variance in your data is captured by each principal component.
Outlier Detection: Outliers are observations that are far from the main cluster of data points. PCA can help detect outliers by looking at the principal components that explain the most variance. Outliers will often show up as extreme values on the principal components.
# combine the transformed data into a matrix
X_projected = X * PCs
1001×6 Matrix{Float64}:
-0.0487902 -0.153135 -0.148671 -0.356575 -0.279929 -3.78442
0.211935 -0.18175 -0.100737 -0.00251267 -0.388454 -3.72006
-0.23903 0.0926737 -0.116734 -0.216991 -0.233949 -2.9678
0.0285837 0.107767 0.323038 0.0465011 -0.0687001 -3.526
0.213082 -0.129324 0.100039 -0.0260301 -0.254558 -3.62457
0.328612 0.229011 -0.0812484 -0.00570254 -0.0334571 -3.30068
0.202881 -0.235805 -0.0870228 -0.216485 -0.0251338 -3.40368
-0.160716 0.175451 -0.714523 0.0737911 -0.0931874 -3.17933
0.174087 -0.202091 -0.208311 0.332084 -0.162269 -3.12894
-0.0367802 -0.448312 0.0196243 -0.104477 0.0478867 -2.99175
0.262592 -0.252778 0.546094 0.141217 0.0398668 -3.25365
0.273004 0.135394 -0.269547 -0.015008 -0.166098 -2.76274
-0.270539 0.0616157 -0.263959 -0.362205 -0.285429 -2.40993
⋮ ⋮
-0.0470848 -0.030787 -0.201338 -0.0622525 0.165128 -2.64711
0.104699 0.449502 -0.158306 -0.330888 -0.42641 -2.68224
0.0352551 0.164155 0.426519 0.268094 -0.0867661 -2.62129
0.167551 0.00183478 -0.25159 0.447013 -0.256971 -2.97553
-0.00887385 0.741711 0.0900781 0.402043 -0.368223 -3.63044
-0.157661 -0.468791 0.353183 -0.0273369 0.396244 -3.56053
-0.579603 0.00371911 0.229176 0.0521323 0.124366 -3.25365
-0.0275662 -0.154348 -0.15489 0.234171 -0.176908 -3.51783
0.208989 0.227316 -0.06076 -0.152126 0.0325288 -3.18125
-0.122118 0.18972 0.0261586 0.127465 0.221471 -3.53085
0.188565 -0.00411731 -0.0359811 -0.16584 0.230353 -3.37312
0.0176501 -0.3496 -0.0424846 -0.320037 -0.0376885 -3.42145
# Scatter Plot the 2-dim data
p1 = scatter(X_projected[:, 1], X_projected[:, 2], label="Projected Data PC 1 and 2")
p2 = scatter(X_projected[:, 3], X_projected[:, 4], label="Projected Data PC 3 and 4")
p3 = scatter(X_projected[:, 5], X_projected[:, 6], label="Projected Data PC 5 and 6")
plot(p1, p2, p3, title="Scatter Plot of Data", xlabel="x", ylabel="y")
Correlation We can already see that the data is highly correlated when we look at the scatter plots of the data and compate this to the transformed data. We can also see this by looking at the correlation matrix of the data:
# correlation matrix of our data
cor(X)
6×6 Matrix{Float64}:
1.0 0.910001 -0.944146 -0.958282 0.848025 0.966623
0.910001 1.0 -0.89085 -0.90364 0.803305 0.913341
-0.944146 -0.89085 1.0 0.93712 -0.828146 -0.944741
-0.958282 -0.90364 0.93712 1.0 -0.839584 -0.958206
0.848025 0.803305 -0.828146 -0.839584 1.0 0.851197
0.966623 0.913341 -0.944741 -0.958206 0.851197 1.0
We can also look at the correlation matrix of the data transformed by PCA:
# correlation matrix of the transformed data
cor(X_projected)
6×6 Matrix{Float64}:
1.0 -4.93641e-17 1.0664e-15 … -3.01778e-15 3.63898e-15
-4.93641e-17 1.0 4.12116e-15 -2.00991e-15 1.90508e-16
1.0664e-15 4.12116e-15 1.0 -5.38044e-15 4.00714e-15
2.97123e-15 -4.68712e-15 -2.03414e-15 6.04787e-15 -2.34687e-15
-3.01778e-15 -2.00991e-15 -5.38044e-15 1.0 -3.79644e-15
3.63898e-15 1.90508e-16 4.00714e-15 … -3.79644e-15 1.0
Dimensionality Reduction We can reduce the two dimensional data to one dimension using PCA. This is done by projecting the data onto the first principal component as we have done above.
reduced_data = X_projected[:, end] # reduce the data to one dimension
# Plot the reduced data
plot(reduced_data, label="Projected Data", title="Data Projected onto First Principal Component", ylabel="Values of the First Principal Component")
We can actually see that we are able to retrieve the true motion of the ball by looking at the first principal component. The first principal component is a linear combination of the original measurements. We can see the coefficients of this linear combination by looking at the first row of the matrix V. We are actually able to do this since the ball is moving in a straight line. If the ball was moving in a more complicated way, we would not be able to recover the true motion of the ball by looking at the first principal component, since the first principal component is a linear combination of the original measurements.
Variance Explained: We can combine the explained variance with the data to see how much of the variance in the data is explained by each principal component.
# Compute the variance explained by each principal component
var_explained = λ ./ sum(λ)
println("Variance explained by each principal component: ", var_explained)
# Plot the variance explained
p1 = bar(1:length(var_explained), var_explained, title="Variance Explained by Each Principal Component", xlabel="Principal Component", ylabel="Var. Explained")
# Plot the projected data
p2 = plot(X_projected[:, end], label="Projected Data", title="Data Projected onto First Principal Component", ylabel="First PC")
p3 = plot(X_projected[:, 1:end-1], label="Projected Data", title="Data Projected onto First Principal Component", ylabel="Other PCs")
# Display the plots
plot(p1, p2, p3, layout=(3, 1))
Variance explained by each principal component:
[0.008998090151899951, 0.009252262115718415, 0.009403756389051223, 0.009865838947458442, 0.010655619146177834, 0.9518244332496941]
Task 1: Spring Mass with more noise in the data#
Now let’s take a look at what happens when we increase the amount of noise in the data. We will increase the standard deviation of the noise from 0.25 to 2.0. This will make the data much more noisy and harder to analyze. This noise could be viewed as inaccurate measurements for example by chaking of the cameras. Load the following data set and see what the PCA tells you about the system:
amplitude = 2.0 # amplitude of the motion
noise_level = 2.0 # noise level of the data, this can be viewed as chaking of the cameras
X, x_true, y_true = generate_dataset(amplitude, noise_level)
([2.0376071225165697 -0.962083242971589 … 2.363037544747357 1.1716711190249451; 4.598024666096147 0.7079133232123543 … -0.17952621728149665 -2.9958825647708323; … ; 3.915694802924903 0.6057542724125207 … 1.5481391984291801 0.9783823138730351; -1.382304022737907 -2.9793570095662494 … 3.5091889666028377 0.0888789494164488], [2.0, 1.9960534568565431, 1.9842294026289558, 1.9645745014573774, 1.9371663222572622, 1.902113032590307, 1.859552971776503, 1.809654104932039, 1.7526133600877272, 1.6886558510040302 … 1.6886558510040273, 1.7526133600877276, 1.809654104932036, 1.8595529717764974, 1.9021130325903046, 1.9371663222572653, 1.964574501457378, 1.9842294026289553, 1.9960534568565433, 2.0], [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 … 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0])
# Scatter Plot the 2-dim data
p1 = scatter(X[:, 1], X[:, 2], label="Camera 1")
p2 = scatter(X[:, 3], X[:, 4], label="Camera 2")
p3 = scatter(X[:, 5], X[:, 6], label="Camera 3")
plot(p1, p2, p3, title="Scatter Plot of Data", xlabel="x", ylabel="y")
Task 2: Spring Mass with deviations in the y direction#
Now the mass is released off-center so as to produce motion in the \(y\) plane as well as the \(x\) direction. Load the following data set and see what the PCA tells you about the system:
amplitude = 2.0 # amplitude of the motion
noise_level = 0.25 # noise level of the data, this can be viewed as chaking of the cameras
amplitude_deviation = 1.0 # deviation in the y direction
X, x_true, y_true = generate_dataset(amplitude, noise_level, amplitude_deviation)
([1.9899851968101687 0.8448505875873327 … 0.5817802590711504 1.9034306857587504; 1.7273260669009047 0.8711166952591963 … 0.45328303536572856 2.0030463120172306; … ; 2.2546672985756233 0.5336864029747601 … 0.48174554157696003 1.8793727305747103; 1.3845796630932057 0.7326964547509097 … 0.7906766375604602 2.0756748936454428], [2.0, 1.9960534568565431, 1.9842294026289558, 1.9645745014573774, 1.9371663222572622, 1.902113032590307, 1.859552971776503, 1.809654104932039, 1.7526133600877272, 1.6886558510040302 … 1.6886558510040273, 1.7526133600877276, 1.809654104932036, 1.8595529717764974, 1.9021130325903046, 1.9371663222572653, 1.964574501457378, 1.9842294026289553, 1.9960534568565433, 2.0], [6.123233995736766e-17, -0.0627905195293134, -0.12533323356430415, -0.1873813145857246, -0.24868988716485463, -0.30901699437494734, -0.368124552684678, -0.4257792915650727, -0.48175367410171505, -0.5358267949789964 … 0.5358267949789975, 0.48175367410171344, 0.42577929156507427, 0.36812455268468985, 0.30901699437495656, 0.24868988716484725, 0.1873813145857278, 0.12533323356430426, 0.06279051952931021, 7.839596456452825e-15])
# Scatter Plot the 2-dim data
p1 = scatter(X[:, 1], X[:, 2], label="Camera 1")
p2 = scatter(X[:, 3], X[:, 4], label="Camera 2")
p3 = scatter(X[:, 5], X[:, 6], label="Camera 3")
plot(p1, p2, p3, title="Scatter Plot of Data", xlabel="x", ylabel="y")