Lec 22: Dimensionality Reduction¶
by Josh Hug (Fall 2019)
# Lec 20: Decision Trees
**by Josh Hug (Fall 2019)**import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
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("hongkong_height_weight.csv")
df.head(5)
sns.scatterplot(data = 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("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.
sns.scatterplot(data = 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.
sns.scatterplot(data = 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.
sns.scatterplot(data = 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)
sns.scatterplot(data = 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.
sns.scatterplot(data = df3, x = "neck", y = "chest", hue="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 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.
from pathlib import Path
from ds100_utils import fetch_and_cache
from datetime import datetime
from IPython.display import display
import yaml
plt.rcParams['figure.figsize'] = (4, 4)
plt.rcParams['figure.dpi'] = 150
sns.set()
base_url = 'https://github.com/unitedstates/congress-legislators/raw/master/'
legislators_path = 'legislators-current.yaml'
f = fetch_and_cache(base_url + legislators_path, legislators_path)
legislators_data = yaml.load(open(f))
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.head(3)
# February 2019 House of Representatives roll call votes
# Downloaded using https://github.com/eyeseast/propublica-congress
votes = pd.read_csv('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.
sns.scatterplot(x='530', y='555', data=vote_pivot);
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
sns.scatterplot(x='530', y='555', data=vote_pivot_jittered);
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
sns.scatterplot(x='530', y='555', data=vote_pivot_labeled_jittered, hue="party");
sns.scatterplot(x='545', y='552', data=vote_pivot_labeled_jittered, hue="party");
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("rectangle_data.csv")
rectangle.head(5)
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 0
transform_3D_to_4D = [[1, 0, 0, 2], [0, 1, 0, 2], [0, 0, 1, 0]]
rectangle_with_perimeter_back = np.array(rectangle_no_perim) @ transform_3D_to_4D
pd.DataFrame(rectangle_with_perimeter_back, columns = ["width", "height", "area", "perimeter"]).head(5)
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.
u*s
Similarly, we can look at vt.
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 0
Similarly, 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.
u * s
Thus, it makes sense that we remove the last column of $U \Sigma$.
u = u[:, 0:3]
s = s[0:3]
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.
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.
(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
Low Rank Approximations of Rectangle Data¶
colnames = list(rectangle.columns)
Suppose we our area and perimeter measurements are slightly noisy.
rectangle_noisy = rectangle.copy()
rectangle_noisy["perimeter"] += np.random.normal(loc = 0, scale = 0.01, size = len(rectangle_noisy["perimeter"]))
rectangle_noisy["area"] += np.random.normal(loc = 0, scale = 0.01, size = len(rectangle_noisy["area"]))
rectangle_noisy.head(5)
u, s, vt = np.linalg.svd(rectangle_noisy, full_matrices = False)
Even on this noisy data, we can recover the exact measurements by multiplying $U \Sigma V^T$.
vt
pd.DataFrame((u @ np.diag(s) @ vt)[0:5, :]).head(5)
Let's take a look at the singular values:
s
Despite not being zero, we can nevertheless again eliminate the last column (dimension) of $U$, last singular value from $\Sigma$ and last row (operation) from $V^T$. We see that we get back reasonable answers, though they're a little off.
pd.DataFrame((u[:, :-1] @ np.diag(s[:-1]) @ vt[:-1, :]),
columns = colnames).head(5)
This is what is known as a "low rank approximation" of the data. Specifically, this is a rank 3 approximation since we're using 3 of the available 4 dimensions. Since the noise was small, this rank 3 approximation is good.
As a warmup to the rest of lecture, let's see what the rank 2 approximation looks like.
Note: In the code below, I've specified which rows/columns that I do want rather than the ones I don't want. That is, instead of writing u[:, :-2], I've written u[:, 0:2].
pd.DataFrame((u[:, 0:2] @ np.diag(s[0:2]) @ vt[0:2, :]),
columns = colnames).head(5)
This code is getting a little hard to read, so let's write a function that computes the rank k approximation of a given DataFrame.
def compute_rank_k_approximation(data, k):
u, s, vt = np.linalg.svd(data, full_matrices = False)
return pd.DataFrame(u[:, 0:k] @ np.diag(s[0:k]) @ vt[0:k, :], columns = data.columns)
pd.DataFrame(compute_rank_k_approximation(rectangle_noisy, 2)).head(5)
We can even look at the rank 1 approximation.
compute_rank_k_approximation(rectangle_noisy, 1).head(5)
In effect, what the SVD process has done is figured out a way to rebuild an approximation of our entire 4D matrix using only 1 column of data. This might seem somewhat shocking.
pd.DataFrame((u[:, 0:1] @ np.diag(s[0:1]) @ vt[0:1, :]),
columns = colnames).head(5)
u, s, vt = np.linalg.svd(rectangle, full_matrices = False)
compute_rank_k_approximation(rectangle_noisy, 1).head(5)