Question 1: Geometric Transformation of Shapes Using Matrix Multiplication

1. Geometric Transformation of Shapes Using Matrix Multiplication

Description of the Code and what it does. In the Code below I create a rectangle using a matrix and transforming it into a shape. I then conclude various transformations to showcase what coyld be done using R. The transformation included Enlarging the Shape 1.5 times, rotating the rectangle 90 degrees and reflecting it on the x-axis.

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.3.3
library(gganimate)
## Warning: package 'gganimate' was built under R version 4.3.3
library(transformr)
## Warning: package 'transformr' was built under R version 4.3.3
rectangle <- matrix(c(0, 0, 4, 4, 0,  
                      0, 2, 2, 0, 0), 
                    nrow=2, byrow=TRUE)

Write a function to apply a transformation from a matrix to shape.This function multiplies a transformation matrix with the shape matrix to get the transformed coordinates.

# Function to apply a transformation matrix
apply_transformation <- function(shape, matrix) {
  return(matrix %*% shape)
}
# Scaling transformation (enlarging by a factor of 1.5)
scaling_matrix <- matrix(c(1.5, 0, 
                           0, 1.5), 
                         nrow=2, byrow=TRUE)

scaled_rectangle <- apply_transformation(rectangle, scaling_matrix)

to_df <- function(shape, label) {
  data.frame(x = shape[1, ], y = shape[2, ], type = label)
}

df_original <- to_df(rectangle, "Original")
df_scaled <- to_df(scaled_rectangle, "Scaled")
df_original_scaled <- rbind(df_original, df_scaled)
ggplot(df_original_scaled, aes(x, y, color = type, group = type)) +
  geom_polygon(fill = NA, size = 1) +
  coord_fixed() +
  theme_minimal() +
  labs(title = "Geometric Transformations of a Square")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

I’ve added a loop to animate the scaling transformation using gganimate. The shape will smoothly scale from 1x to 2x over 10 steps, and the animation will be saved as “scaling.gif”. Its important to note that PDF did not support the gganimate package so please to see the animation visit the project Rpubs link from here

scaling_animation <- function(shape, steps = 10) {
  scale_factors <- seq(1, 1.5, length.out = steps)  
  all_shapes <- do.call(rbind, lapply(scale_factors, function(factor) {
    scale_matrix <- matrix(c(factor, 0, 
                              0, factor), 
                            nrow=2, byrow=TRUE)
    shape_transformed <- apply_transformation(shape, scale_matrix)
    to_df(shape_transformed, paste("Scale Factor:", round(factor, 2)))
  }))
  return(all_shapes)
}

df_animation <- scaling_animation(rectangle)

p <- ggplot(df_animation, aes(x, y, group = type, color = type)) +
  geom_polygon(fill = NA, size = 1) +
  coord_fixed() +
  theme_minimal() +
  labs(title = "Scaling Animation: {frame_time}") +
  transition_states(type, transition_length = 2, state_length = 1)

animate(p, nframes = 10, renderer = gifski_renderer("scaling.gif"))

# Rotation Animation (90 degrees counterclockwise)
rotation_animation <- function(shape) {
  rotation_matrix <- matrix(c(0, -1, 
                              1, 0), #we can change the rotaion matrix to switch directions of the rotation. 
                            nrow=2, byrow=TRUE)
  rotated_shape <- apply_transformation(shape, rotation_matrix)
  df_rotated <- to_df(rotated_shape, "Rotated 90 Degrees")
  return(df_rotated)
}

df_rotated <- rotation_animation(rectangle)

ggplot(rbind(to_df(rectangle, "Original"), df_rotated), aes(x, y, color = type, group = type)) +
  geom_polygon(fill = NA, size = 1) +
  coord_fixed() +
  theme_minimal() +
  labs(title = "Before and After Rotation (90 Degrees)")

