Dimensionality Reduction¶
This notebook is a revised version of the oustanding lecture notebook by Professor Hug.
Imports¶
As with other notebooks we will use the same set of standard imports.
import numpy as np
import pandas as pd
np.random.seed(23)
import plotly.offline as py
import plotly.express as px
import plotly.graph_objects as go
from plotly.subplots import make_subplots
import plotly.figure_factory as ff
import cufflinks as cf
cf.set_config_file(offline=True, sharing=False, theme='ggplot');
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn import datasets
from sklearn.model_selection import train_test_split
from sklearn.model_selection import cross_val_score
from sklearn.linear_model import LinearRegression
from sklearn.linear_model import LogisticRegression, LogisticRegressionCV
Visualizing Relationships¶
We often create visualizations in order to facilitate exploratory data analysis. For example, we might create scatterplots to explore the relationship between pairs of variables in a dataset. Below, we see weight vs. height.
df = pd.read_csv("data/hongkong_height_weight.csv")
df
px.scatter(df, x="Height", y="Weight")
The dataset below gives the "percentage body fat, age, weight, height, and ten body circumference measurements" for 252 men.
http://jse.amstat.org/v4n1/datasets.johnson.html
For simplicity, we read in only 8 of the provided attributes, yielding the given dataframe.
#http://jse.amstat.org/datasets/fat.txt
df3 = pd.read_fwf("data/fat.dat.txt", colspecs = [(17, 21), (23, 29), (35, 37),
(39, 45), (48, 53), (73, 77),
(80, 85), (88, 93)], header=None, names = ["% fat", "density", "age", "weight", "height", "neck", "chest", "abdomen"])
df3.head()
We see that percentage fat and density in g/cm^3 are almost completely redundant.
px.scatter(df3, x = "% fat", y = "density")
By contrast, while there is a strong correlation between neck and chest measurements, the resulting data is very 2 dimensional.
px.scatter(df3, x = "neck", y = "chest")
Age and height show a small correlation as peolpe seem to get slightly smaller with greater age in this dataset.
px.scatter(df3, x = "age", y = "height")
We note that there is one outlier where a person is slightly less than 29.5 inches tall. While there are a extraordinarily small number of adult males who are less than three feet tall, reflection on the rest of the data from this observation suggest that this was simply an error.
df3.query("height < 40")
df3 = df3.drop(41)
px.scatter(df3, x = "age", y = "height")
We can try to visualize more than 2 attributes at once, but the relationships displayed in e.g. the color and dot size space are much harder for human readers to see. For example, above we saw that density and % fat are almost entirely redundant, but this relationship is impossible to see when comparing the colors and dot sizes.
px.scatter(df3, x = "neck", y = "chest", color="density", size = "% fat")
Seaborn gives us the ability to create a matrix of all possible pairs of variables. This is can be useful, though even with only 8 variables it's still difficult to fully digest.
sns.pairplot(df3)
We can also do this in plotly:
fig = px.scatter_matrix(df3)
fig.update_traces(diagonal_visible=False)
Notice that when you interact with the above plotly plot you actually are interacting with all the plots at once. Try selecting the outlier point.
We should note that despite the very strong relationship between % fat and density, the numerical rank of the data matrix is still 8. For the rank to be 7, we'd need the data to be almost exactly on a line. We'll talk about techniques to reduce the dimensionality over the course of this lecture and the next.
np.linalg.matrix_rank(df3)
House of Representatives Voting Data¶
Next, let's consider voting data from the house of representatives in the U.S. during the month of September 2019. In this example, our goal will be to try to find clusters of representatives who vote in similar ways. For example, we might expect to find that Democrats and Republicans vote similarly to other members of their party.
import yaml
with open("data/legislators-current.yaml", "r") as f:
legislators_data = yaml.safe_load(f)
from datetime import datetime
def to_date(s):
return datetime.strptime(s, '%Y-%m-%d')
legs = pd.DataFrame(
columns=['leg_id', 'first', 'last', 'gender', 'state', 'chamber', 'party', 'birthday'],
data=[[x['id']['bioguide'],
x['name']['first'],
x['name']['last'],
x['bio']['gender'],
x['terms'][-1]['state'],
x['terms'][-1]['type'],
x['terms'][-1]['party'],
to_date(x['bio']['birthday'])] for x in legislators_data])
legs
# February 2019 House of Representatives roll call votes
# Downloaded using https://github.com/eyeseast/propublica-congress
votes = pd.read_csv('data/votes.csv')
votes = votes.astype({"roll call": str})
votes.head()
votes.merge(legs, left_on='member', right_on='leg_id').sample(5)
def was_yes(s):
if s.iloc[0] == 'Yes':
return 1
else:
return 0
vote_pivot = votes.pivot_table(index='member',
columns='roll call',
values='vote',
aggfunc=was_yes,
fill_value=0)
print(vote_pivot.shape)
vote_pivot.head()
vote_pivot.shape
This data has 441 observations (members of the House of Representatives including the 6 non-voting representatives) and 41 dimensions (votes). While politics is quite polarized, none of these columns are linearly dependent as we note below.
np.linalg.matrix_rank(vote_pivot)
Suppose we want to find clusters of similar voting behavior. We might try by reducing our data to only two dimensions and looking to see if we can identify clear patterns. Let's start by looking at what votes were most controversial.
np.var(vote_pivot, axis=0).sort_values(ascending = False)
We see that roll call 548 had very little variance. According to http://clerk.house.gov/evs/2019/roll548.xml, this bill was referring to the 2019 Whistleblower Complaint about President Trump and Ukraine. The full text of the house resolution for this roll call can be found at https://www.congress.gov/bill/116th-congress/house-resolution/576/text:
(1) the whistleblower complaint received on August 12, 2019, by the Inspector General of the Intelligence Community shall be transmitted immediately to the Select Committee on Intelligence of the Senate and the Permanent Select Committee on Intelligence of the House of Representatives; and
(2) the Select Committee on Intelligence of the Senate and the Permanent Select Committee on Intelligence of the House of Representatives should be allowed to evaluate the complaint in a deliberate and bipartisan manner consistent with applicable statutes and processes in order to safeguard classified and sensitive information.
We see that 421 congresspeople voted for this resolution, and 12 did not vote for this resolution. 2 members answered "present" but did not vote no, and 10 did not vote at all. Clearly, a scatterplot involving this particular dimension isn't going to be useful.
vote_pivot['548'].value_counts()
By contrast, we saw high variance for most of the other roll call votes. Most them had variances near 0.25, which is the maximum possible for a variable which can take on values 0 or 1. Let's consider the two highest variance variables, shown below:
vote_pivot['555'].value_counts()
vote_pivot['530'].value_counts()
Let's use these as our two dimensions for our scatterplot and see what happens.
px.scatter(vote_pivot, x='530', y='555')
By adding some random noise, we can get rid of the overplotting.
vote_pivot_jittered = vote_pivot.copy()
vote_pivot_jittered.loc[:, '515':'555'] += np.random.random(vote_pivot_jittered.loc[:, '515':'555'].shape) * 0.3;
px.scatter(vote_pivot_jittered, x='530', y='555')
We can also look at this data labeled by party.
vote_pivot_labeled = vote_pivot.reset_index().merge(legs, left_on='member', right_on='leg_id').set_index('member')
vote_pivot_labeled.head(5)
vote_pivot_labeled_jittered = vote_pivot_labeled.copy()
vote_pivot_labeled_jittered.loc[:, '515':'555'] += np.random.random(vote_pivot_labeled_jittered.loc[:, '515':'555'].shape) * 0.3;
px.scatter(vote_pivot_labeled_jittered, x='530', y='555', color="party",
color_discrete_sequence=["#EF553B", "#00CC96", "#636EFA"])
px.scatter(vote_pivot_labeled_jittered, x='545', y='552', color="party",
color_discrete_sequence=["#EF553B", "#00CC96", "#636EFA"])
We see that considering only two votes does seem to do a pretty good job of telling Republicans from Democrats. We'll see in the next lecture how we can do even better using a technique called "Principle Component Analysis" or PCA.
Before we can get there, we'll need to spend some time reviewing key linear algebra principles.
Matrix Operations in Python¶
age_and_height = np.array([[182, 28], [399, 30], [725, 33]])
M = np.array([[1, 0, 0], [0, 1, 1/12]])
age_and_height @ M
Singular Value Decomposition Experiment¶
Manual Decomposition¶
In the table below, we have the width, height, area, and perimeter of a rectangle stored in a dataframe.
rectangle = pd.read_csv("data/rectangle_data.csv")
rectangle
Naturally the perimeter is just the sum of 2x the width and 2x the height. Thus, if we create a new dataframe that has only the width, height, and area...
rectangle_no_perim = rectangle[["width", "height", "area"]]
rectangle_no_perim.head(5)
... then we can recover the perimeter by multiplying this matrix by
1 0 0 2
0 1 0 2
0 0 1 0transform_3D_to_4D = [[1, 0, 0, 2], [0, 1, 0, 2], [0, 0, 1, 0]]
We can recover the dataframe by multiplying by the transformation matrix and renaming the columns:
dfprod = (rectangle_no_perim @ transform_3D_to_4D)
dfprod.columns = ["width", "height", "area", "perimeter"]
dfprod
Singular Value Decomposition Example¶
Singular value decomposition is a numerical technique to (among other things) automatically uncover such redundancies. Given an input matrix X, SVD will return $U\Sigma$ and $V^T$ such that $ X = U \Sigma V^T $.
u, s, vt = np.linalg.svd(rectangle, full_matrices = False)
As we did before with our manual decomposition, we can recover our original rectangle data by multiplying the three return values of this function back together.
pd.DataFrame(u * s @ vt).head(4)
The two key pieces of the decomposition are $U\Sigma$ and $V^T$, which we can think of for now as analogous to our 'data' and 'transformation operation' from our manual decomposition earlier.
Let's start by looking at $U\Sigma$, which we can compute with the Python code u*s.
pd.DataFrame(u*s)
Similarly, we can look at vt.
pd.DataFrame(vt)
The automatic decomposition returned by the svd function looks quite different than what we got when we manually decomposed our data into "data" and "operations". That is, vt is a bunch of seemingly arbitrary numbers instead of the rather simple:
1 0 0 2
0 1 0 2
0 0 1 0Similarly, if we look at the shape of $U\Sigma$ and $V^T$ we see that they are bigger than in our manual decomposition. Specifically $U\Sigma$ still has 4 columns, meaning that each observation is 4 dimensional. Furthermore, rather than our transformation operation $V^T$ being 3x4, it's 4x4 rows tall, meaning that it maps 4 dimensional inputs back to 4 dimensions.
This seems problematic, because our goal of using SVD was to find a transformation operation that takes 3D inputs and maps them up to 4 dimensions.
Luckily, if we look carefully at $U\Sigma$, we see that the last attribute of each observation is very close to 0.
pd.DataFrame(u * s)
Thus, it makes sense that we remove the last column of $U \Sigma$.
u = u[:, 0:3]
s = s[0:3]
pd.DataFrame(u * s)
Similarly, because the observations are now 3D, we should remove the last row of $V^T$, since we want to use $V^T$ to map our now 3D data into 4D (rather than expecting 4D input data).
vt = vt[0:3, :]
After removing the redundant portions of $U\Sigma$ and $V^T$, we can verify that multiplying them together again yields our original array.
pd.DataFrame(u * s @ vt)
Singular Value Decomposition Experiment¶
The reasons that $U \Sigma$ and $V^T$ look so different than the results of our manual decomposition are a consequence of how singular value decomposition. Specifically, given X, SVD computers $U$, $\Sigma$, and $V$ such that:
- $\Sigma$ is a diagonal matrix containing the singular values of X.
- $U$ and $V$ are matrices whose columns form an orthonormal set.
- $X = U \Sigma V^T$
That is, there are an infinite number of matrices such that $X = AB$. In our example above, we created A and B through manual calculation and insight (recognizing that perimeter was 2w + 2h). SVD computes them automatically, and the results have the specific properties enumerated above.
u, s, vt = np.linalg.svd(rectangle, full_matrices = False)
The Singular Values Matrix $\Sigma$¶
The middle result returned by the SVD process is a diagonal matrix consisting of the singular values. We won't say a lot about what they mean today, but we will note that if the matrix is of rank r, then the first r singular values will be non-zero, and the rest will be zero.
s
Python returns the singular values in a slightly funny format (as a list). To get them into the correct form, we can use diag(s).
np.diag(s)
If we do this, then we can compute the original matrix using only the matrix multiplication operator. That is, rather than writing u * s @ vt, we can do u @ np.diag(s) @ vt.
(u @ np.diag(s) @ vt)[0:5, :]
The fact that the last singular value is zero is why we were able to remove the last dimension of $U$ and the last operation in $V^T$. That is, since the last column of $U \Sigma$ is always zero, we can just remove it completely.
pd.DataFrame(u @ np.diag(s))
(u[:, :-1] @ np.diag(s[:-1]) @ vt[:-1, :])[0:5, :]
$U$ and $V^T$¶
Let's try verifying that $U$ is orthonormal. If so, any column dot producted with itself should be 1, and any columns dot producted with any other should be zero.
np.dot(u[:, 0], u[:, 0])
np.dot(u[:, 1], u[:, 1])
np.dot(u[:, 2], u[:, 2])
for i in range(0, u.shape[1]):
print(f"row {i} dot producted with itself is {np.dot(u[:, i], u[:, i])}")
for i in range(0, u.shape[1]):
for j in range(i + 1, u.shape[1]):
print(f"row {i} dot producted with {j} is {np.dot(u[:, i], u[:, j])}")
Let's now look at $V$.
v = vt.T
for i in range(0, v.shape[1]):
print(f"row {i} of V dot producted with itself is {np.dot(v[:, i], v[:, i])}")
for i in range(0, v.shape[1]):
for j in range(i + 1, v.shape[1]):
print(f"row {i} of V dot producted with {j} is {np.dot(v[:, i], v[:, j])}")
We can also see that the transpose of u (and v) is also the inverse of u (and v).
u.T @ u
v.T @ v