Library

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.0.4     
## ── 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
library(ggplot2)
library(Matrix)
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(pracma)
## 
## Attaching package: 'pracma'
## 
## The following objects are masked from 'package:Matrix':
## 
##     expm, lu, tril, triu
## 
## The following object is masked from 'package:purrr':
## 
##     cross
library(tibble)
library(imager)
## Loading required package: magrittr
## 
## Attaching package: 'magrittr'
## 
## The following objects are masked from 'package:pracma':
## 
##     and, mod, or
## 
## The following object is masked from 'package:purrr':
## 
##     set_names
## 
## The following object is masked from 'package:tidyr':
## 
##     extract
## 
## 
## Attaching package: 'imager'
## 
## The following object is masked from 'package:magrittr':
## 
##     add
## 
## The following object is masked from 'package:stringr':
## 
##     boundary
## 
## The following object is masked from 'package:dplyr':
## 
##     where
## 
## The following object is masked from 'package:tidyr':
## 
##     fill
## 
## The following objects are masked from 'package:stats':
## 
##     convolve, spectrum
## 
## The following object is masked from 'package:graphics':
## 
##     frame
## 
## The following object is masked from 'package:base':
## 
##     save.image

1. Geometric Transformation of Shapes Using Matrix Multiplication

# Build a square centered at the origin
square <- tibble(
  x = c(-1, 1, 1, -1, -1),
  y = c(-1, -1, 1, 1, -1)
)

plot_shape <- function(df, title = "Shape") {
  ggplot(df, aes(x, y)) +
    geom_path(linewidth = 1) +
    geom_point(size = 1.5) +
    coord_fixed(xlim = c(-4, 4), ylim = c(-4, 4)) +
    theme_minimal(base_size = 12) +
    labs(title = title, x = "x", y = "y")
}
plot_shape(square, "Original square")

# Define transformations
scaling <- matrix(c(1.6, 0,
              0, 0.7), nrow = 2, byrow = TRUE)

# Rotation by theta = 30 degrees
theta <- pi/6
rotation <- matrix(c(cos(theta), -sin(theta),
                     sin(theta), cos(theta)), nrow = 2, byrow = TRUE)

# Reflection across the x-axis
reflection <- matrix(c(1, 0, 
                       0, -1), nrow = 2, byrow = TRUE)

# Apply M to every (x, y)
apply_transform <- function(df, M) {
  pts <- as.matrix(df[, c("x", "y")])
  tr <- t(M %*% t(pts))
  tibble(x = tr[,1], y = tr[,2])
}

sq_scaled <- apply_transform(square, scaling)
sq_rotated <- apply_transform(square, rotation)
sq_reflected <- apply_transform(square, reflection)

# Compare with original 
p_scale <- plot_shape(
  bind_rows(
    square %>% mutate(what = "original"),
    sq_scaled %>% mutate(what = "scaled (1.6x, 0.7y)")
  ),
  "Scaling vs Original"
) + aes(color = what) + scale_color_brewer(palette = "Set1") + theme(legend.position = "bottom")

p_rot <- plot_shape(
  bind_rows(
    square %>% mutate(what = "original"),
    sq_rotated %>% mutate(what = "rotate (30 degree)")
  ),
  "Rotation vs Original"
) + aes(color = what) + scale_color_brewer(palette = "Set1") + theme(legend.position = "bottom")

p_refl <- plot_shape(
  bind_rows(
    square %>% mutate(what = "original"),
    sq_reflected %>% mutate(what = "reflected across x-axis")
  ),
  "Reflection vs Original"
) + aes(color = what) + scale_color_brewer(palette = "Set1") + theme(legend.position = "bottom")

p_scale; p_rot; p_refl

# Animation Frames

angles <- seq(0, pi/2, length.out = 6)
scales <- seq(0.6, 1.6, length.out = 6)

