library(HSAUR2)Loading required package: toolsdata(USairpollution)Activity 3.1 - SVD
SUBMISSION INSTRUCTIONS
- Render to html
- Publish your html to RPubs
- Click the Publish button
- Select RPubs
- Follow the prompts for creating an account on RPubs and publishing your solution
- Submit a link to your published solutions
Problem 1
Reconsider the US air pollution data set:
A)
Perform singular value decomposition of this data matrix. Then create the matrix \(D\). Describe what this matrix looks like.
X_raw <- as.matrix(USairpollution)
X <- scale(X_raw, center = TRUE, scale = FALSE)
dim(X) # n x p[1] 41 7n <- nrow(X)
p <- ncol(X)# Singular Value Decomposition
sv <- svd(X)
U <- sv$u
d <- sv$d
V <- sv$v
# Build the matrix
D <- diag(d)
D [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 5053.601 0.0000 0.0000 0.00000 0.00000 0.00000 0.00000
[2,] 0.000 769.7282 0.0000 0.00000 0.00000 0.00000 0.00000
[3,] 0.000 0.0000 167.5661 0.00000 0.00000 0.00000 0.00000
[4,] 0.000 0.0000 0.0000 90.55428 0.00000 0.00000 0.00000
[5,] 0.000 0.0000 0.0000 0.00000 68.32413 0.00000 0.00000
[6,] 0.000 0.0000 0.0000 0.00000 0.00000 21.96092 0.00000
[7,] 0.000 0.0000 0.0000 0.00000 0.00000 0.00000 7.61237D is a diagonal matrix whose diagonal entries are the non-negative singular values of X, sorted from largest to smallest. Off-diagonals are zero.
B)
Verify that \(X=UDV^T\) by plotting all the entries of \(X\) versus all the entries of \(UDV^T\) with the 0/1 line.
X_hat <- U %*% diag(d) %*% t(V)
plot(as.vector(X), as.vector(X_hat),
xlab = "Entries of X", ylab = "Entries of U D V^T",
main = "Check: X vs U D V^T")
abline(0, 1, lwd = 2, col = "gray20")
C)
Consider low-dimensional approximations of the data matrix. What is the fewest number of dimensions required to yield a correlation between the entries of \(X\) and \(\tilde X\) of at least 0.9?
d <- sv$d
k_needed <- which(cumsum(d^2) / sum(d^2) >= 0.81)[1]
k_needed[1] 1D)
Find \(\Sigma\), the covariance matrix of this data set. Then perform eigen-decomposition of this matrix. Verify that
- The eigenvectors of \(\Sigma\) equal the columns of \(V\)
- The eigenvalues of \(\Sigma\) equal the diagonals of \(D^2/(n-1)\)
# Covariance + eigen checks
Sigma <- cov(X_raw)
eg <- eigen(Sigma)
Sigma SO2 temp manu popul wind
SO2 550.947561 -73.560671 8527.7201 6711.9945 3.1753049
temp -73.560671 52.239878 -773.9713 -262.3496 -3.6113537
manu 8527.720122 -773.971341 317502.8902 311718.8140 191.5481098
popul 6711.994512 -262.349634 311718.8140 335371.8939 175.9300610
wind 3.175305 -3.611354 191.5481 175.9301 2.0410244
precip 15.001799 32.862988 -215.0199 -178.0529 -0.2185311
predays 229.929878 -82.426159 1968.9598 645.9860 6.2143902
precip predays
SO2 15.0017988 229.92988
temp 32.8629884 -82.42616
manu -215.0199024 1968.95976
popul -178.0528902 645.98598
wind -0.2185311 6.21439
precip 138.5693840 154.79290
predays 154.7929024 702.59024# Eigenvalues
val_ok <- isTRUE(all.equal(eg$values, (sv$d^2)/(n-1), tolerance = 1e-8))
# Eigenvectors
S <- sign(diag(t(eg$vectors) %*% sv$v)); S[S == 0] <- 1
vec_ok <- max(abs(eg$vectors %*% diag(S) - sv$v)) < 1e-8
c(eigenvalues_match = val_ok, eigenvectors_match = vec_ok) eigenvalues_match eigenvectors_match
TRUE TRUE Yes, they match correctly.
Problem 2
In this problem we explore how “high” a low-dimensional SVD approximation of an image has to be before you can recognize it.
.Rdata objects are essentially R workspace memory snapshots that, when loaded, load any type of R object that you want into your global environment. The command below, when executed, will load three objects into your memory: mysteryU4, mysteryD4, and mysteryV4. These are the first \(k\) vectors and singular values of an SVD I performed on a 700-pixels-tall \(\times\) 600-pixels-wide image of a well-known villain.
load('Data/mystery_person_k4.Rdata')A)
Write a function that takes SVD ingredients u, d and v and renders the \(700 \times 600\) image produced by this approximation using functions from the magick package. Use your function to determine whether a 4-dimensional approximation to this image is enough for you to tell who the mystery villain is. Recall that you will likely need to rescale your recomposed approximation so that all pixels are in [0,1].
library(magick)Linking to ImageMagick 6.9.13.29
Enabled features: cairo, fontconfig, freetype, heic, lcms, pango, raw, rsvg, webp
Disabled features: fftw, ghostscript, x11show_svd_image <- function(u, d, v) {
X_approx <- u %*% diag(d) %*% t(v)
X_scaled <- (X_approx - min(X_approx)) / (max(X_approx) - min(X_approx))
img_raster <- as.raster(X_scaled, max = 1)
plot(as.raster(img_raster))
invisible(X_scaled)
}
show_svd_image(mysteryU4, mysteryD4, mysteryV4)
B)
I’m giving you slightly higher-dimensional approximations (\(k=10\) and \(k=50\), respectively) in the objects below:
load('Data/mystery_person_k10.Rdata')
load('Data/mystery_person_k50.Rdata')
show_svd_image <- function(u, d, v) {
X_approx <- u %*% diag(d) %*% t(v)
X_scaled <- (X_approx - min(X_approx)) / (max(X_approx) - min(X_approx))
plot(as.raster(X_scaled, max = 1))
}
show_svd_image(mysteryU10, mysteryD10, mysteryV10)
show_svd_image(mysteryU50, mysteryD50, mysteryV50)
Create both of the images produced by these approximations. At what point can you tell who the mystery villain is?
At k=10, I can see the face taking shape and hints of the eyes, nose, and mouth, but it’s still not identifiable. At k=50, the features and background are fully clear, and I can recognize the villain. So k=50 is the point where recognition kicks in.
C)
How many numbers need to be stored in memory for each of the following:
- A full \(700\times 600\) image?
700 × 600 = 420,000 numbers SVD approx- 700k + 600k + k= (700+600+1)k = 1301k numbers
- A 4-dimensional approximation?
k = 4: 1301 × 4 = 5,204 numbers
- A 10-dimensional approximation?
k = 10: 1301 × 10 = 13,010 numbers
- A 50-dimensional approximation?
k = 50: 1301 × 50 = 65,050 numbers