1. Introduction

In this exercise, you will implement the K-means clustering algorithm and apply it to find clusters of proteins with similar patterns of presence in blood and cerebrospinal fluid. In a second part, you will use principal component analysis to find a low-dimensional representation of the expression of 692 human genes in a set of 96 RNA-seq samples from different tissues.

To get started with the exercise, you will need to download the starter code and unzip its contents to the directory where you wish to complete the exercise. You will be required to modify several parts of the code, by following the instructions in this assignment. Files needed for this practical include:


2. K-means clustering

In this this exercise, you will implement the K-means algorithm and use it to identify groups of proteins that have similar patterns of presence in blood and cerebrospinal fluid.


2.1. Implementing K-means

The K-means algorithm is a method to automatically cluster similar data examples together. Concretely, you are given a training set \(x^{(i)} \dots x^{(m)}\) (where \(x(i) \in \mathbb{R}^n\)), and want to group the data into a few cohesive “clusters”. The intuition behind K-means is an iterative procedure that starts by guessing the initial centroids, and then refines this guess by repeatedly assigning examples to their closest centroids and then recomputing the centroids based on the assignments.

The K-means algorithm is as follows:

# Initialize centroids
centroids <- kMeansInitCentroids(X, K)
for (iter in 1:iterations){
  # Cluster assignment step: Assign each data point to the
  # closest centroid. idx(i) corresponds to c^(i), the index
  # of the centroid assigned to example i
  idx <- findClosestCentroids(X, centroids)
  
  # Move centroid step: Compute means based on centroid
  # assignments
  centroids <- computeMeans(X, idx, K)
}

The inner-loop of the algorithm repeatedly carries out two steps: i) Assigning each training example \(x^{(i)}\) to its closest centroid, and ii) Recomputing the mean of each centroid using the points assigned to it. The K-means algorithm will always converge to some final set of means for the centroids. Note that the converged solution may not always be ideal and depends on the initial setting of the centroids. Therefore, in practice the K-means algorithm is usually run a few times with different random initializations. One way to choose between these different solutions from different random initializations is to choose the one with the lowest cost function value. You will implement the two phases of the K-means algorithm separately in the next sections.


2.1.1. Finding closest centroids

In the “cluster assignment” phase of the K-means algorithm, the algorithm assigns every training example \(x^{(i)}\) to its closest centroid, given the current positions of centroids. Specifically, for every example \(i\) we set \(c^{(i)} := j\) that minimizes \(|| x^{(i)} - \mu_j ||^2\) where \(c^{(i)}\) is the index of the centroid that is closest to \(x^{(i)}\), and \(\mu_j\) is the position (value) of the \(j\)’th centroid. Note that \(c^{(i)}\) corresponds to idx[i] in the starter code. Your task is to complete the code in findClosestCentroids.R. This function takes the data matrix \(X\) and the locations of all centroids inside centroids and should output a one-dimensional array idx that holds the index (a value in \(\{1, \dots, K\}\), where K is total number of centroids) of the closest centroid to every training example. You can implement this using a loop over every training example and every centroid. Once you have completed findClosestCentroids.R, the following code will run your function. You should see the output [1, 3, 2] corresponding to the centroid assignments for the first 3 examples.

# Load the dataset 
load("data1.RData")

# Select an initial set of centroids
K <- 3 # 3 Centroids
initial_centroids <- matrix(c(3, 3, 6, 2, 8, 5), K, 2,byrow = TRUE)

# Find the closest centroids for the examples using the
# initial_centroids
source("findClosestCentroids.R")
idx <- findClosestCentroids(X, initial_centroids)


2.1.2. Computing centroid means

Given assignments of every point to a centroid, the second phase of the algorithm recomputes, for each centroid, the mean of the points that were assigned to it. Specifically, for every centroid \(k\) we set \[ \mu_k := \frac{1}{|C_k|}\sum\limits_{i \in C_k}x^{(i)} \]

where \(C_k\) is the set of examples that are assigned to centroid \(k\). Concretely, if two examples say \(x^{(3)}\) and \(x^{(5)}\) are assigned to centroid \(k = 2\), then you should update \(\mu_2 = \frac{1}{2}(x^{(3)}+x^{(5)})\). You should now complete the code in computeCentroids.R. You can implement this function using a loop over the centroids. You can also use a loop over the examples; but if you can use a vectorized implementation that does not use such a loop, your code may run faster.

