Eigen Analysis

Example from Image Compression

## library need to open a jpg file
library(jpeg)
## Download from the author's website
myurl <- "http://polytopes.net/Tora_Sleeping.JPG"
z <- tempfile()
download.file(myurl,z,mode="wb")
Kitty <- readJPEG(z)
## dimension
d <- dim(Kitty)
r <- Kitty[,,1]
g <- Kitty[,,2]
b <- Kitty[,,3]
library(mvtnorm)
library(ggplot2)
Warning: package 'ggplot2' was built under R version 4.2.3
library(matlib)
## Standard deviation
sigma <- matrix(c(4,3,3,4), ncol = 2, nrow = 2) 
## Mean
mu <- c(5, 5)
n <- 1000
set.seed(123)
x <- rmvnorm(n = n, mean = mu, sigma = sigma)
d <- data.frame(x)
p2 <- ggplot(d, aes(x = X1, y = X2)) + geom_point(alpha = .5) + geom_density_2d()
p2

A <- matrix(c(2, 1, 1, 2), 2, 2)
eigen(A)
eigen() decomposition
$values
[1] 3 1

$vectors
          [,1]       [,2]
[1,] 0.7071068 -0.7071068
[2,] 0.7071068  0.7071068
### Practical Applications
y <- t(x) %*% x
eigen(y)
eigen() decomposition
$values
[1] 58081.296  1069.938

$vectors
          [,1]       [,2]
[1,] 0.7067305 -0.7074829
[2,] 0.7074829  0.7067305
prcomp(d, center = FALSE)
Standard deviations (1, .., p=2):
[1] 7.624922 1.034896

Rotation (n x k) = (2 x 2):
          PC1        PC2
X1 -0.7067305 -0.7074829
X2 -0.7074829  0.7067305

Diagonalization

## Example 150
Q <- matrix(c(1, 1, -1, 1), 2, 2)
D <- matrix(c(3, 0, 0, 1), 2, 2)
Q %*% D %*% inv(Q)
     [,1] [,2]
[1,]    2    1
[2,]    1    2
## Example 151
A <- matrix(c(2, -1, 1, -1, -3, 1, 3, 3, 3), 3, 3)
eigen(A)
eigen() decomposition
$values
[1]  4.397382 -3.824458  1.427076

$vectors
          [,1]       [,2]       [,3]
[1,] 0.7437078  0.2542122  0.8674476
[2,] 0.1624484  0.9508909 -0.3981322
[3,] 0.6484668 -0.1765859 -0.2983713
## Example 152
A <- matrix(c(1, 0, 3, -2, 2, -3, 1, -3, -2), 3, 3)
eigen(A)
eigen() decomposition
$values
[1]  4.50349224 -3.56576496  0.06227272

$vectors
           [,1]        [,2]       [,3]
[1,]  0.5276592 -0.01504018 -0.7715961
[2,] -0.6521968 -0.47442002 -0.5343413
[3,]  0.5442565 -0.88017012 -0.3451359
### Practical Applications
A <- matrix(c(1/4, 1/5, 1/3, 1/6, 1/4, 1/5, 1/3, 1/3, 1/4, 2/5, 1/6, 1/3, 1/4, 1/5, 1/6, 1/6), 4, 4)
p <- eigen(A)$vectors
D <- diag(eigen(A)$values)
p %*% D %*% solve(p)
             [,1]         [,2]         [,3]         [,4]
[1,] 0.2500000+0i 0.2500000+0i 0.2500000+0i 0.2500000+0i
[2,] 0.2000000-0i 0.2000000+0i 0.4000000+0i 0.2000000+0i
[3,] 0.3333333+0i 0.3333333-0i 0.1666667-0i 0.1666667+0i
[4,] 0.1666667+0i 0.3333333+0i 0.3333333+0i 0.1666667+0i
ini <- c(1/4, 1/4, 1/4, 1/4)
ini %*% p %*% D^(10000) %*% solve(p)
             [,1]        [,2]         [,3]         [,4]
[1,] 0.2435175+0i 0.276212+0i 0.2841037+0i 0.1961669+0i

PCA and eigen analysis (Principal component analysis - Wikipedia)

Tora Data

## library need to open a jpg file
library(jpeg)
## Download from the author's website
myurl <- "http://polytopes.net/Tora_Sleeping.JPG"
z <- tempfile()
download.file(myurl,z,mode="wb")
Kitty <- readJPEG(z)
## dimension
d <- dim(Kitty)
r <- Kitty[,,1]
g <- Kitty[,,2]
b <- Kitty[,,3]

kitty.r.pca <- prcomp(r, center = FALSE)
kitty.g.pca <- prcomp(g, center = FALSE)
kitty.b.pca <- prcomp(b, center = FALSE)