# Rotation Animation with Incremental Angle Change
rotation_animation <- function(shape, steps = 10) {
  angles <- seq(0, pi/2, length.out = steps)  # Rotate from 0 to 90 degrees
  all_shapes <- do.call(rbind, lapply(angles, function(angle) {
    rotation_matrix <- matrix(c(cos(angle), -sin(angle), 
                                sin(angle), cos(angle)), 
                              nrow=2, byrow=TRUE)
    shape_transformed <- apply_transformation(shape, rotation_matrix)
    to_df(shape_transformed, paste("Rotation Angle:", round(angle * 180 / pi, 2)))
  }))
  return(all_shapes)
}

df_rotated_animation <- rotation_animation(rectangle)

p_rot <- ggplot(df_rotated_animation, aes(x, y, group = type, color = type)) +
  geom_polygon(fill = NA, size = 1) +
  coord_fixed() +
  theme_minimal() +
  labs(title = "Rotation Animation: {frame_time}") +
  transition_states(type, transition_length = 2, state_length = 1)

animate(p_rot, nframes = 10, renderer = gifski_renderer("rotation.gif"))

# Reflection across Y-axis Animation
reflection_animation <- function(shape) {
  reflection_matrix <- matrix(c(-1, 0, 
                                 0, 1), 
                               nrow=2, byrow=TRUE)
  reflected_shape <- apply_transformation(shape, reflection_matrix)
  df_reflected <- to_df(reflected_shape, "Reflected Across Y-axis")
  return(df_reflected)
}

df_reflected <- reflection_animation(rectangle)

ggplot(rbind(to_df(rectangle, "Original"), df_reflected), aes(x, y, color = type, group = type)) +
  geom_polygon(fill = NA, size = 1) +
  coord_fixed() +
  theme_minimal() +
  labs(title = "Before and After Reflection (Y-axis)")

reflection_animation <- function(shape, steps = 10) {
  reflection_factors <- seq(1, -1, length.out = steps)  # Flip gradually
  all_shapes <- do.call(rbind, lapply(reflection_factors, function(factor) {
    reflection_matrix <- matrix(c(factor, 0, 
                                   0, 1), 
                                 nrow=2, byrow=TRUE)
    shape_transformed <- apply_transformation(shape, reflection_matrix)
    to_df(shape_transformed, paste("Reflection Factor:", round(factor, 2)))
  }))
  return(all_shapes)
}

df_reflected_animation <- reflection_animation(rectangle)

p_reflect <- ggplot(df_reflected_animation, aes(x, y, group = type, color = type)) +
  geom_polygon(fill = NA, size = 1) +
  coord_fixed() +
  theme_minimal() +
  labs(title = "Reflection Animation: {frame_time}") +
  transition_states(type, transition_length = 2, state_length = 1)

animate(p_reflect, nframes = 10, renderer = gifski_renderer("reflection.gif"))

Question 2: Matrix Properties and Decomposition

In the following script, I aimed to demonstrate how Singular Value Decomposition (SVD) for image compression. It begins by reading a grayscale image and converting it into a numerical matrix. The function svd_image_compression() applies SVD to decompose the image matrix into its singular values. The image is then reconstructed using only the top k singular values and their corresponding vectors, effectively compressing the image while retaining important features. The effect of different values of k (5, 20, 50, 100) on image quality is visualized, where lower k values result in greater compression but reduced clarity. This method is useful for reducing image storage while maintaining a balance between quality and file size.

# Function to perform SVD-based image compression
svd_image_compression <- function(image_matrix, k) {
  # Perform Singular Value Decomposition
  svd_result <- svd(image_matrix)
  
  # Keep only the top k singular values and corresponding vectors
  compressed_matrix <- svd_result$u[, 1:k] %*% diag(svd_result$d[1:k]) %*% t(svd_result$v[, 1:k])
  
  return(compressed_matrix)
}

library(jpeg)
library(ggplot2)

# Read image and convert to grayscale matrix
image <- readJPEG("C:/Users/PC/Documents/grayscale.jpg")
gray_image <- 0.2989 * image[,,1] + 0.5870 * image[,,2] + 0.1140 * image[,,3]  # Convert to grayscale

# Apply SVD compression for different values of k
k_values <- c(15, 75)
compressed_images <- lapply(k_values, function(k) svd_image_compression(gray_image, k))