Once you have completed the code in computeCentroids.R, the code below will run your function and output the centroids after the first step of K- means.

#  Compute means based on the closest centroids found in the previous part.
source("computeCentroids.R")
centroids <- computeCentroids(X, idx, K)

centroids
##          [,1]     [,2]
## [1,] 2.428301 3.157924
## [2,] 5.813503 2.633656
## [3,] 7.119387 3.616684


2.2. K-means on example dataset

After you have completed the two functions (findClosestCentroids and computeCentroids), the next step is to run the K-means algorithm on our 2D dataset to help you understand how K-means works. Your functions will be called from inside the runKmeans.R script. Take a look at the script to understand how it works. Notice that the code calls the two functions you implemented in a loop. Then the K-means code will produce a visualization that steps you through the progress of the algorithm at each iteration. Notice how each step of the K-means algorithm changes the centroids and cluster assignments.

#  Compute means based on the closest centroids found in the previous part.
source("runkMeans.R")
max_iters <- 8
kMeans <- runkMeans(X, initial_centroids, max_iters, TRUE)
centroids <- kMeans$centriods
idx <- kMeans$idx


2.3. Random initialization

The initial assignments of centroids for the example dataset were just designed so that you would see the same plot that is shown above. In practice, a good strategy for initializing the centroids is to select random examples from the training set. In this part of the exercise, you should complete the function kMeansInitCentroids.R to randomly select K examples from the dataset, which will become your initial centroids. Make sure that you don’t select the same example twice.

#  Compute means based on the closest centroids found in the previous part.
source("kMeansInitCentroids.R")
initial_centroids <- kMeansInitCentroids(X, K)
kMeans <- runkMeans(X, initial_centroids, max_iters, TRUE)
centroids <- kMeans$centroids
idx <- kMeans$idx


3. Principal components analysis

In this exercise, you will use principal component analysis (PCA) to perform dimensionality reduction. You will first experiment with an example 2D dataset to get intuition on how PCA works, and then use it on a bigger RNA-seq dataset that contains the expression of 692 human genes measured in 96 samples from 3 different tissues.


3.1. Example dataset

To help you understand how PCA works, you will first start with a 2D dataset which has one direction of large variation and one of smaller variation. The code below will plot the training data. In this part of the exercise, you will visualize what happens when you use PCA to reduce the data from 2D to 1D. In practice, you might want to reduce data from 10,000 to 10 dimensions, say; but using lower dimensional data in this example allows us to visualize the algorithms better.

load("data2.RData")
plot(X, col='blue', asp = 1)


3.2. Implementing PCA

In this part of the exercise, you will implement PCA. PCA consists of two computational steps: First, you compute the covariance matrix of the data. Second, you use R svd function to compute the eigenvectors \(u_1, u_2, \dots, u_n\). These will correspond to the principal components of variation in the data. Before using PCA, it is important to first center the variables by subtracting the mean values, and scaling each dimension so that they are in the same range. In the code below this step has been performed for you using the featureScaling.R function. After centering and scaling the data, you can run PCA to compute the principal components. You task is to complete the code in pca.R to compute the principal components of the dataset. As stated above, first, you should compute the covariance matrix of the data, which is given by:

\[ \Sigma = \frac{1}{m}X^TX \]

where \(X\) is the data matrix with examples in rows, and \(m\) is the number of examples. Note that \(\Sigma\) is a \(n \times n\) matrix and not the summation operator. After computing the covariance matrix, you can perform singular value decomposition (SVD) on it to compute the principal components. In R, you can do it using the function svd. The SVD output will contain three objects: \(d\), \(U\) and \(V\). \(U\) will contain the principal components and \(d\) will be the diagonal elements of a diagonal matrix.

Once you have completed pca.R, the following code will run PCA on the example dataset and plot the corresponding principal components found. The script will also output the top principal component (eigenvector) found, and you should expect to see an output of about [-0.707, -0.707]. (It is possible that R may instead output the negative of this, since \(u_1\) and \(-u_1\) are equally valid choices for the first principal component.)

