Activity 3.1 - SVD

SUBMISSION INSTRUCTIONS

  1. Render to html
  2. Publish your html to RPubs
  1. Submit a link to your published solutions

Problem 1

Reconsider the US air pollution data set:

library(HSAUR2)
Loading required package: tools
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.2     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.1.0     
── 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 errors
data(USairpollution)

A)

Perform singular value decomposition of this data matrix. Then create the matrix \(D\). Describe what this matrix looks like.

components <- svd(USairpollution)
U <- components$u
D <- components$d
V <- components$v

D_Matrix <-diag(D)

head(D_Matrix)
         [,1]     [,2]    [,3]     [,4]     [,5]    [,6] [,7]
[1,] 7051.949   0.0000   0.000  0.00000  0.00000  0.0000    0
[2,]    0.000 931.1211   0.000  0.00000  0.00000  0.0000    0
[3,]    0.000   0.0000 540.463  0.00000  0.00000  0.0000    0
[4,]    0.000   0.0000   0.000 92.70909  0.00000  0.0000    0
[5,]    0.000   0.0000   0.000  0.00000 85.23724  0.0000    0
[6,]    0.000   0.0000   0.000  0.00000  0.00000 52.9465    0

The D matrix looks like a square pxp matrix that only has values in the diagonal. Everything else contains values of 0.

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 <- (U %*% diag(D) %*% t(V)) %>% round(1)

k <- 7
(Xtilde2 <- U[,1:k] %*% diag(D[1:k]) %*% t(V[,1:k])) %>% round(1)
      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
 [1,]   46 47.6   44  116  8.8 33.4  135
 [2,]   11 56.8   46  244  8.9  7.8   58
 [3,]   24 61.5  368  497  9.1 48.3  115
 [4,]   47 55.0  625  905  9.6 41.3  111
 [5,]   11 47.1  391  463 12.4 36.1  166
 [6,]   31 55.2   35   71  6.5 40.7  148
 [7,]  110 50.6 3344 3369 10.4 34.4  122
 [8,]   23 54.0  462  453  7.1 39.0  132
 [9,]   65 49.7 1007  751 10.9 35.0  155
[10,]   26 51.5  266  540  8.6 37.0  134
[11,]    9 66.2  641  844 10.9 35.9   78
[12,]   17 51.9  454  515  9.0 12.9   86
[13,]   17 49.0  104  201 11.2 30.8  103
[14,]   35 49.9 1064 1513 10.1 31.0  129
[15,]   56 49.1  412  158  9.0 43.4  127
[16,]   10 68.9  721 1233 10.8 48.2  103
[17,]   28 52.3  361  746  9.7 38.7  121
[18,]   14 68.4  136  529  8.8 54.5  116
[19,]   14 54.5  381  507 10.0 37.0   99
[20,]   13 61.0   91  132  8.2 48.5  100
[21,]   30 55.6  291  593  8.3 43.1  123
[22,]   10 61.6  337  624  9.2 49.1  105
[23,]   10 75.5  207  335  9.0 59.8  128
[24,]   16 45.7  569  717 11.8 29.1  123
[25,]   29 43.5  699  744 10.6 25.9  137
[26,]   18 59.4  275  448  7.9 46.0  119
[27,]    9 68.3  204  361  8.4 56.8  113
[28,]   31 59.3   96  308 10.6 44.7  116
[29,]   14 51.5  181  347 10.9 30.2   98
[30,]   69 54.6 1692 1950  9.6 39.9  115
[31,]   10 70.3  213  582  6.0  7.0   36
[32,]   61 50.4  347  520  9.4 36.2  147
[33,]   94 50.0  343  179 10.6 42.7  125
[34,]   26 57.8  197  299  7.6 42.6  115
[35,]   28 51.0  137  176  8.7 15.2   89
[36,]   12 56.7  453  716  8.7 20.7   67
[37,]   29 51.1  379  531  9.4 38.8  164
[38,]   56 55.9  775  622  9.5 35.9  105
[39,]   29 57.3  434  757  9.3 38.9  111
[40,]    8 56.6  125  277 12.7 30.6   82
[41,]   36 54.0   80   80  9.0 40.2  114
plot(X, Xtilde2); 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?

k <- 2
(Xtilde2 <- U[,1:k] %*% diag(D[1:k]) %*% t(V[,1:k])) %>% round(1)
       [,1]  [,2]   [,3]   [,4] [,5] [,6]  [,7]
 [1,]   3.5  23.8   17.2  151.5  3.6 15.2  42.1
 [2,]   5.9  35.8   45.7  244.5  5.5 22.8  63.7
 [3,]  17.3  37.8  352.3  516.8  6.2 24.6  73.1
 [4,]  30.4  63.4  629.6  900.0 10.5 41.3 123.5
 [5,]  17.1  32.9  363.0  498.0  5.5 21.5  64.8
 [6,]   2.4  20.0    0.5  115.4  3.0 12.7  35.1
 [7,] 131.5  80.9 3371.0 3334.4 17.1 56.6 208.0
 [8,]  18.2  20.7  433.8  488.4  3.8 13.8  44.9
 [9,]  34.5  -2.1  962.8  807.5  1.1  0.1  15.2