# Plot original and compressed images
par(mfrow = c(1, length(k_values) + 1))
image(gray_image, col = gray.colors(255), main = "Original Image")
for (i in 1:length(k_values)) {
  image(compressed_images[[i]], col = gray.colors(255), main = paste("Compressed Image (k=", k_values[i], ")", sep=""))
}

Question 3: Matrix Rank, Properties, and Eigenspace

  1. Determine the Rank of the Given Matrix:

The rank of a matrix represents the number of linearly independent rows or columns. A full-rank matrix means all rows are linearly independent, whereas a lower rank indicates dependencies. In this case, the computed rank tells us whether there are linear dependencies among the rows or columns of A. Since the matrix A has 4 rows and 4 columns, a rank of 4 means all rows and columns are linearly independent. we can conclude that A is invertible (non-singular).

# Calculate the Rank of a Given Matrix
A <- matrix(c(2, 4, 1, 3,
              -2, -3, 4, 1,
              5, 6, 2, 8,
              -1, -2, 3, 7), nrow=4, byrow=TRUE)

# Compute the rank for A
rank_A <- qr(A)$rank
print(paste("The rank of matrix A is:", rank_A))
## [1] "The rank of matrix A is: 4"
  1. Matrix Rank Boundaries:
# Rank Boundaries Analysis
m <- 5  # Number of rows
n <- 3  # Number of columns
B <- matrix(runif(m * n), nrow = m, ncol = n)  # Random m x n matrix

# Maximum and Minimum Rank
max_rank <- min(m, n)  # Maximum rank is the smaller of m and n
min_rank <- 1  # Minimum rank for a non-zero matrix is at least 1
rank_B <- qr(B)$rank
print(paste("Maximum possible rank:", max_rank))
## [1] "Maximum possible rank: 3"
print(paste("Minimum possible rank:", min_rank))
## [1] "Minimum possible rank: 1"
print(paste("Computed rank of B:", rank_B))
## [1] "Computed rank of B: 3"
# Demonstrate that the rank equals the dimension of row space
C <- matrix(c(1, 2, 3,
              4, 5, 6,
              7, 8, 9), nrow=3, byrow=TRUE)  # Example of dependent rows
rank_C <- qr(C)$rank
print(paste("The rank of matrix C is:", rank_C))
## [1] "The rank of matrix C is: 2"

min(m,n)=3.Given an m x n matrix where m > n, the rank is at most n since there cannot be more than n linearly independent columns. The minimum rank is 1 for a non-zero matrix. This script demonstrates these boundaries by computing the rank of a random m x n matrix. It shows that the rank corresponds to the dimension of the row space, meaning the number of linearly independent rows.

  1. Rank and Row Reduction:
B_matrix <- matrix(c(2, 5, 7,
                     4, 10, 14,
                     1, 2.5, 3.5), nrow=3, byrow=TRUE)
rank_B_matrix <- qr(B_matrix)$rank
rank_B_matrix
## [1] 1

Matrix B has dependent rows since the second row is exactly 2 times the first row, and the third row is a scaled version of the first row. This means B is rank-deficient. The computed rank shows how many linearly independent rows exist in B. Performing row reduction helps reveal these dependencies by reducing the matrix to its row echelon form.

  1. Compute the Eigenvalues and Eigenvectors:
# Define matrix A
A <- matrix(c(3, 1, 2,
              0, 5, 4,
              0, 0, 2), nrow=3, byrow=TRUE)

# Compute eigenvalues and eigenvectors
eigen_result <- eigen(A)

# Extract eigenvalues and eigenvectors
eigenvalues <- eigen_result$values
eigenvectors <- eigen_result$vectors

# Print the results
print("Eigenvalues of A:")
## [1] "Eigenvalues of A:"
print(eigenvalues)
## [1] 5 3 2
print("Eigenvectors of A:")
## [1] "Eigenvectors of A:"
print(eigenvectors)
##           [,1] [,2]       [,3]
## [1,] 0.4472136    1 -0.3713907
## [2,] 0.8944272    0 -0.7427814
## [3,] 0.0000000    0  0.5570860
# Verify if eigenvectors are linearly independent
rank_eigenvectors <- qr(eigenvectors)$rank
print(paste("Rank of eigenvectors matrix:", rank_eigenvectors))
## [1] "Rank of eigenvectors matrix: 3"

