Linear Algebra & Applications

2201NSC

Week 10 - Singular Value Decomposition

Problem 1

Use the SVD to give a basis and dimension for the column, row, and null spaces of:

(a) $\,A = \begin{pmatrix} 1 & 3 & 2 \\ 2 & 1 & 4 \\ 4 & 7 & 8 \end{pmatrix}.$ (b) $\,A = \begin{pmatrix} 1 & 0 & 2 & 1 \\ 0 & 1 & 1 & 4 \\ 0 & 2 & -1 & -1 \end{pmatrix}.$






Problem 1 (cont.)

(a) $\,A = \begin{pmatrix} 1 & 3 & 2 \\ 2 & 1 & 4 \\ 4 & 7 & 8 \end{pmatrix}.$

If \(A=U\Sigma V^{T}\) is a SVD with singular values \(\sigma_1\ge\sigma_2\ge\cdots\ge 0\), then

  • $\mathrm{Col}(A)=\operatorname{span}\{u_1,\dots,u_r\}$ where $u_i$ are the columns of $U$ corresponding to nonzero singular values (an orthonormal basis).
  • $\mathrm{Row}(A)=\operatorname{span}\{v_1,\dots,v_r\}$ where $v_i$ are the columns of $V$ corresponding to nonzero singular values (an orthonormal basis).
  • $\mathrm{Null}(A)=\operatorname{span}\{v_{r+1},\dots,v_n\}$ (right singular vectors for zero singular values).

Here \(r=\operatorname{rank}(A)\).


Problem 2 (cont.) Python code
						
import numpy as np

# Define the matrix
A = np.array([[1, 3, 2],
              [2, 1, 4],
              [4, 7, 8]], dtype=float)

# Compute the full SVD
U, s, Vt = np.linalg.svd(A)

# Round results to 6 decimal places
U = np.round(U, 6)
s = np.round(s, 6)
Vt = np.round(Vt, 6)

print("Matrix A:")
print(A)

print("\nLeft singular vectors (U):")
print(U)

print("\nSingular values (s):")
print(s)

print("\nRight singular vectors (V^T):")
print(Vt)

Problem 2 (cont.) Python code
						
# Rank of A = number of nonzero singular values
tol = 1e-10
rank = np.sum(s > tol)
print("\nRank of A:", rank)

# Column space basis (first 'rank' columns of U)
col_space = U[:, :rank]
print("\nBasis for column space:")
print(np.round(col_space, 6))
print("Dimension of column space:", col_space.shape[1])

# Row space basis (first 'rank' rows of V^T, equivalently first 'rank' cols of V)
row_space = Vt[:rank, :].T
print("\nBasis for row space:")
print(np.round(row_space, 6))
print("Dimension of row space:", row_space.shape[1])

# Null space basis (last n-rank columns of V)
null_space = Vt[rank:, :].T
print("\nBasis for null space:")
print(np.round(null_space, 6))
print("Dimension of null space:", null_space.shape[1])

Run Python online here: trinket.io/python3 

						
import numpy as np
# Define the matrix
A = np.array([ [1, 3, 2], [2, 1, 4], [4, 7, 8] ], dtype=float)
U, s, Vt = np.linalg.svd(A) # Compute the full SVD
# Round results to 6 decimal places
U = np.round(U, 6)
s = np.round(s, 6)
Vt = np.round(Vt, 6)
print("Matrix A:")
print(A)
print("\nLeft singular vectors (U):")
print(U)
print("\nSingular values (s):")
print(s)
print("\nRight singular vectors (V^T):")
print(Vt)
# Rank of A = number of nonzero singular values
tol = 1e-10
rank = np.sum(s > tol)
print("\nRank of A:", rank)
# Column space basis (first 'rank' columns of U)
col_space = U[:, :rank]
print("\nBasis for column space:")
print(np.round(col_space, 6))
print("Dimension of column space:", col_space.shape[1])
# Row space basis (first 'rank' rows of V^T, equivalently first 'rank' cols of V)
row_space = Vt[:rank, :].T
print("\nBasis for row space:")
print(np.round(row_space, 6))
print("Dimension of row space:", row_space.shape[1])
# Null space basis (last n-rank columns of V)
null_space = Vt[rank:, :].T
print("\nBasis for null space:")
print(np.round(null_space, 6))
print("Dimension of null space:", null_space.shape[1])

Problem 2

