library(jpeg)
library(ggplot2)

Introduction

This project demonstrates dimensionality reduction using Principal Component Analysis (PCA) in a practical setting: image compression. An image can be represented as a matrix of pixel intensities. PCA/SVD allows us to approximate this matrix using a smaller number of components.

Load and preprocess the image

img <- readJPEG("redukcja.jpg")   
dim(img)
## [1] 1688 3000    3
is_rgb <- (length(dim(img)) == 3 && dim(img)[3] == 3)
if (!is_rgb) {
  img <- array(rep(img, 3), dim = c(dim(img), 3))
}
R <- img[,,1]
G <- img[,,2]
B <- img[,,3]

summary(as.vector(R))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.01176 0.05098 0.13481 0.17255 1.00000
summary(as.vector(G))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.04706 0.11765 0.17634 0.23137 1.00000
summary(as.vector(B))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.00000 0.07059 0.13725 0.17944 0.23137 1.00000
par(mar=c(1,1,2,1))
plot.new()
rasterImage(img, 0, 0, 1, 1)
title("Original image (RGB)")

The image is treated as a 3-channel RGB array. To perform dimensionality reduction, we apply SVD-based PCA separately to each color channel (R, G, B) and then combine the reconstructed channels back into a color image.

PCA/SVD reconstruction

We use Singular Value Decomposition (SVD), which is closely related to PCA for matrix data.

reconstruct_svd <- function(X, k) {
  sv <- svd(X)
  U <- sv$u[, 1:k, drop = FALSE]
  D <- diag(sv$d[1:k], nrow = k, ncol = k)
  V <- sv$v[, 1:k, drop = FALSE]
  U %*% D %*% t(V)
}

mse <- function(a, b) mean((a - b)^2)
psnr <- function(a, b, maxval = 1) {
  err <- mse(a, b)
  if (err == 0) return(Inf)
  10 * log10(maxval^2 / err)
}
reconstruct_rgb <- function(img_rgb, k) {
  Rr <- reconstruct_svd(img_rgb[,,1], k)
  Gr <- reconstruct_svd(img_rgb[,,2], k)
  Br <- reconstruct_svd(img_rgb[,,3], k)

  rec <- array(0, dim = dim(img_rgb))
  rec[,,1] <- Rr
  rec[,,2] <- Gr
  rec[,,3] <- Br
  rec[rec < 0] <- 0
  rec[rec > 1] <- 1
  rec
}

Reconstructions for different k

We reconstruct the image for several values of k and evaluate reconstruction quality using MSE and PSNR.

kmax <- min(dim(img)[1], dim(img)[2])
ks <- c(5, 10, 25, 50, 100, 150)
ks <- ks[ks <= kmax]

metrics_rgb <- data.frame(k = ks, MSE = NA_real_, PSNR = NA_real_)

par(mfrow=c(2,3), mar=c(1,1,2,1))
plot.new(); rasterImage(img, 0, 0, 1, 1); title("Original (RGB)")

for (i in seq_along(ks)) {
  k <- ks[i]
  rec <- reconstruct_rgb(img, k)

  # MSE averaged across channels
  mse_total <- mean(c(mse(img[,,1], rec[,,1]),
                      mse(img[,,2], rec[,,2]),
                      mse(img[,,3], rec[,,3])))
  metrics_rgb$MSE[i] <- mse_total
  metrics_rgb$PSNR[i] <- 10 * log10(1 / mse_total)

  plot.new()
  rasterImage(rec, 0, 0, 1, 1)
  title(paste("Reconstruction k =", k))
}

par(mfrow=c(1,1))

metrics_rgb
##     k         MSE     PSNR
## 1   5 0.011054724 19.56452
## 2  10 0.008059600 20.93687
## 3  25 0.005344079 22.72127
## 4  50 0.003812189 24.18826
## 5 100 0.002507637 26.00735
## 6 150 0.001813946 27.41376

Quality vs number of components

ggplot(metrics_rgb, aes(x = k, y = MSE)) +
  geom_line() + geom_point() +
  labs(title = "RGB reconstruction error (MSE) vs number of components",
       x = "Number of components (k)",
       y = "MSE")

ggplot(metrics_rgb, aes(x = k, y = PSNR)) +
  geom_line() + geom_point() +
  labs(title = "RGB reconstruction quality (PSNR) vs number of components",
       x = "Number of components (k)",
       y = "PSNR (dB)")

Compression ratio

n <- dim(img)[1]
m <- dim(img)[2]

metrics_rgb$original_values <- 3 * n * m
metrics_rgb$stored_values_k <- 3 * (n*metrics_rgb$k + metrics_rgb$k + m*metrics_rgb$k)
metrics_rgb$compression_ratio <- metrics_rgb$stored_values_k / metrics_rgb$original_values

metrics_rgb
##     k         MSE     PSNR original_values stored_values_k compression_ratio
## 1   5 0.011054724 19.56452        15192000           70335       0.004629739
## 2  10 0.008059600 20.93687        15192000          140670       0.009259479
## 3  25 0.005344079 22.72127        15192000          351675       0.023148697
## 4  50 0.003812189 24.18826        15192000          703350       0.046297393
## 5 100 0.002507637 26.00735        15192000         1406700       0.092594787
## 6 150 0.001813946 27.41376        15192000         2110050       0.138892180

The results illustrate a clear quality–compression trade-off. For example, increasing k from 5 to 50 improves PSNR from 19.56 dB to 24.19 dB, while the storage ratio increases from 0.46% to 4.63%. Further increasing k to 150 improves PSNR to 27.41 dB, but the storage ratio rises to 13.89%, indicating diminishing returns.

psnr_min <- metrics_rgb$PSNR[1]
psnr_max <- max(metrics_rgb$PSNR)
target <- psnr_min + 0.95 * (psnr_max - psnr_min)

k_star <- metrics_rgb$k[min(which(metrics_rgb$PSNR >= target))]
target
## [1] 27.02129
k_star
## [1] 150

Using a simple heuristic (achieving 95% of the PSNR improvement), the recommended k is the smallest value that meets the target PSNR, which balances quality and compression.

ggplot(metrics_rgb, aes(x = k, y = compression_ratio)) +
  geom_line() +
  geom_point() +
  labs(title = "Approximate storage ratio vs number of components (RGB)",
       x = "Number of components (k)",
       y = "Stored values / original values")

Conclusion

PCA/SVD provides a compact representation of the image. Increasing k improves reconstruction quality (lower MSE, higher PSNR), but with diminishing returns. This illustrates the key idea of dimensionality reduction: representing high-dimensional data using fewer components while preserving as much information as possible.

References

Image source: https://wallpapercave.com/wp/wp12680196.jpg