frames <- map2(angles, scales, function(ang, sc) {
  Rf <- matrix(c(cos(ang), -sin(ang),
                 sin(ang), cos(ang)), 2, 2, byrow = TRUE)
  Sf <- matrix(c(sc, 0,
               0, sc), 2, 2, byrow = TRUE)
  
  M <- Rf %*% Sf
  deg <- round(ang * 180 / pi)
  apply_transform(square, M) |>
    dplyr::mutate(
      # plotmath: theta==72*degree ~ "," ~ s==1.4
      frame = sprintf("theta==%d*degree~\",\"~s==%.1f", deg, sc)
    )
})

frames_df <- dplyr::bind_rows(frames)
ggplot() +
  geom_path(data = square, aes(x, y), color = "grey65") +
  geom_path(data = frames_df, aes(x, y, group = frame)) +
  facet_wrap(~ frame, ncol = 2) +
  coord_fixed(xlim = c(-4,4), ylim = c(-4,4)) +
  theme_minimal(base_size = 11) +
  labs(title = "Static Fames of Rotation + Scaling", x = "x", y = "y")

2. Matrix Properties and Decomposition

# Prove that AB not equal to BA
mA <- matrix(c(0,1,
              0,0),2,2, byrow = TRUE)
mB <- matrix(c(0,0,
              1,0),2,2, byrow = TRUE)

mAB <- mA %*% mB
mBA <- mB %*% mA
mAB; mBA
##      [,1] [,2]
## [1,]    1    0
## [2,]    0    0
##      [,1] [,2]
## [1,]    0    0
## [2,]    0    1
# Prove A^T A is always symmetric
set.seed(1)
A <- matrix(rnorm(9), 3, 3)
AtA <- t(A) %*% A
identical(AtA, t(AtA))
## [1] TRUE
max(abs(AtA - t(AtA)))
## [1] 0
# Prove that the determinant of A^T A is non-negative 
AtA_vals <- eigen(AtA, symmetric = TRUE)$values
AtA_det <- det(AtA)
AtA_vals
## [1] 3.5239523 1.5898832 0.4515081
AtA_det
## [1] 2.529652
# Singular Value Decomposition and Image Compression 
img_col <- load.image("border collie.jpg")
img_gray <- grayscale(img_col)
par(mfrow = c(1, 2))
plot(img_col, main = "Original Colored Image")
plot(img_gray, main = "Original Grayscale Image")

a_img <- as.array(img_gray)[,,1,1]
a_img <- a_img / max(a_img)
dim(a_img)
## [1] 736 490
m <- nrow(a_img); n <- ncol(a_img)

# SVD
sv_img <- svd(a_img)
energy <- cumsum(sv_img$d^2) / sum(sv_img$d^2)
rankmax <- min(m, n)

# Reconstruct using top k singular values/vectors
svd_reconstruct_k <- function(k) {
  k <- max(1, min(k, length(sv_img$d)))        # clamp k to valid range
  U <- sv_img$u[, 1:k, drop = FALSE]
  D <- diag(sv_img$d[1:k], k)
  V <- sv_img$v[, 1:k, drop = FALSE]
  U %*% D %*% t(V)
}

k_vals <- c(5, 20, 50, n)

par(mfrow = c(2, 3), mar = c(3,3,2,1))
plot(img_gray, main = "Original Grayscale (imager)")

image(a_img[, ncol(a_img):1], col = gray.colors(256),
      main = "Original Matrix view", asp = m/n)

for (k in k_vals) {
  Rk <- svd_reconstruct_k(k)
  image(Rk[, ncol(Rk):1], col = gray.colors(256),
        main = paste0("Reconstruction, k=", k,
                      "\nVar kept ~ ", sprintf("%.1f%%", 100*energy[k])),
        asp = m/n)
}

par(mfrow = c(1,1))

# Metrics table 
fro <- function(M) sqrt(sum(M^2))
metrics <- data.frame(
  k = k_vals,
  retained_variance = energy[k_vals],
  frobenius_error   = sapply(k_vals, function(k) fro(a_img - svd_reconstruct_k(k))),
  storage_original  = m*n,
  storage_rank_k    = m*k_vals + k_vals + n*k_vals, 
  compression_ratio = (m*n) / (m*k_vals + k_vals + n*k_vals)
)
metrics

3. Matrix Rank, Properties, and Eigenspace

