Homework01
Kaylie Evans
2025-09-28
Introduction
The following RMD contains CUNY SPS DATA 605 Fall 2025 context for the Homework 01 assignment. This project has the following 4 sections:
1. Geometric Transformation of Shapes Using Matrix Multiplication
2. Matrix Properties and Decomposition
3. Matrix Rank, Properties, and Eigenspace
4. Project: Eigenfaces from the LFW (Labeled Faces in the Wild) Dataset
Data Preparation
Load Libraries
library(tidyverse)
library(jpeg)
library(httr)
library(pracma)
library(ggplot2)
library(gridExtra)
library(reticulate)1. Geometric Transformation of Shapes Using Matrix Multiplication
Task: Create a simple shape (like a square or triangle) using point plots in R. Implement R code to apply different transformations (scaling, rotation, reflection) to the shape by left multiplying a transformation matrix by each of the point vectors. Demonstrate these transformations through animated plots.
Create triangle
I am creating a triangle with a 2×N matrix, where each column is a vertex. The last point equals the first, so the plot closes the shape.
# Define a triangle with 3 points
triangle <- matrix(c(0, 1, 0.5,
0, 0, 1),
nrow = 2, byrow = TRUE)
# Close the triangle
triangle <- cbind(triangle, triangle[,1])
# Plot the original triangle
plot(triangle[1,], triangle[2,], type = "l", lwd = 2,
xlim = c(-2, 2), ylim = c(-2, 2), asp = 1,
main = "Triangle")
Scaling transformation
Here, I multiply the shape by the scaling matrix. The triangle stretches horizontally (×2) and shrinks vertically (×0.5).
# Define a scaling matrix
scale_matrix <- matrix(c(2, 0,
0, 0.5), nrow = 2, byrow = TRUE)
# Apply scaling
triangle_scaled <- scale_matrix %*% triangle
# Plot both original and scaled shapes
plot(triangle[1,], triangle[2,], type = "l", lwd = 2,
xlim = c(-2, 2), ylim = c(-2, 2), asp = 1,
main = "Triangle with Scaling")
lines(triangle_scaled[1,], triangle_scaled[2,], col = "salmon", lwd = 2)
legend("topright", legend = c("Original", "Scaled"),
col = c("black", "salmon"), lty = 1, lwd = 2)
Rotation transformation
I use this formula for the rotation of a 2D point by an angle \(\theta\) around the origin:
\[ R(\theta) = \begin{bmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{bmatrix} \]
# Rotation by 45 degrees
theta <- pi / 4
rotation_matrix <- matrix(c(cos(theta), -sin(theta),
sin(theta), cos(theta)), nrow = 2, byrow = TRUE)
# Apply rotation
triangle_rotated <- rotation_matrix %*% triangle
# Plot
plot(triangle[1,], triangle[2,], type = "l", lwd = 2,
xlim = c(-2, 2), ylim = c(-2, 2), asp = 1,
main = "Triangle with Rotation")
lines(triangle_rotated[1,], triangle_rotated[2,], col = "darkgoldenrod1", lwd = 2)
legend("topright", legend = c("Original", "Rotated"), col = c("black", "darkgoldenrod1"), lty = 1, lwd = 2)
Reflection transformation
The reflection matrix flips the y-coordinates while keeping x the same. This mirrors the triangle across the horizontal axis.
# Reflect across x-axis (flip y values)
reflect_x_matrix <- matrix(c(1, 0,
0, -1), nrow = 2, byrow = TRUE)
# Apply reflection
triangle_reflected_x <- reflect_x_matrix %*% triangle
# Plot
plot(triangle[1,], triangle[2,], type = "l", lwd = 2,
xlim = c(-2, 2), ylim = c(-2, 2), asp = 1,
main = "Triangle Reflected Across X-axis")
lines(triangle_reflected_x[1,], triangle_reflected_x[2,], col = "darkseagreen3", lwd = 2)
legend("topright", legend = c("Original", "Reflected (X-axis)"),
col = c("black", "darkseagreen3"), lty = 1, lwd = 2)
Animation of transformation
I rotate the triangle in small steps (0 → 2π). Each frame shows the shape in a new orientation, creating an animation when run sequentially. Scroll quickly through the below graphs to see the animation.
# Set number of frames
n_frames <- 16 # choose small number like 12, 20, 30
angles <- seq(0, 2*pi, length.out = n_frames)
# Create animation through iterated plots
for (angle in angles) {
R <- matrix(c(cos(angle), -sin(angle),
sin(angle), cos(angle)), nrow = 2, byrow = TRUE)
rotated <- R %*% triangle
# Plot for each triangle
plot(rotated[1,], rotated[2,], type = "l", lwd = 2,
xlim = c(-2, 2), ylim = c(-2, 2), asp = 1,
main = sprintf("Rotation: %.2f rad", angle))
# Pause between each animation
Sys.sleep(0.08)
}
2. Matrix Properties and Decomposition
Task: Write an R function that performs Singular Value Decomposition (SVD) on a grayscale image (which can be represented as a matrix). Use this decomposition to compress the image by keeping only the top k singular values and their corresponding vectors. Demonstrate the effect of different values of k on the compressed image’s quality. You can choose any open-access grayscale image that is appropriate for a professional program.
Proofs + Properties
Prove that AB <> BA
Matrix multiplication is generally not commutative. Computing AB and BA show different results. This numerically demonstrates the property that AB<>BA in most cases. There are special cases where they are equal, like identity matrices.
# Define two matrices A and B
A <- matrix(c(1, 2,
0, 1), nrow = 2, byrow = TRUE)
B <- matrix(c(3, 0,
4, 5), nrow = 2, byrow = TRUE)
# Compute AB and BA
AB <- A %*% B
BA <- B %*% A
# AB = A multiplied by B
print(AB)## [,1] [,2]
## [1,] 11 10
## [2,] 4 5# BA = B multiplied by A
print(BA)## [,1] [,2]
## [1,] 3 6
## [2,] 4 13Prove that A^T A is always symmetric
Part of the defining traits of transposed matrices are symmetry with the original matrix. Here, I get A^T A, transpose it, and compare that to the original. They are symmetrical.
# Define a random matrix
set.seed(64)
A <- matrix(rnorm(9), nrow = 3)
# Compute A^T A
ATA <- t(A) %*% A
# Check symmetry
all.equal(ATA, t(ATA))## [1] TRUEProve that det(A^T A) is greater than or equal to 0
The matrix A^T A is positive semidefinite. Its determinant is the product of its eigenvalues, and since eigenvalues are non-negative, the determinant cannot be negative. The code computes it for an example and results in positive 26.9, which is greater than 0. Sometimes it’s exactly 0 if A is rank-deficient.
# Generate a random matrix
set.seed(64)
A <- matrix(rnorm(9), nrow = 3)
# Compute determinant of A^T A
det_ATA <- det(t(A) %*% A)
det_ATA## [1] 26.93264Decomposition
Read an image
I have a picture of my cat, Astro, saved in GitHub. I download the image, make sure it is fully grayscale, and show the image, which is stored as a matrix of pixel intensities between 0 and 1.
# URL to raw image file on my github
url <- "https://raw.githubusercontent.com/evanskaylie/DATA605/main/IMG_2111.jpg"
# Download and read the image
resp <- httr::GET(url)
if (resp$status_code != 200) stop("Could not fetch image")
img_bin <- httr::content(resp, "raw")
# Write to a temp file
tmpfile <- tempfile(fileext = ".jpg")
writeBin(img_bin, tmpfile)
# Read the image with jpeg
astro_img <- readJPEG(tmpfile)
# Convert to grayscale if there is any RGB there
if (length(dim(astro_img)) == 3) {
gray_astro <- 0.299*astro_img[,,1] + 0.587*astro_img[,,2] + 0.114*astro_img[,,3]
} else {
gray_astro <- astro_img
}
# Display the grayscale image
image(1:ncol(gray_astro), 1:nrow(gray_astro), t(gray_astro[nrow(gray_astro):1, ]),
col = gray((0:255)/255),
xlab = "", ylab = "", main = "Small box. Chunky cat.",
axes = FALSE, asp = 1)
Perform SVD
SVD breaks a matrix into orthogonal directions and importance weights. Keeping the largest singular values allows me to compress the image while preserving the main structure.
# Singular Value Decomposition
svd_img <- svd(gray_astro)
U <- svd_img$u
D <- diag(svd_img$d)
V <- svd_img$vCompress the image
I am keeping only the first k singular values and vectors. This approximates the image with reduced complexity. Then, I visualize the image, starting with the original and then getting more blurry.
# Function to reconstruct with top k singular values
reconstruct_astro <- function(U, D, V, k) {
U_k <- U[, 1:k]
D_k <- D[1:k, 1:k]
V_k <- V[, 1:k]
U_k %*% D_k %*% t(V_k)
}
# Visualize compression
par(mfrow = c(2,2))
# Original
image(1:ncol(gray_astro), 1:nrow(gray_astro), t(gray_astro[nrow(gray_astro):1, ]),
col = gray((0:255)/255),
main = "Original",
axes = FALSE, asp = 1)
# Compressed with different k
for (k in c(50, 20, 5)) {
astro_k <- reconstruct_astro(U, D, V, k)
image(t(apply(astro_k, 2, rev)),
col = gray((0:255)/255),
main = paste("k =", k),
axes = FALSE, asp = 1)
}
3. Matrix Rank, Properties, and Eigenspace
Task: This prompt includes determining the rank of a given matrix, matrix rank boundaries, rank & row reduction, computing the eigenvalues & eigenvectors, and diagonalization of a matrix.
Determine the rank of the given matrix
If rank < number of rows, some rows are linear combinations of others. If rank < number of columns, some columns are linear combinations of others. Rank gives the dimension of the row space or column space, which is the maximum number of linearly independent rows/columns.
# Define the given matrix A
A <- matrix(c(2, 4, 1, 3,
-2, -3, 4, 1,
5, 6, 2, 8,
-1, -2, 3, 7), nrow = 4, byrow = TRUE)
# Find the rank with qr decomposition
rank_A <- qr(A)$rank
# Find linear dependent status for rows and columms
cat("The rank of A is:", rank_A, "\n")## The rank of A is: 4if (rank_A < nrow(A)) {
cat("There are linear dependencies among the rows.\n")
} else {
cat("All rows are linearly independent.\n")
}## All rows are linearly independent.if (rank_A < ncol(A)) {
cat("There are linear dependencies among the columns.\n")
} else {
cat("All columns are linearly independent.\n")
}## All columns are linearly independent.Matrix rank boundaries
I generate a random matrix with m > n dimensions, 7x3. qr(A)$rank gives the actual rank of the matrix, which confirms the rank cannot exceed n (3) and must be at least 1 if the matrix is non-zero.
# Example: 7 x 3 matrix (m > n)
set.seed(64)
A <- matrix(rnorm(21), nrow = 7, ncol = 3)
# Get the rank
rank_A <- qr(A)$rank
rank_A## [1] 3# Rank ranges
cat("Maximum possible rank:", min(nrow(A), ncol(A)), "\n")## Maximum possible rank: 3cat("Minimum possible rank (for non-zero matrices): 1 \n")## Minimum possible rank (for non-zero matrices): 1The matrix below has numbers chosen for this example. Rows 1 and 2 are dependent (row 2 = 3 * row 1). Row 3 is independent. Again, qr(A)$rank gives the actual rank of the matrix. Then, when I check the rows/columns, the number of linearly independent rows = rank (row space dimension) and the number of linearly independent columns = rank (column space dimension).
# Example matrix
B <- matrix(c(1, 2, 3,
3, 6, 9,
0, 1, 1), nrow = 3, byrow = TRUE)
# Get the rank
rank_B <- qr(B)$rank
rank_B## [1] 2# Rank dimensions for rows and columns are equal
cat("Rank equals dimension of row space:", rank_B, "\n")## Rank equals dimension of row space: 2cat("Rank equals dimension of column space:", rank_B, "\n")## Rank equals dimension of column space: 2Rank and row reduction
I get qr(B)$rank which returns the number of linearly independent rows/columns. This gives a quick numeric rank.
# Define the given matrix B
B <- matrix(c(2, 5, 7,
4, 10, 14,
1, 2.5, 3.5), nrow = 3, byrow = TRUE)
# Rank using QR decomposition
rank_B <- qr(B)$rank
rank_B## [1] 1Here, rref(B) converts B into reduced row echelon form. Non-zero rows correspond to independent rows. There is only 1 non-zero row, so rank = 1 makes sense.
# Row reduce
B_ref <- rref(B)
B_ref## [,1] [,2] [,3]
## [1,] 1 2.5 3.5
## [2,] 0 0.0 0.0
## [3,] 0 0.0 0.0Special properties: Matrix B is a rank-deficient matrix because its rank is 1, which is less than the number of rows or columns. This indicates that most of the rows and columns are linearly dependent: specifically, the second row is exactly twice the first row, and the third row is half of the first row. Similarly, the columns are also dependent, with the second column being 2.5 times the first column and the third column being 3.5 times the first column. Because the rank is less than the matrix dimension, the determinant ofB is zero, making it a singular matrix that cannot be inverted. Essentially, all the information in B is contained in a single independent row (or column), which also means the matrix has a low-rank structure and could be represented as an outer product of a row vector with appropriate scaling factors. These properties highlight that B is highly redundant, with multiple rows and columns providing no new independent information.
Compute the eigenvalues and eigenvectors
# Define the given matrix A
A <- matrix(c(3, 1, 2,
0, 5, 4,
0, 0, 2), nrow = 3, byrow = TRUE)
A## [,1] [,2] [,3]
## [1,] 3 1 2
## [2,] 0 5 4
## [3,] 0 0 2The characteristic polynomial is obtained from det(A - λI) = 0.
For this matrix:
A -λI =
\[\begin{bmatrix} 3-\lambda & 1 & 2 \\ 0 & 5-\lambda & 4 \\ 0 & 0 & 2-\lambda \end{bmatrix}\]Because A - λI is upper triangular, the determinant is the product of the diagonal entries:
det(A - λI) = (3 - λ)(5 - λ)(2 - λ)
Then, I set this equal to 0. The eigenvalues are then λ1 = 3, λ2 = 5 λ3 = 2.
The next step is to compute the eigenvectors and eigenvalues
# Compute eigenvalues and eigenvectors
eig <- eigen(A)
# Eigenvalues
eig$values## [1] 5 3 2# Eigenvectors
eig$vectors## [,1] [,2] [,3]
## [1,] 0.4472136 1 -0.3713907
## [2,] 0.8944272 0 -0.7427814
## [3,] 0.0000000 0 0.5570860The columns of eig$vectors correspond to eigenvectors of the eigenvalues in eig$values.
For clarity, we can label them:
λ = 3 → eigenvector: eig$vectors[, which(eig$values==3)]
λ = 5 → eigenvector: eig$vectors[, which(eig$values==5)]
λ = 2 → eigenvector: eig$vectors[, which(eig$values==2)]
Next, I verify linear independence. The determinant of the eigenvector matrix is non-zero, so the eigenvectors are linearly independent.
# Verify linear independence
det(eig$vectors)## [1] -0.4982729We know the fact that all eigenvalues are distinct, which guarantees that their corresponding eigenvectors are linearly independent.
Diagonalization of matrix
I am using the matrix A defined above. The determinant of the eigenvector matrix is non-zero so eigenvectors are linearly independent. All eigenvalues are distinct, which guarantees diagonalizability. Therefore, matrix A can be diagonalized!
# Check if it can be diagonalized
eig <- eigen(A)
eig$values## [1] 5 3 2eig$vectors## [,1] [,2] [,3]
## [1,] 0.4472136 1 -0.3713907
## [2,] 0.8944272 0 -0.7427814
## [3,] 0.0000000 0 0.5570860# Check linear independence by determinant
det(eig$vectors)## [1] -0.4982729The next step is to construct the diagonal matrix D and matrix P. D = diag(λ1, λ2, λ3). P is a matrix with the eigenvectors of A as its columns. P^-1 A P = D is the diagonalization equation.
# Diagonal matrix of eigenvalues
D <- diag(eig$values)
D## [,1] [,2] [,3]
## [1,] 5 0 0
## [2,] 0 3 0
## [3,] 0 0 2# Matrix of eigenvectors
P <- eig$vectors
P## [,1] [,2] [,3]
## [1,] 0.4472136 1 -0.3713907
## [2,] 0.8944272 0 -0.7427814
## [3,] 0.0000000 0 0.5570860# Verify diagonalization: P^-1 A P = D
round(solve(P) %*% A %*% P, 6)## [,1] [,2] [,3]
## [1,] 5 0 0
## [2,] 0 3 0
## [3,] 0 0 2The eigenvectors of matrix A show directions in space that do not change direction when multiplied by A. The eigenvalues tell us how much vectors in those directions are stretched or shrunk. For this matrix, the eigenvalues 3, 5, and 2 mean that vectors along these three directions are stretched by those amounts. Because all the eigenvalues are positive, the matrix does not rotate the vectors. It only stretches them along the eigenvector directions.
4. Project: Eigenfaces from the LFW (Labeled Faces in the Wild) Dataset
Task: Using the LFW (Labeled Faces in the Wild) dataset, build and visualize eigenfaces that account for 80% of the variability in the dataset. The LFW dataset is a well-known dataset containing thousands of labeled facial images, available for academic research.
Download the LFW dataset and preprocess images
I am using a Python emulator to get the LFW dataset into my RMD. For the PCA, I want to make sure each person we look at has enough data, so I am taking only people that have 20+ faces in their system. The images are grayscale and converted into an R array. Each image is then flattened into a 1D vector, and all vectors are stacked into a matrix X_flat where each row represents an image and each column corresponds to a pixel. This format is ready for PCA analysis, and the output confirms the number of images, their dimensions, and the size of the flattened matrix.
# Using my system Python environment, import Python libraries
sklearn <- import("sklearn")
np <- import("numpy")
plt <- import("matplotlib.pyplot")
# Fetch LFW dataset using sklearn
lfw <- sklearn$datasets$fetch_lfw_people(
min_faces_per_person = as.integer(20), # force integer
resize = 0.5,
color = FALSE
)
# Convert images to R array
X_r <- py_to_r(lfw$images) # convert to R array
n_samples <- dim(X_r)[1]
h <- dim(X_r)[2]
w <- dim(X_r)[3]
# Flatten images
X_flat <- matrix(X_r, nrow = n_samples, ncol = h*w)
# Results
cat("Number of images:", n_samples, "\nImage size:", h, "x", w, "\nFlattened size:", dim(X_flat)[2], "\n")## Number of images: 3023
## Image size: 62 x 47
## Flattened size: 2914Apply PCA
I apply Principal Component Analysis (PCA) to the flattened image matrix X_flat to identify the main patterns of variation across all faces. PCA finds a set of orthogonal vectors, called principal components or eigenfaces, that capture the directions of greatest variability in the dataset. I compute how many components are needed to explain 80% of the total variance, which allows me to reduce dimensionality while retaining most of the information.
# Using my system Python environment, use sklearn
PCA <- sklearn$decomposition$PCA
# Keep components explaining 80% variance
pca <- PCA(n_components = 0.80, svd_solver = "full")
X_pca <- pca$fit_transform(X_flat)
# Results
cat("Number of principal components to explain 80% variance:", pca$n_components_, "\n")## Number of principal components to explain 80% variance: 39Visualize Eigenfaces
The first few eigenfaces are visualized to show there is a general face shape for each of the displayed images, including chins/necks, eyes, nose, mouth, and some cheek shapes.
# Convert to R array if not already
components_r <- py_to_r(pca$components_)
# Number of principal components
n_pc <- nrow(components_r)
# Reshape each row into h x w matrix
eigenfaces_r <- array(components_r, dim = c(n_pc, h, w))
# Function to plot first few eigenfaces
par(mfrow = c(2,5), mar = c(1,1,2,1))
for (i in 1:10) {
image(t(apply(eigenfaces_r[i,,], 2, rev)), col = gray((0:255)/255), axes = FALSE,
main = paste("PC", i))
}
Here, I reconstruct some images using only these principal components and compare them to the original images. This demonstrates that most facial features can be preserved using a reduced set of eigenfaces, illustrating PCA’s effectiveness in compressing and summarizing complex image data. Faces are always a bit unsettling when generated like this. Naive generating systems tend to give an uncanny valley look at best, and not even close to the original at worse. Either way, here you go!
# Reconstruct first few images using top PCs
n_images <- 5
# R indexes from 1 to n_images
X_reconstructed <- pca$inverse_transform(X_pca[1:n_images, , drop = FALSE])
# Convert R integers to Python integers using as.integer()
py2int <- function(x) as.integer(x)
par(mfrow = c(2, n_images), mar = c(1,1,2,1))
for (i in 1:n_images) {
# Original
image(t(apply(matrix(X_flat[i,], nrow=h, ncol=w), 2, rev)),
col = gray((0:255)/255), axes = FALSE, main = "Original")
# Reconstructed
image(t(apply(matrix(X_reconstructed[i,], nrow=h, ncol=w), 2, rev)),
col = gray((0:255)/255), axes = FALSE, main = "Reconstructed")
}