kitty.rgb.pca <- list(kitty.r.pca, kitty.g.pca, kitty.b.pca)

num <- c(2, 50, 100, 1000, 2000)

for (i in num) {
     pca.img <- sapply(kitty.rgb.pca, function(j) {
         compressed.img <- j$x[,1:i] %*% t(j$rotation[,1:i])
     }, simplify = 'array')
     writeJPEG(pca.img, paste('tora_compressed_', 
     round(i,0), '_components.jpg', sep = ''))
 }

Image Compression

library(jpeg)
path = "C:/Users/Subhajit_code/My Drive/MTH209A/Lectures/L6"
setwd(path)

tour = readJPEG(paste(path,"paris_tourist.jpg", sep = "/"))
no_tour = readJPEG(paste(path,"paris_no_tourist.jpg", sep = "/"))
new = readJPEG(paste(path,"rome.jpg", sep = "/"))
dim(tour); dim(no_tour); dim(new)
[1] 459 715   3
[1] 458 715   3
[1] 500 750   3
rno_tour = no_tour[, ,1]
gno_tour = no_tour[, ,2]
bno_tour = no_tour[, ,3]
no_tour_rgb = list(rno_tour,gno_tour,bno_tour)

prtour = prcomp(rno_tour, center = FALSE)
pgtour = prcomp(gno_tour, center = FALSE)
pbtour = prcomp(bno_tour, center = FALSE)

tour.pca = list(prtour, pgtour, pbtour)

library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
# Create a dataframe for easier plotting
df = data.frame(scheme = rep(c("R", "G", "B"), each = 458), index = rep(1:458,  3), var = c(prtour$sdev^2, pgtour$sdev^2, pbtour$sdev^2))

df %<>% group_by(scheme) %>%
  mutate(propvar = 100*var/sum(var)) %>%
  mutate(cumsum = cumsum(propvar)) %>%
  ungroup()

# scree plot
library(ggplot2)
#relevel to make it look nicer
#df$scheme = factor(df$scheme,levels(df$scheme)[c(3,2,1)])

df %>% ggplot(aes(x = index, y = propvar, fill = scheme)) + 
    geom_bar(stat="identity") + 
    labs(title="Screeplot of Principal Component", x ="Principal Component", 
         y="% of Variance") + geom_line()  + 
    scale_x_continuous(limits = c(0,30)) +
    facet_wrap(~scheme)
Warning: Removed 1284 rows containing missing values (`position_stack()`).
Warning: Removed 3 rows containing missing values (`geom_bar()`).
Warning: Removed 1284 rows containing missing values (`geom_line()`).

df %>% ggplot(aes( x = index, y = cumsum, fill = scheme)) + 
    geom_bar(stat="identity") + 
    labs(title="Cumulative Proportion of Variance Explained Principal Component", x="Principal Component", 
         y="Cumulative % of Variance") + geom_line() + 
    scale_x_continuous(limits = c(0,30)) +
    facet_wrap(~scheme)
Warning: Removed 1284 rows containing missing values (`position_stack()`).
Warning: Removed 3 rows containing missing values (`geom_bar()`).
Warning: Removed 1284 rows containing missing values (`geom_line()`).

################################################################
# This is the number of desired PCs
#pcnum = c(2,30,200, 300)
pcnum = (1:5)*10

compress <- function(trained_rgb_pca, newrgb, pcnum, dims){
   r_rotate = trained_rgb_pca[[1]]$rotation
   g_rotate = trained_rgb_pca[[2]]$rotation
   b_rotate = trained_rgb_pca[[3]]$rotation
   
   r = newrgb[[1]]
   g = newrgb[[2]]
   b = newrgb[[3]]
   
  pred_r = (r %*% r_rotate)[,1:pcnum] %*% t(r_rotate[,1:pcnum])
  pred_g = (g %*% g_rotate)[,1:pcnum] %*% t(g_rotate[,1:pcnum])
  pred_b = (b %*% b_rotate)[,1:pcnum] %*% t(b_rotate[,1:pcnum])
  
  pred.pca = list(pred_r, pred_g, pred_b)
  pred.array = array(as.numeric(unlist(pred.pca)),dim = dims)
  return(pred.array)
}

for(i in pcnum){
    pca.img = compress(tour.pca, no_tour_rgb, pcnum = i, dims = c( 458, 715,3))
     writeJPEG(pca.img, paste('no_tourist_compressed_', 
     round(i,0), '_components.jpg', sep = ''))
}

Exercises

  1. Is the following matrix invertible?

  2. Diagonalize the matrix given above.

  3. Create compressed images of the “paris_tourist.jpg” and “rome.jpg” datasets.