Download the Octave code $\texttt{magic_approx_SVD.m}$ from Learning@Griffith.

  1. Calculate the rank 1, 2, and 3 approximations to the $10\times 10$ magic matrix, $A$.
  2. For each of the approximations, calculate the level of compression achieved (assume that the elements of the magic square are independent, otherwise we're going to have trouble).
  3. Prove, by hand, that there is no $2\times 2$ magic square.



Problem 2 (cont.) Octave code

						
pkg load linear-algebra  % Load the linear algebra package

A = magic(10);           % Generate a 10x10 magic square matrix
r = 1;                   % Set the rank for the approximation

% Compute the rank-r SVD approximation
[u, s, v] = svds(A, r);  % Compute the truncated SVD 
A_approx = u * s * v';   % Reconstruct the rank-r

% Find the largest & smallest elements across both A and A_approx
amax = max(max(max(A)), max(max(A_approx)));
amin = min(min(min(A)), min(min(A_approx)));

% Rescale matrices to the range [0, 63] for grayscale image display
A_color = (A - amin) / (amax - amin) * 63;
A_approx_color = (A_approx - amin) / (amax - amin) * 63;

Problem 2 (cont.) Octave code

						
% Display the original matrix as a grayscale image
figure(1); 
colormap(bone);       % Set the colormap to 'bone'
image(A_color);       % Display the rescaled matrix as an image
title('Original A');  % Add a title

% Display the rank-r approximation as a grayscale image
figure(2); 
colormap(bone);
image(A_approx_color);
title('Rank-1 Approximation');

Run Octave online here: octave-online.net

👉 Compute
compression:
 [m,n]=size(A)
$1-r\left(\dfrac{1-m+n}{mn}\right)$ 

						
pkg load linear-algebra  % Load the linear algebra package

A = magic(10);           % Generate a 10x10 magic square matrix
r = 1;                   % Set the rank for the approximation

% Compute the rank-r SVD approximation
[u, s, v] = svds(A, r);  % Compute the truncated SVD 
A_approx = u * s * v';   % Reconstruct the rank-r

% Find the largest & smallest elements across both A and A_approx
amax = max(max(max(A)), max(max(A_approx)));
amin = min(min(min(A)), min(min(A_approx)));

% Rescale matrices to the range [0, 63] for grayscale image display
A_color = (A - amin) / (amax - amin) * 63;
A_approx_color = (A_approx - amin) / (amax - amin) * 63;

% Display the original matrix as a grayscale image
figure(1); 
colormap(bone);       % Set the colormap to 'bone'
image(A_color);       % Display the rescaled matrix as an image
title('Original A');  % Add a title

% Display the rank-r approximation as a grayscale image
figure(2); 
colormap(bone);
image(A_approx_color);
title('Rank-1 Approximation');

Problem 3

Download the Octave code from $\texttt{SVD_image_compress.m}$ from Learning@Griffith.

  1. Calculate the rank $10$, $50$, and $100$ approximations to the image $\texttt{Lambs.jpg}.$
  2. For each of the approximations, calculate the level of compression achieved.
  3. Try compressing $\texttt{note_screenshot.png}.$ What sort of artefacts do you encounter?



Problem 3 (cont.) Python code

						
import numpy as np
from scipy.sparse.linalg import svds  # Truncated SVD

import matplotlib.pyplot as plt # Library for plotting
from PIL import Image  # For manipulating and saving images  
import requests        # Access images from the web
from io import BytesIO # Treat raw byte data

# Rank for approximation
r = 50

# Download the image from URL
url = "https://esc-maths.github.io/assets/imgs/Lambs.jpg"
response = requests.get(url)
a = np.array(Image.open(BytesIO(response.content))) # numpy array

# Extract color channels as float
R = a[:, :, 0].astype(float)
G = a[:, :, 1].astype(float)
B = a[:, :, 2].astype(float)

Note: To access the other image use the link: https://esc-maths.github.io/assets/imgs/note.png


Problem 3 (cont.) Python code

						
# Function for rank-r SVD approximation
def svd_approx(X, r):
    U, S, Vt = svds(X, k=r)
    # Reverse order: svds returns smallest → largest
    S = np.diag(S[::-1])
    # Now reorder the corresponding singular vectors
    U = U[:, ::-1]
    Vt = Vt[::-1, :]
    return U @ S @ Vt

# Apply to each channel
R_approx = svd_approx(R, r)
G_approx = svd_approx(G, r)
B_approx = svd_approx(B, r)

# Reconstruct image
a_out = np.stack([R_approx, G_approx, B_approx], axis=2)
a_out = np.clip(a_out, 0, 255).astype(np.uint8)

Problem 3 (cont.)

						
# Plot side by side (like montage in Octave)
fig, axes = plt.subplots(1, 2, figsize=(15, 10))

axes[0].imshow(a)
axes[0].set_title("Original")
axes[0].axis("off")

axes[1].imshow(a_out)
axes[1].set_title(f"Rank-{r} Approximation")
axes[1].axis("off")

plt.show()

Run Python online here: trinket.io/python3

👉 Compute
compression:
 m, n, _ = a.shape
$1-r\left(\dfrac{1-m+n}{mn}\right)$ 

						
import numpy as np
from scipy.sparse.linalg import svds
import matplotlib.pyplot as plt
from PIL import Image
import requests
from io import BytesIO
r = 50 # Rank
url = "https://esc-maths.github.io/assets/imgs/Lambs.jpg"
response = requests.get(url)
a = np.array(Image.open(BytesIO(response.content))) # numpy array
R = a[:, :, 0].astype(float)
G = a[:, :, 1].astype(float)
B = a[:, :, 2].astype(float)
# Function for rank-r SVD approximation
def svd_approx(X, r):
    U, S, Vt = svds(X, k=r)
    # Reverse order: svds returns smallest → largest
    S = np.diag(S[::-1])
    U = U[:, ::-1] # Now reorder the corresponding singular vectors
    Vt = Vt[::-1, :]
    return U @ S @ Vt
# Apply to each channel
R_approx = svd_approx(R, r)
G_approx = svd_approx(G, r)
B_approx = svd_approx(B, r)
# Reconstruct image
a_out = np.stack([R_approx, G_approx, B_approx], axis=2)
a_out = np.clip(a_out, 0, 255).astype(np.uint8)
# Plot side by side (like montage in Octave)
fig, axes = plt.subplots(1, 2, figsize=(15, 10))
axes[0].imshow(a)
axes[0].set_title("Original")
axes[0].axis("off")
axes[1].imshow(a_out)
axes[1].set_title(f"Rank-{r} Approximation")
axes[1].axis("off")
plt.show()