library(HSAUR2)Loading required package: toolsdata(USairpollution)Activity 3.1 - SVD
SUBMISSION INSTRUCTIONS
- Render to html
- Publish your html to RPubs or Posit Connect
- Submit a link to your published solutions
Problem 1
A)
Consider plots of \(UD_i\) for four different \(D_i\):
Match the \(D\) to the plot.
- \(D_1 = \begin{pmatrix} 3 & 0 \\ 0 & 2 \end{pmatrix}\) = B
- \(D_2 = \begin{pmatrix} 3 & 0 \\ 0 & 1 \end{pmatrix}\) = C
- \(D_3 = \begin{pmatrix} 1 & 0 \\ 0 & 0.5 \end{pmatrix}\) = D
- \(D_4 = \begin{pmatrix} 2 & 0 \\ 0 & 1 \end{pmatrix}\) = A
B)
Consider two different \(D\) matrices:
\[D_1 = \begin{pmatrix} 4 & 0\\ 0 & 1 \end{pmatrix}\] \[D_2 = \begin{pmatrix} 2 & 0\\ 0 & 1 \end{pmatrix}\]
and two different \(V\) matrices:
\[V_1 = \begin{pmatrix} 0.894 & -0.447\\ 0.447 & 0.894 \end{pmatrix}\] \[V_2 = \begin{pmatrix} -0.447 & 0.894\\ 0.894 & 0.447 \end{pmatrix}\]
Consider the following matrix product:
\[UD_i V_j^T, i\in \{1,2\}, j \in \{1,2\}\]
In other words, \(U\) is multiplied by \(D_i V_j^T\) for all four combinations of \(D_i\) and \(V_j\). Match the product to the plot.
A = U * D1 V1
B = U * D1 V2
C = U * D2 V1
D = U * D2 V2
Problem 2
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.
library(tidyverse)── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr 1.2.1 ✔ readr 2.2.0
✔ forcats 1.0.1 ✔ stringr 1.6.0
✔ ggplot2 4.0.3 ✔ tibble 3.3.1
✔ lubridate 1.9.5 ✔ tidyr 1.3.2
✔ purrr 1.2.2
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorsdecompose_air <- svd(as.matrix(USairpollution))
U <- decompose_air$u
D <- decompose_air$d
V <- decompose_air$v
(D_mat <- diag(D)) %>% round(1) [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 7051.9 0.0 0.0 0.0 0.0 0.0 0.0
[2,] 0.0 931.1 0.0 0.0 0.0 0.0 0.0
[3,] 0.0 0.0 540.5 0.0 0.0 0.0 0.0
[4,] 0.0 0.0 0.0 92.7 0.0 0.0 0.0
[5,] 0.0 0.0 0.0 0.0 85.2 0.0 0.0
[6,] 0.0 0.0 0.0 0.0 0.0 52.9 0.0
[7,] 0.0 0.0 0.0 0.0 0.0 0.0 10.17 by 7 matrix with values along the diagonal in descending order based on value. Each diagonal entry is a singular value, and it measures how much the data spreads out along one particular direction.
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 <- as.matrix(USairpollution)
Xrecomposed <- U %*% diag(D) %*% t(V)
plot(X, Xrecomposed); abline(0, 1)
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?
approx_cor <- \(k) {
Xtilde <- U[,1:k, drop = FALSE] %*% diag(D[1:k], nrow = k) %*% t(V[,1:k, drop = FALSE])
return(cor(as.vector(X), as.vector(Xtilde)))
}
(cor_table <- tibble(k = 1:7)
%>% mutate(correlation = map_dbl(k, approx_cor))
)# A tibble: 7 × 2
k correlation
<int> <dbl>
1 1 0.986
2 2 0.997
3 3 1.000
4 4 1.000
5 5 1.000
6 6 1.000
7 7 1 K=1 works w/ correlation 0.986
D)
Find \(\Sigma\), the covariance matrix of this data set. Then perform eigen-decomposition of this matrix. Also, mean-center (but don’t unit-scale) \(X\): use scale(X, center = TRUE, scale = FALSE). SVD this mean-centered data matrix. Verify that:
- The eigenvectors of \(\Sigma\) equal the columns of \(V\);
- The eigenvalues of \(\Sigma\) equal the diagonals of \(D^2/(n-1)\).
n <- nrow(X)
(Sigma <- cov(X)) %>% round(1) SO2 temp manu popul wind precip predays
SO2 550.9 -73.6 8527.7 6712.0 3.2 15.0 229.9
temp -73.6 52.2 -774.0 -262.3 -3.6 32.9 -82.4
manu 8527.7 -774.0 317502.9 311718.8 191.5 -215.0 1969.0
popul 6712.0 -262.3 311718.8 335371.9 175.9 -178.1 646.0
wind 3.2 -3.6 191.5 175.9 2.0 -0.2 6.2
precip 15.0 32.9 -215.0 -178.1 -0.2 138.6 154.8
predays 229.9 -82.4 1969.0 646.0 6.2 154.8 702.6eigenSigma <- eigen(Sigma)
(X_centered <- X
%>% scale(center = TRUE, scale = FALSE)
) %>% head SO2 temp manu popul wind precip
Albany 15.9512195 -8.1634146 -419.09756 -492.6098 -0.6439024 -3.4090244
Albuquerque -19.0487805 1.0365854 -417.09756 -364.6098 -0.5439024 -28.9990244
Atlanta -6.0487805 5.7365854 -95.09756 -111.6098 -0.3439024 11.5709756
Baltimore 16.9512195 -0.7634146 161.90244 296.3902 0.1560976 4.5409756
Buffalo -19.0487805 -8.6634146 -72.09756 -145.6098 2.9560976 -0.6590244
Charleston 0.9512195 -0.5634146 -428.09756 -537.6098 -2.9439024 3.9809756
predays
Albany 21.097561
Albuquerque -55.902439
Atlanta 1.097561
Baltimore -2.902439
Buffalo 52.097561
Charleston 34.097561svdX <- svd(X_centered)
eigenSigma$vectors %>% round(3) [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] 0.017 -0.100 0.209 0.959 -0.152 -0.054 -0.027
[2,] -0.001 0.026 -0.072 -0.110 -0.478 -0.853 -0.164
[3,] 0.697 -0.710 -0.067 -0.073 -0.010 -0.002 0.001
[4,] 0.717 0.693 0.057 0.049 0.011 0.003 0.000
[5,] 0.000 -0.001 0.005 -0.015 0.025 0.177 -0.984
[6,] 0.000 0.001 0.266 -0.163 -0.833 0.453 0.064
[7,] 0.003 -0.069 0.934 -0.185 0.233 -0.184 -0.019svdX$v %>% round(3) [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,] -0.017 0.100 0.209 -0.959 0.152 -0.054 -0.027
[2,] 0.001 -0.026 -0.072 0.110 0.478 -0.853 -0.164
[3,] -0.697 0.710 -0.067 0.073 0.010 -0.002 0.001
[4,] -0.717 -0.693 0.057 -0.049 -0.011 0.003 0.000
[5,] 0.000 0.001 0.005 0.015 -0.025 0.177 -0.984
[6,] 0.000 -0.001 0.266 0.163 0.833 0.453 0.064
[7,] -0.003 0.069 0.934 0.185 -0.233 -0.184 -0.019eigenSigma$values %>% round(3)[1] 638471.962 14812.038 701.960 205.002 116.705 12.057 1.449(svdX$d^2/(n-1)) %>% round(3)[1] 638471.962 14812.038 701.960 205.002 116.705 12.057 1.449The eigenvalues of Σ exactly equal D²/(n − 1). The eigenvectors of Σ equal the columns of V up to sign: some columns are multiplied by −1, but since an eigenvector and its negative define the same axis, they represent the same directions.
Problem 3: Spooky scary zombie professor
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 certain WSU professor who has turned into a zombie. How low-dimensional an approximation of the full image can you get away with, and still identify the professor?
load('Data/zombie_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 has been zombified. 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, x11render_svd <- \(u, d, v) {
recomposed <- u %*% diag(d, nrow = length(d)) %*% t(v)
R <- max(recomposed) - min(recomposed)
recomposed_scaled <- (recomposed - min(recomposed))/R
recomposed_image <- (recomposed_scaled
%>% as.raster
%>% image_read
)
return(recomposed_image)
}
render_svd(mysteryU4, mysteryD4, mysteryV4)
4-dimensional approximation allows us to begin to make out the first person (maybe someone with a white beard?? This still has quite a bit of work to do to be a clear picture.
B)
I’m giving you slightly higher-dimensional approximations (\(k=10\) and \(k=50\), respectively) in the objects below:
load('Data/zombie_k10.Rdata')
load('Data/zombie_k50.Rdata')
render_svd(mysteryU10, mysteryD10, mysteryV10)
render_svd(mysteryU50, mysteryD50, mysteryV50)
Create both of the images produced by these approximations. At what point can you tell who has been zombified?
K = 10, still is quite grainy and is tough to make out. K = 50 gives essentially a clear picture.
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
A 4-dimensional approximation?
700 * 4 + 4 + 4 * 600 = 5,204
A 10-dimensional approximation?
700 × 10 + 10 + 10 × 600 = 13,010
A 50-dimensional approximation?
700 × 50 + 50 + 50 × 600 = 65,050