[10,]  16.2  56.9  260.0  548.2  9.0 36.6 104.9
[11,]  29.3  50.8  643.7  839.4  8.6 33.3 101.6
[12,]  19.2  26.4  443.0  528.5  4.6 17.4  54.9
[13,]   6.3  26.2   86.1  223.8  4.1 16.8  47.7
[14,]  50.9  96.7 1088.6 1481.9 16.2 63.2 190.8
[15,]  11.3 -15.3  362.9  221.1 -1.8 -9.2 -19.7
[16,]  38.8  98.3  746.5 1199.8 16.0 63.6 186.9
[17,]  22.1  73.6  368.8  736.7 11.7 47.4 136.3
[18,]  13.5  71.1  138.3  526.0 11.0 45.4 127.3
[19,]  17.7  35.8  369.9  520.6  6.0 23.4  70.1
[20,]   4.6  18.8   65.6  163.8  2.9 12.1  34.2
[21,]  17.7  60.9  289.1  595.9  9.6 39.2 112.5
[22,]  19.2  58.8  336.5  624.1  9.4 37.9 109.7
[23,]  11.0  36.3  183.1  364.5  5.8 23.4  67.3
[24,]  25.5  43.8  559.9  728.1  7.4 28.7  87.7
[25,]  28.5  32.6  680.9  766.8  5.9 21.8  70.5
[26,]  14.5  42.3  261.2  465.3  6.8 27.3  79.3
[27,]  11.4  38.9  186.7  382.2  6.2 25.0  71.8
[28,]   8.3  44.1   83.6  324.4  6.8 28.2  78.9
[29,]  10.6  37.2  171.2  359.4  5.9 23.9  68.6
[30,]  71.6  81.4 1712.1 1924.8 14.8 54.4 176.3
[31,]  15.9  62.2  233.5  555.9  9.8 39.9 113.6
[32,]  17.4  45.1  331.6  541.4  7.3 29.2  85.6
[33,]  10.4  -4.0  301.8  233.9 -0.2 -2.1  -1.1
[34,]  10.0  30.3  177.0  324.6  4.8 19.5  56.6
[35,]   6.4  17.3  118.4  200.1  2.8 11.2  32.7
[36,]  23.2  53.9  462.5  703.5  8.8 35.0 103.6
[37,]  18.2  44.1  358.0  558.0  7.2 28.6  84.3
[38,]  27.5   4.6  745.0  660.4  1.8  4.0  22.6
[39,]  23.7  65.8  438.8  751.3 10.6 42.5 123.9
[40,]   8.1  33.2  115.3  288.9  5.2 21.3  60.5
[41,]   3.4  13.3   49.9  118.9  2.1  8.5  24.3
cor(as.vector(X), as.vector(Xtilde2))
[1] 0.9965504

The fewest number of dimensions required to have a correlation of at least 0.9 is 2 dimensions.

D)

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)\)
cov_fun <- cov(USairpollution)
head(cov_fun)
               SO2        temp        manu       popul        wind       precip
SO2     550.947561  -73.560671   8527.7201   6711.9945   3.1753049   15.0017988
temp    -73.560671   52.239878   -773.9713   -262.3496  -3.6113537   32.8629884
manu   8527.720122 -773.971341 317502.8902 311718.8140 191.5481098 -215.0199024
popul  6711.994512 -262.349634 311718.8140 335371.8939 175.9300610 -178.0528902
wind      3.175305   -3.611354    191.5481    175.9301   2.0410244   -0.2185311
precip   15.001799   32.862988   -215.0199   -178.0529  -0.2185311  138.5693840
          predays
SO2     229.92988
temp    -82.42616
manu   1968.95976
popul   645.98598
wind      6.21439
precip  154.79290
X <- as.matrix(USairpollution)
X_centered <- scale(X, center = TRUE, scale = FALSE)

n <- nrow(X)

manual_cov <- t(X_centered) %*% X_centered / (n - 1)
head(manual_cov)
               SO2        temp        manu       popul        wind       precip
SO2     550.947561  -73.560671   8527.7201   6711.9945   3.1753049   15.0017988
temp    -73.560671   52.239878   -773.9713   -262.3496  -3.6113537   32.8629884
manu   8527.720122 -773.971341 317502.8902 311718.8140 191.5481098 -215.0199024
popul  6711.994512 -262.349634 311718.8140 335371.8939 175.9300610 -178.0528902
wind      3.175305   -3.611354    191.5481    175.9301   2.0410244   -0.2185311
precip   15.001799   32.862988   -215.0199   -178.0529  -0.2185311  138.5693840
          predays
