1.1.1 Finding closest centroids

ex7data2 <- read.csv("data/ex7data/ex7data2.csv")

g1 <- ggplot(ex7data2, aes(V1, V2)) + geom_point()
g1

This part is fairly straightforward, but it’s broken into parts to make it easier to debug. minEuDist takes a coordinate and a vector of centroid coordinates and returns the index of the nearest one.

minEuDist <- function(xi, centroids){
    euDist <- apply(centroids, 1, function(j){dist(rbind(j, xi))})
    return(which.min(euDist))
}

initial_centroids <- matrix(c(3, 3, 6, 2, 8, 5), 3, 2, TRUE)
minEuDist(ex7data2[15,], initial_centroids)

## [1] 1

findClosestCentroids <- function(X, centroids){
    apply(X, 1, function(i){minEuDist(i, centroids)})
}

findClosestCentroids applies minEuDist over an array.

We should see the output [1 3 2] corresponding to the centroid assignments for the first 3 examples.

idx <- findClosestCentroids(ex7data2, initial_centroids)
idx[1:3]

## [1] 1 3 2

1.1.2 Computing centroid means

The other part of K-means clustering is updating the value of the centroids. This function takes the array of points, their initial assignment to a centroid, and the positions of the centroids. The data are then subsetted by centroid, and that centroid is moved to the new average value of those points that are assigned to it.

computeCentroids <- function(X, idx, K){
    newCentroids <- K
    for(k in 1:nrow(K)){
        Xk <- X[k == idx,]
        newCentroids[k, ] <- apply(Xk, 2, mean)
    }
    return(newCentroids)
}

computeCentroids(ex7data2, idx, initial_centroids)

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

1.4.1 \(K\)-means on pixels

bird_small <- t(read.csv("data/ex7data/bird_small.csv"))
birdInit <- randomInit(bird_small, 16)

kMeansBird <- kMeans(bird_small, birdInit, 15)

finalCentroid <- kMeansBird[nrow(kMeansBird),(ncol(kMeansBird)-16*3+1):ncol(kMeansBird)]
finalCentroid <- round(matrix(finalCentroid, nrow = 16, ncol = 3, byrow = TRUE))
compressedImage <- matrix(kMeansBird[nrow(kMeansBird), 1:(128*128)])

If we take compressedImage and finalCentroid together, this provides all the information needed. To reconstruct the image, we map it back into standard RGB values. Note that these are uncompressed.

r <- sapply(compressedImage, function(x){finalCentroid[x, 1]})
g <- sapply(compressedImage, function(x){finalCentroid[x, 2]})
b <- sapply(compressedImage, function(x){finalCentroid[x, 3]})

col <- rgb(r, g, b, maxColorValue = 255)
dim(col) <- c(128, 128)

grid.raster(col, interpolate=FALSE)

1.4.2 Optional (ungraded) exercise: Use your own image

2.2.1 Plotting the eigenvectors of the scatterplot

The eigenvectors are centered at the mean of the data, and extend to the product of the singular values and left singular vectors. I again got some help from kaleko and he multiplied the endpoint of the vectors by 1.5, which I’ve done here as well. I’m not sure why that’s the case, but it appears to match the figure in the assignment sheet. Regardless, for dimensionality reduction, only the direction of the vector matters.

We plot these figures with the same scale on both axes so the eigenvectors appear perpendicular.

ex71means <- apply(ex7data1, 2, mean)

g4 <- g3 + geom_segment(x = ex71means[1], 
                 xend = ex71means[1] + 1.5 * ex71svd$d[1] * ex71svd$u[1,1], 
                 y = ex71means[2], 
                 yend = ex71means[2] + 1.5 * ex71svd$d[1] * ex71svd$u[1,2]) +
    geom_segment(x = ex71means[1], 
                 xend = ex71means[1] + 1.5 * ex71svd$d[2] * ex71svd$u[2,1], 
                 y = ex71means[2], 
                 yend = ex71means[2] + 1.5 * ex71svd$d[2] * ex71svd$u[2,2])
g4

2.3.1 Projecting the data onto the principal components

projectData <- function(X, U, K){
    X %*% U[, 1:K]
}

We should see a value of about 1.481.

Z <- projectData(ex71normd, ex71svd$u, 1)
Z[1]

## [1] 1.481274

2.3.2 Reconstructing an approximation of the data

recoverData <- function(Z, U, K){
    Z %*% t(U[, 1:K])
}

We should see a value of about [-1.047 -1.047].

X_rec <- recoverData(Z, ex71svd$u, 1)
X_rec[1,]

## [1] -1.047419 -1.047419

2.3.3 Visualizing the projections

m <- as.data.frame(ex71normd)
m$pc <- "Original"

n <- as.data.frame(X_rec)
n$pc <- "Projected"

ex71projected <- rbind(m, n)
rm(m)
rm(n)

ex71projected$pc <- as.factor(ex71projected$pc)
ex71projected$match <- 1:50

g5 <- ggplot(ex71projected, aes(V1, V2, group = match)) +
    geom_point(aes(color = pc)) +
    geom_line(linetype = "dotdash") +
    theme(legend.position="none")
g5

2.4.1 PCA on Faces

facesScaled <- scale(ex7faces)
facesSvd <- svd(cov(facesScaled))
facesProjectd <- projectData(X = facesScaled, U = facesSvd$u, K = 100)

2.4.2 Dimensionality Reduction

facesRecovered <- recoverData(facesProjectd, facesSvd$u, 100)
plotFaces(facesRecovered)