The eigenvalues are the solutions to det(A - λI) = 0, forming the characteristic polynomial. The eigenvectors correspond to these eigenvalues and define the directions of transformation. If the rank of the eigenvector matrix equals the number of columns (full rank), the eigenvectors are linearly independent. If not, it means some eigenvectors are dependent, usually occurring when the matrix is defective.

# Determine if matrix A can be diagonalized

# Define matrix A
A <- matrix(c(3, 1, 2,
              0, 5, 4,
              0, 0, 2), nrow=3, byrow=TRUE)

# Compute eigenvalues and eigenvectors
eigen_result <- eigen(A)

# Extract eigenvalues and eigenvectors
eigenvalues <- eigen_result$values
eigenvectors <- eigen_result$vectors

# Check if the matrix is diagonalizable
# A matrix is diagonalizable if the eigenvectors form a basis (i.e., they are linearly independent)
rank_eigenvectors <- qr(eigenvectors)$rank
is_diagonalizable <- (rank_eigenvectors == ncol(A))

# Print results
print("Eigenvalues of A:")
## [1] "Eigenvalues of A:"
print(eigenvalues)
## [1] 5 3 2
print("Eigenvectors of A:")
## [1] "Eigenvectors of A:"
print(eigenvectors)
##           [,1] [,2]       [,3]
## [1,] 0.4472136    1 -0.3713907
## [2,] 0.8944272    0 -0.7427814
## [3,] 0.0000000    0  0.5570860
print(paste("Is matrix A diagonalizable?", is_diagonalizable))
## [1] "Is matrix A diagonalizable? TRUE"
# If diagonalizable, construct the diagonal matrix D and matrix P of eigenvectors
if (is_diagonalizable) {
  D <- diag(eigenvalues)  # Diagonal matrix of eigenvalues
  P <- eigenvectors       # Matrix of eigenvectors
  P_inv <- solve(P)       # Inverse of P
  A_reconstructed <- P %*% D %*% P_inv  # Verify diagonalization
  
  print("Diagonal matrix D:")
  print(D)
  
  print("Matrix P (eigenvectors):")
  print(P)
  
  print("Reconstructed A (should be close to original A):")
  print(A_reconstructed)
}
## [1] "Diagonal matrix D:"
##      [,1] [,2] [,3]
## [1,]    5    0    0
## [2,]    0    3    0
## [3,]    0    0    2
## [1] "Matrix P (eigenvectors):"
##           [,1] [,2]       [,3]
## [1,] 0.4472136    1 -0.3713907
## [2,] 0.8944272    0 -0.7427814
## [3,] 0.0000000    0  0.5570860
## [1] "Reconstructed A (should be close to original A):"
##      [,1] [,2] [,3]
## [1,]    3    1    2
## [2,]    0    5    4
## [3,]    0    0    2

A matrix is diagonalizable if it has a full set of linearly independent eigenvectors. The diagonal matrix D contains the eigenvalues, and P contains the eigenvectors. The transformation P * D * P^(-1) reconstructs the original matrix A. Geometric Interpretation: Eigenvectors represent directions in which A stretches or shrinks vectors in R^3. Eigenvalues tell how much the corresponding eigenvector is scaled. If all eigenvalues are distinct, A has a full basis of eigenvectors and is diagonalizable. If some eigenvalues are repeated but eigenvectors are still independent, A is still diagonalizable. If A lacks enough independent eigenvectors, it is not diagonalizable and requires Jordan form instead.

Question 4: Project: Eigenfaces from the LFW (Labeled Faces in the Wild) Dataset