SO2     229.92988
temp    -82.42616
manu   1968.95976
popul   645.98598
wind      6.21439
precip  154.79290
SVD_US <- svd(X_centered)
eigen_sigma <- eigen(cov_fun)

eigen_sigma$vectors
              [,1]         [,2]         [,3]        [,4]         [,5]
[1,]  0.0168607518 -0.099835625  0.208775573  0.95883106 -0.152191203
[2,] -0.0011417794  0.025814390 -0.071600745 -0.11014784 -0.477854201
[3,]  0.6968327936 -0.710249079 -0.067182201 -0.07319788 -0.009643654
[4,]  0.7170284512  0.692912523  0.056666935  0.04906669  0.010735457
[5,]  0.0004067530 -0.001011680  0.005386606 -0.01506609  0.025401917
[6,] -0.0004336922  0.001225937  0.265807619 -0.16261712 -0.832729325
[7,]  0.0028836950 -0.069155051  0.934279828 -0.18459052  0.232812295
             [,6]          [,7]
[1,] -0.053952911 -2.704138e-02
[2,] -0.852534945 -1.640507e-01
[3,] -0.002153023  1.136208e-03
[4,]  0.002915751 -5.682006e-05
[5,]  0.176541256 -9.838347e-01
[6,]  0.453206155  6.376788e-02
[7,] -0.183568735 -1.891454e-02
SVD_US$v
              [,1]         [,2]         [,3]        [,4]         [,5]
[1,] -0.0168607518  0.099835625  0.208775573 -0.95883106  0.152191203
[2,]  0.0011417794 -0.025814390 -0.071600745  0.11014784  0.477854201
[3,] -0.6968327936  0.710249079 -0.067182201  0.07319788  0.009643654
[4,] -0.7170284512 -0.692912523  0.056666935 -0.04906669 -0.010735457
[5,] -0.0004067530  0.001011680  0.005386606  0.01506609 -0.025401917
[6,]  0.0004336922 -0.001225937  0.265807619  0.16261712  0.832729325
[7,] -0.0028836950  0.069155051  0.934279828  0.18459052 -0.232812295
             [,6]          [,7]
[1,] -0.053952911 -2.704138e-02
[2,] -0.852534945 -1.640507e-01
[3,] -0.002153023  1.136208e-03
[4,]  0.002915751 -5.682006e-05
[5,]  0.176541256 -9.838347e-01
[6,]  0.453206155  6.376788e-02
[7,] -0.183568735 -1.891454e-02
eigen_sigma$values
[1] 6.384720e+05 1.481204e+04 7.019599e+02 2.050019e+02 1.167047e+02
[6] 1.205705e+01 1.448704e+00
SVD_US$d^2/(n-1)
[1] 6.384720e+05 1.481204e+04 7.019599e+02 2.050019e+02 1.167047e+02
[6] 1.205705e+01 1.448704e+00

They both match up!

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.

library(magick)
Linking to ImageMagick 6.9.13.29
Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
Disabled features: fontconfig, x11
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].

recompose_villan <- mysteryU4 %*% diag(mysteryD4) %*% t(mysteryV4)
R <- max(recompose_villan)-min(recompose_villan)
recomposed_scaled_villan <- (recompose_villan-min(recompose_villan))/R
(recomposed_scaled_villan
  %>% as.raster
  %>% image_read
 )

The four dimensional approximation to this image is not enough for me to determine who this well know villain is.

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')

Create both of the images produced by these approximations. At what point can you tell who the mystery villain is?

recompose_villan_10 <- mysteryU10 %*% diag(mysteryD10) %*% t(mysteryV10)

R <- max(recompose_villan_10)-min(recompose_villan_10)
recomposed_scaled_villan_10 <- (recompose_villan_10-min(recompose_villan_10))/R

(recomposed_scaled_villan_10
  %>% as.raster
  %>% image_read
 )

recompose_villan_50 <- mysteryU50 %*% diag(mysteryD50) %*% t(mysteryV50)

R <- max(recompose_villan_50)-min(recompose_villan_50)
recomposed_scaled_villan_50 <- (recompose_villan_50-min(recompose_villan_50))/R

(recomposed_scaled_villan_50
  %>% as.raster
  %>% image_read
 )

I believe that this well know villain is Todd Iverson.

C)

How many numbers need to be stored in memory for each of the following:

  • A full \(700\times 600\) image?
700*600
[1] 420000
  • A 4-dimensional approximation?
(700*4)+4+(4*600)
[1] 5204
  • A 10-dimensional approximation?
(700*10)+10+(10*600)
[1] 13010
  • A 50-dimensional approximation?
(700*50)+50+(50*600)
[1] 65050