# A tibble: 60,000 × 785
X0 X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 X11 X12
<dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 0 0 0 0 0 0 0 0 0 0 0 0 0
2 0 0 0 0 0 0 0 0 0 0 0 0 0
3 0 0 0 0 0 0 0 0 0 0 0 0 0
4 0 0 0 0 0 0 0 0 0 0 0 0 0
5 0 0 0 0 0 0 0 0 0 0 0 0 0
6 0 0 0 0 0 0 0 0 0 0 0 0 0
7 0 0 0 0 0 0 0 0 0 0 0 0 0
8 0 0 0 0 0 0 0 0 0 0 0 0 0
9 0 0 0 0 0 0 0 0 0 0 0 0 0
10 0 0 0 0 0 0 0 0 0 0 0 0 0
# ℹ 59,990 more rows
# ℹ 772 more variables: X13 <dbl>, X14 <dbl>, X15 <dbl>, X16 <dbl>, X17 <dbl>,
# X18 <dbl>, X19 <dbl>, X20 <dbl>, X21 <dbl>, X22 <dbl>, X23 <dbl>,
# X24 <dbl>, X25 <dbl>, X26 <dbl>, X27 <dbl>, X28 <dbl>, X29 <dbl>,
# X30 <dbl>, X31 <dbl>, X32 <dbl>, X33 <dbl>, X34 <dbl>, X35 <dbl>,
# X36 <dbl>, X37 <dbl>, X38 <dbl>, X39 <dbl>, X40 <dbl>, X41 <dbl>,
# X42 <dbl>, X43 <dbl>, X44 <dbl>, X45 <dbl>, X46 <dbl>, X47 <dbl>, …Unsupervised Learning
ACTL3142 & ACTL5110 Statistical Machine Learning for Risk Applications
Some of the figures in this presentation are taken from "An Introduction to Statistical Learning, with applications in R" (Springer, 2013) with permission from the authors: G. James, D. Witten, T. Hastie and R. Tibshirani
K-means clustering: Local Optima
- The algorithm finds a local rather than a global optimum
- Results depend on initial centroids used
- Important to run the algorithm multiple times and select the best solution (minimum within-cluster variation)
What is the right value of K?
MNIST Dataset
The MNIST dataset.
The data is flattened
An image turned into a vector.
Preparation
Take just a fraction of the data, and make a plotting function.
# Separate the features and the labels
x <- as.matrix(train_df[, -785])
y <- train_df$label
# Split the data into train/validation and test sets
set.seed(88)
test_indices <- sample(1:nrow(x),
size = 0.2*nrow(x))
x_test <- x[test_indices,]
y_test <- y[test_indices]
x_train_val <- x[-test_indices,]
y_train_val <- y[-test_indices]
train_indices <- sample(1:nrow(x_train_val),
size = 0.75*nrow(x_train_val))
x_train <- x_train_val[train_indices,]
y_train <- y_train_val[train_indices]
x_val <- x_train_val[-train_indices,]
y_val <- y_train_val[-train_indices]# Convert to NumPy arrays
train_df = pd.read_csv("mnist_train.csv")
x = train_df.iloc[:, :-1].to_numpy()
y = train_df.iloc[:, -1].to_numpy().flatten()
# Split the data into train/validation and test sets
x_train_val, x_test, y_train_val, y_test = train_test_split(x, y, test_size=0.2, random_state=0)
x_train, x_val, y_train, y_val = train_test_split(x_train_val, y_train_val, test_size=0.25, random_state=0)
x_train = x_train[:1000]
y_train = y_train[:1000]
# Plot a single digit
def plot_digit(digit, ax):
digit_matrix = 1 - digit.reshape((28, 28))
ax.imshow(digit_matrix, cmap='gray',
vmin=0, vmax = 1,
interpolation='none')
ax.set_xticks([])
ax.set_yticks([]) 
Plotting the data
There are 10 natural clusters
K-means Clustering on MNIST I
The within-cluster variation is
K-means Clustering on MNIST II
The within-cluster variation is
K-means Clustering on MNIST III
The within-cluster variation is
Elbow method MNIST I
Code
KMeans(n_clusters=20)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
KMeans(n_clusters=20)
Elbow method MNIST II
Code
KMeans(n_clusters=200)In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.
KMeans(n_clusters=200)
Hierarchical Clustering: The dendogram
Choose a max distance I
Choose a max distance II
Choose a max distance III
Same Data, Different Linkage
Average, complete, and single linkage applied to an example data set. Average and complete linkage tend to yield more balanced clusters.
PCA Motivation: Data compression I
Code
set.seed(2)
x1 <- 1:10
x2 <- 1:10 + runif(n = 10)*4
dataPCA <- data.frame(x1 = x1 - mean(x1), x2 = x2 - mean(x2), col = factor(1:10))
# Plot simple PCA
PCA0 <- ggplot(dataPCA) + geom_point(aes(x = x1, y = x2, colour = col),
shape = 16, size = 3) +
theme_bw() +
theme(axis.text=element_blank(),
axis.ticks=element_blank(),
legend.position = "none",
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line = element_line(arrow = arrow())) +
coord_equal() +
labs(x = expression("Math Skill - " * x[1]), y = expression("Math Enjoyment - " * x[2])) +
xlim(-6,6) + ylim(-6,6)
PCA0
PCA Motivation: Data compression II
PCA Motivation: Data compression III
Code
dataPCA <- dataPCA %>% mutate(z1 = pca.out$x[,1],
z2 = pca.out$x[,2],
proj1 = z1/sqrt(pca.out$rotation[1,1]^2 +
pca.out$rotation[2,1]^2)*pca.out$rotation[1,1],
proj2 = z1/sqrt(pca.out$rotation[1,1]^2 +
pca.out$rotation[2,1]^2)*pca.out$rotation[2,1])
(PCA2 <- PCA1 + geom_segment(data = dataPCA, aes(x = proj1, y= proj2, xend = x1, yend = x2)))
PCA Motivation: Data compression IV
Code
PROJ1 <- ggplot(dataPCA) + geom_point(aes(x = z1, y = 1, colour = col),
shape = 16, size = 3) +
theme_bw() +
theme(axis.text=element_blank(),
axis.ticks=element_blank(),
legend.position = "none",
panel.border = element_blank(),
panel.grid = element_blank(),
axis.line.x = element_line(arrow = arrow()),
axis.line.y = element_line(arrow = arrow(), colour = "white"),
axis.title.y = element_text(colour = "white")) +
coord_equal() +
labs(x = expression("Math Aptitude - " * z[1]), y = "asfdasf") +
ylim(0.5,1.5) + xlim(-1 -6*(pca.out$rotation[1,1] + pca.out$rotation[2,1]),
6*(pca.out$rotation[1,1] + pca.out$rotation[2,1]))
PROJ1
2 principal component directions I
Ninety observations simulated in three dimensions. The observations are displayed in color for ease of visualization. The first two principal component directions span the plane that best fits the data. The plane is positioned to minimize the sum of squared distances to each point.
2 principal component directions II
The first two principal component score vectors give the coordinates of the projection of the 90 observations onto the plane.
How many principal components to use?
Left: a scree plot depicting the proportion of variance explained by each of the four principal components in the USArrests data. Right: the cumulative proportion of variance explained by the four principal components in the USArrests data.
PCA on MNIST
Code
([<matplotlib.axis.XTick object at 0x76c0137ff090>, <matplotlib.axis.XTick object at 0x76c0137fecd0>, <matplotlib.axis.XTick object at 0x76c013a1e8d0>, <matplotlib.axis.XTick object at 0x76c013acb0d0>, <matplotlib.axis.XTick object at 0x76c013b08850>], [Text(0, 0, '0'), Text(200, 0, '200'), Text(400, 0, '400'), Text(600, 0, '600'), Text(800, 0, '800')])
Scree plot
A scree plot is a line plot of the eigenvalues of principal components on the y-axis and the factors on the x-axis.
The scree plot is used to determine the number of factors to retain in a principal components analysis.
The point at which the line starts to level off is the number of factors to retain.
PCA on MNIST: Cumulative
Code
Autoencoder
An autoencoder takes an observation, maps it to a latent space via an encoder module, then decodes it back to an output with the same dimensions via a decoder module.
Schematic of an autoencoder.
Code
# Perform PCA on non-zero variance columns
non_zero_var_cols <- which(col_std_devs != 0)
pca <- prcomp(x_train[, non_zero_var_cols], center = TRUE, scale. = TRUE)
pca_train_full <- predict(pca, newdata = x_train[, non_zero_var_cols])
indices <- c()
for (digit in 0:9) {
for (i in seq_along(y_train)) {
if (y_train[i] == digit) {
indices <- c(indices, i)
break
}
}
}
pca_autoencoder <- function(pca_train_full, num_components) {
pca_train <- pca_train_full[, 1:num_components]
# Inverse transform from the PCA reduced space back to the original space
pca_reconstructed <- pca_train %*% t(pca$rotation[, 1:num_components])
pca_reconstructed <- scale(pca_reconstructed, center = -pca$center, scale = FALSE)
# Create a full reconstructed dataset with zero-variance columns
pca_reconstructed_full <- matrix(0, nrow = nrow(x_train), ncol = 784)
pca_reconstructed_full[, non_zero_var_cols] <- pca_reconstructed
pca_reconstructed_full <- pmin(pmax(pca_reconstructed_full, 0), 1)
pca_reconstructed_full
}
plot_reconstructions <- function(x_train, pca_reconstructed_full, indices) {
par(mfrow = c(2, 5))
for (index in indices[1:5]) {
plot_digit(x_train[index, ])
}
for (index in indices[1:5]) {
plot_digit(pca_reconstructed_full[index,])
}
}PCA on MNIST: Reconstructed with 400 PCs
PCA on MNIST: Reconstructed with 100 PCs
PCA on MNIST: Reconstructed with 25 PCs
PCA on MNIST: Reconstructed with 10 PCs
PCA on MNIST: Reconstructed with 5 PCs
PCA on MNIST: Reconstructed with 400 PCs II
PCA on MNIST: Reconstructed with 100 PCs II
PCA on MNIST: Reconstructed with 25 PCs II
PCA on MNIST: Reconstructed with 10 PCs II
PCA on MNIST: Reconstructed with 5 PCs II
Just pull out the 4s and the 7s
Code
# Pull out the 4s and 7s
idxs <- which(y_train %in% c(4,7))
y_train <- y_train[idxs] == 7
x_train <- x_train[idxs,]
# Prepare validation data
idxs_val <- which(y_val %in% c(4, 7))
y_val <- y_val[idxs_val] == 7
x_val <- x_val[idxs_val,]
# Prepare validation data
idxs_test <- which(y_test %in% c(4, 7))
y_test <- y_test[idxs_test] == 7
x_test <- x_test[idxs_test,]
col_std_devs <- apply(x_train, 2, sd)
pr_comp_train <- prcomp(x_train)
# Print the first 10 digits in x_train
par(mfrow = c(2,5))
for (i in 1:10) {
plot_digit(x_train[i,])
}
Code
# Pull out the 4s and 7s
idxs = np.where(np.isin(y_train, [4, 7]))
y_train = y_train[idxs] == 7
x_train = x_train[idxs]
# Do the same for x_val, y_val, x_test, y_test
idxs_val = np.where(np.isin(y_val, [4, 7]))
y_val = y_val[idxs_val] == 7
x_val = x_val[idxs_val]
idxs_test = np.where(np.isin(y_test, [4, 7]))
y_test = y_test[idxs_test] == 7
x_test = x_test[idxs_test]
col_std_devs = np.std(x_train, axis=0)
pca = PCA()
pr_out = pca.fit(x_train)Compression
“A photograph, which used to be a pattern of pigment on a sheet of chemically coated paper, is now a string of numbers, each one representing the brightness and color of a pixel. An image captured on a 4-megapixel camera is a list of 4 million numbers-no small commitment of memory for the device shooting the picture. But these numbers are highly correlated with each other. If one pixel is bright green, the next one over is likely to be as well. The actual information contained in the image is much less than 4 million numbers’ worth-and it’s precisely this fact that makes it possible to have compression, the critical mathematical technology that allows images, videos, music, and text to be stored in much smaller spaces than you’d think.” Jordan Ellenberg