# Matrix A
A_rank <- matrix(c(
  2,  4, 1, 3,
 -2, -3, 4, 1,
  5,  6, 2, 8,
 -1, -2, 3, 7
), nrow = 4, byrow = TRUE)

A_rank
##      [,1] [,2] [,3] [,4]
## [1,]    2    4    1    3
## [2,]   -2   -3    4    1
## [3,]    5    6    2    8
## [4,]   -1   -2    3    7
qr(A_rank)$rank      
## [1] 4
Matrix::rankMatrix(A_rank)
## [1] 4
## attr(,"method")
## [1] "tolNorm2"
## attr(,"useGrad")
## [1] FALSE
## attr(,"tol")
## [1] 8.881784e-16

The rank is 4, because rank equals the maximum number of linearly independent rows/columns, rank = 4 for a 4x4 matrix means all rows and columns are linearly independent. Therefore, no nontrivial linear dependencies exist among the rows or among the columns

Rank boundaries

Problem a: Given a non-zero m x n matrix with m > n, determine the maximum and minimum possible rank, assuming that the matrix is non-zero

m <- 6; n <- 3
c(min_rank = 1, max_rank = n)
## min_rank max_rank 
##        1        3

problem b: Prove that the rank of a matrix equals the dimension of its row space (or column space). Provide an example to illustrate the concept.

Row operations covert A to row-echelon form without changing the row space dimension. The number of pivot rows in that echelon form equals the rank, so rank A = dim(row space). Since rank A = rank(A⊤), the column space of A has the same dimension.

C_ex <- rbind(c(1,2,3),
              c(2,4,6),
              c(1,1,1))
pracma::rref(C_ex)   
##      [,1] [,2] [,3]
## [1,]    1    0   -1
## [2,]    0    1    2
## [3,]    0    0    0
qr(C_ex)$rank 
## [1] 2

Rank and Row Reduction for B

B_rank <- rbind(
  c(2,   5,  7),
  c(4,  10, 14),
  c(1, 2.5, 3.5)
)

qr(B_rank)$rank
## [1] 1
pracma::rref(B_rank)
##      [,1] [,2] [,3]
## [1,]    1  2.5  3.5
## [2,]    0  0.0  0.0
## [3,]    0  0.0  0.0

Each row is a scalar multiple of the first row, so all rows lie on the same 1-dimensional line in row space.

Eigenvalues and Eigenvectors

A_tri <- rbind(
  c(3, 1, 2),
  c(0, 5, 4),
  c(0, 0, 2)
)

eig <- eigen(A_tri)
eig$values      
## [1] 5 3 2
eig$vectors     
##           [,1] [,2]       [,3]
## [1,] 0.4472136    1 -0.3713907
## [2,] 0.8944272    0 -0.7427814
## [3,] 0.0000000    0  0.5570860
qr(eig$vectors)$rank  
## [1] 3

Eigenvalues are 3, 5, 2. They are distinct, therefore the corresponding eigenvectors are linearly independent.

Diagonalization

P <- eig$vectors

D <- diag(eig$values)

norm(A_tri %*% P - P %*% D, type = "F")
## [1] 0
A_rebuilt <- P %*% D %*% solve(P)

max(abs(A_tri - A_rebuilt))  
## [1] 8.881784e-16

Eigenvectors mark directions in R3 that the transformation doesn’t rotate, only re-scale. For A tri, the scaling factors are 3, 5, and 2, so the transformation stretches or compresses space by those amounts along the corresponding eigenvector axes.

4. Project: Eaigenfaces from the LFW (Labeled Faces in the Wild) Dataset

Context: Eigenfaces are a popular application of Principal Component Analysis (PCA) in computer vision. They are used for face recognition by finding the principal components (eigenvectors) of the covariance matrix of a set of facial images. These principal components represent the “eigenfaces” that can be combined to approximate any face in the data set.

Task: Using the LFW (Labeled Faces in the Wild) data set, build and visualize eigenfaces that account for 80% of the variability in the data set. The LFW data set is a well-known data set containing thousands of labeled facial images, available for academic research.