source("featureScaling.R")
source("pca.R")

#  Before running PCA, it is important to first normalize X
fS <- featureScaling(X)
Xnorm <- fS$Xnorm
mu <- fS$mu
sigma <- fS$sigma

#  Run PCA
pca_out <- pca(Xnorm)
U <- pca_out$U
D <- pca_out$D

#  Compute mu, the mean of the each feature

#  Draw the eigenvectors centered at mean of data. These lines show the
#  directions of maximum variations in the dataset.
plot(X, col='blue', asp = 1)
lines(rbind(t(mu), mu + D[1,1]*t(U[,1])), lwd = 2)
lines(rbind(t(mu), mu + D[2,2]*t(U[,2])), lwd = 2)


3.3. Dimensionality reduction with PCA

After computing the principal components, you can use them to reduce the feature dimension of your dataset by projecting each example onto a lower dimensional space, \(x^{(i)} \to z^{(i)}\) (e.g., projecting the data from 2D to 1D). In this part of the exercise, you will use the eigenvectors returned by PCA and project the example dataset into a 1-dimensional space. In practice, if you were using a learning algorithm such as linear regression or a neural network, you could now use the projected data instead of the original data. By using the projected data, you can train your model faster as there are less dimensions in the input.


3.3.1. Projecting the data onto the principal components

You should now complete the code in projectData.R. Specifically, you are given a dataset \(X\), the principal components \(U\), and the desired number of dimensions to reduce to \(K\). You should project each example in \(X\) onto the top \(K\) components in \(U\). Note that the top \(K\) components in \(U\) are given by the first \(K\) columns of \(U\), that is U_reduce = U[, 1:K]. Once you have completed the code in projectData.R, the code below will project the first example onto the first dimension and you should see a value of about 1.481 (or possibly -1.481, if you got \(-u_1\) instead of \(u_1\)).

#  Project the data onto K = 1 dimension
source("projectData.R")
K <- 1
Z <- projectData(Xnorm, U, K)
cat(sprintf('Projection of the first example: %f\n', Z[1]))
## Projection of the first example: 1.481274


3.3.2. Reconstructing an approximation of the data

After projecting the data onto the lower dimensional space, you can approximately recover the data by projecting them back onto the original high dimensional space. Your task is to complete recoverData.R to project each example in \(Z\) back onto the original space and return the recovered approximation in \(X_rec\). Once you have completed the function, you should recover an approximation of the first example about [-1.047, -1.047].

#  Reconstruct an approximation of the data
source("recoverData.R")
X_rec  <- recoverData(Z, U, K)
cat(sprintf('Approximation of the first example: %f %f\n', X_rec[1,1], X_rec[1,2]))
## Approximation of the first example: -1.047419 -1.047419


3.3.3. Visualizing the projections

After completing both projectData and recoverData, the following code will now perform both the projection and approximate reconstruction to show how the projection affects the data. The original data points are indicated with the blue circles, while the projected data points are indicated with the red circles. The projection only retains the information in the direction given by \(u_1\).

plot(Xnorm, col='blue', ylim = c(-2.5, 2.5), asp = 1)
points(X_rec, col = 'red')
for (i in 1:nrow(Xnorm)){
  lines(rbind(Xnorm[i, ], X_rec[i, ]))
}


3.4. Visualizing RNA-seq data with PCA

Now that you have implemented the PCA algorithm, use it to plot a low-dimensional (3D) representation of the expression of 692 human genes (692D) in a set of 96 RNA-seq samples from 3 different tissues (brain, liver and heart). Visualizing a dataset in a low-dimensional space often helps to understand the structure of the data (e.g. presence/absence of a batch effect, samples that come from different experiments, tissues or conditions, etc.).

# Load dataset
load("data3.RData")

# Perform PCA
# That's yours ;)

# Plot in 3D, colored by tissue
library(plotly)
plot_ly(x = Z[,1], y = Z[,2], z = Z[,3], 
        color = rownames(X), 
        marker = list(size = 5)) %>% 
  add_markers()  %>% 
  layout(scene = list(xaxis = list(title = "Z[,1]"), 
                      yaxis = list(title = "Z[,2]"),
                      zaxis = list(title = "Z[,3]")))