In R, I used the reticulate package to interface with Python and access the LFW (Labeled Faces in the Wild) dataset from sklearn.datasets. First, I imported the sklearn.datasets module using import("sklearn.datasets"), allowing us to call Python functions within R. I then used fetch_lfw_people() to download the dataset, setting min_faces_per_person = 70 to filter for individuals with at least 70 images and resizing the images to 40% of their original size (resize = 0.25) to reduce computational complexity. Finally, the image data, originally in NumPy format, was converted into an R data frame using as.data.frame(lfw_people$data), making it compatible with R functions for further processing, such as Principal Component Analysis (PCA).

I first begin by loading the LFW dataset, specifically using the deep-funneled images, which are pre-aligned for better recognition accuracy. We ensure that only individuals with at least 70 images are included to maintain a high-quality dataset.

library(reticulate)
## Warning: package 'reticulate' was built under R version 4.3.3
library(ggplot2)
library(jpeg)
library(stats)

# Import Python's sklearn module
sklearn <- import("sklearn.datasets")

# Load LFW dataset with deep-funneled images
lfw_people <- sklearn$fetch_lfw_people(min_faces_per_person = as.integer(70), resize = 0.4)


# Convert to a dataframe
lfw_data <- as.data.frame(lfw_people$data)
# Extract image data and shape
image_data <- lfw_people$images  # Corrected extraction of images
image_shape <- dim(image_data)  # Ensure correct dimensions
print(image_shape)  # Verify if the shape is correct
## [1] 1288   50   37
# Perform PCA using prcomp in R
flattened_data <- matrix(image_data, nrow=image_shape[1], ncol=image_shape[2] * image_shape[3])
pca_result <- prcomp(flattened_data, center=TRUE, scale=FALSE)

# Compute the number of principal components required to explain 80% variance
explained_variance <- cumsum(pca_result$sdev^2) / sum(pca_result$sdev^2)
k_80 <- which(explained_variance >= 0.80)[1]
print(paste("Number of components required to retain 80% variability:", k_80))
## [1] "Number of components required to retain 80% variability: 34"
# Function to reconstruct images using the top k eigenfaces
reconstruct_images <- function(images, pca, k) {
  reconstructed_images <- lapply(1:nrow(images), function(i) {
    image <- as.numeric(images[i,])  # Ensure numeric format
    image_pca <- predict(pca, newdata=as.data.frame(matrix(image, nrow=1)))[,1:k]
    reconstructed <- image_pca %*% t(pca$rotation[,1:k])
    return(matrix(reconstructed, nrow=image_shape[2], ncol=image_shape[3]))  # Reshape correctly
  })
  return(reconstructed_images)
}

# Select multiple test images and reconstruct them
num_images <- 5  # Number of images to reconstruct
original_images <- lapply(1:num_images, function(i) {
  matrix(flattened_data[i,], nrow=image_shape[2], ncol=image_shape[3])
})
reconstructed_images <- reconstruct_images(flattened_data[1:num_images,], pca_result, k_80)

# Compute Mean Squared Error (MSE) for multiple images
mse_values <- sapply(1:num_images, function(i) mean((original_images[[i]] - reconstructed_images[[i]])^2))
# Improve image clarity by increasing resolution and contrast
normalize_image <- function(img) {
  img <- img - min(img)  # Shift values to start at 0
  img <- img / max(img)  # Scale between 0 and 1
  return(img)
}

par(mfrow=c(num_images,2), mar=c(2,2,2,2))
for (i in 1:num_images) {
  orig <- normalize_image(original_images[[i]])
  recon <- normalize_image(reconstructed_images[[i]])
  image(t(apply(orig, 2, rev)), col=gray.colors(256, start=0, end=1), main=paste("Original Image", i))
  image(t(apply(recon, 2, rev)), col=gray.colors(256, start=0, end=1), main=paste("Reconstructed Image", i))
}

num_eigenfaces <- 10 # Display first 10 eigenfaces
par(mfrow=c(2,5), mar=c(2,2,2,2))
for (i in 1:num_eigenfaces) {
  eigenface <- matrix(pca_result$rotation[,i], nrow=image_shape[2], ncol=image_shape[3])
  image(t(apply(eigenface, 2, rev)), col=gray.colors(256, start=0, end=1), main=paste("Eigenface", i))
}