MATH307Su26 - Low Rank Approximations of Images

# Low Rank Approximations of Images > [!abstract] The project in one line > Work the singular value decomposition by hand on small matrices, confirm every number in your computational environment, then decompose an image the same way and reconstruct a recognizable picture from a small percentage of its singular values, seeing exactly what was lost (the discarded layers). ## Introduction A digital image is just a grid of numbers. A [grayscale](https://en.wikipedia.org/wiki/Grayscale) photo is a matrix whose entry $A_{ij}$ is the brightness of [pixel](https://en.wikipedia.org/wiki/Pixel) $(i,j)$, typically an integer from $0$ (black) to $255$ (white). A $512\times512$ image is therefore a matrix with $262{,}144$ entries. Storing every one of them is expensive, but real images are far from random: neighboring pixels are correlated, edges and regions repeat, and large smooth areas carry little independent information. That redundancy is exactly what linear algebra can expose and exploit. The tool is the [singular value decomposition](https://en.wikipedia.org/wiki/Singular_value_decomposition) (SVD), which writes *any* matrix as a sum of simple rank-one layers, ordered from most to least important. Keeping only the first few layers gives a [low-rank approximation](https://en.wikipedia.org/wiki/Low-rank_approximation): a nearby matrix that captures the dominant structure of the image while discarding the fine detail that costs the most to store. The following shows this on a picture of 8-bit Link, where the number in the upper right is the percentage of the singular values kept. The final frame of this sweep (cleaned up) is [`link_gray.png`](data/link_gray.png), the very image you will decompose below. ![Rank sweep on 8-bit Link; the counter is the percentage of singular values kept](Media/link.gif) ## Project Description The project moves in four passes. First, work the decomposition by hand on small matrices, confirm every number computationally, and read the geometry out of the results (orthonormal directions, stretch factors). Second, mirror the same decomposition on the 8-bit Link sprite: compute its SVD, watch the singular values decay, and recover low-order approximations. Third, generalize to color on a photograph of Link's plush (an RGB image is three matrices, so the same tool runs three times). Finally, do it on three pictures of your own, chosen for three different spectral personalities. **Foundation task:** Perform the hand calculations in the Mathematical Background (a symmetric eigendecomposition and a complete $3\times2$ SVD, run all the way through dropping a layer), then confirm every number in your language of choice and state the geometric conclusions. **Application task:** Mirror the decomposition on [`link_gray.png`](data/link_gray.png): compute the SVD, plot the decay and the energy captured, reconstruct at $0.2\%$, $1\%$, $4\%$, and $10\%$ of the singular values, confirm the $100\%$ reconstruction to machine precision, and report the energy captured and average pixel error at each percentage. Then generalize to [`link_plush.png`](data/link_plush.png): decompose each RGB channel, compare the three spectra to each other and to the sprite's, and reconstruct in color. **Key deliverable:** The same study on three images you choose with three spectral personalities (one that should compress easily, one with strong variance across the color channels, one that should resist compression): each spectrum presented and interpreted, your own justified percentages per image displayed alongside $100\%$, storage accounting, and a cross-image compressibility ranking defended from the spectra. Two companions support the milestone forms: the [[MATH307Su26 - Low Rank Approximations of Images (Milestone Map)|Milestone Map]] says which artifact from this handout answers each form field, and the [[MATH307Su26 - Low Rank Approximations of Images (Verification Table)|Verification Table]] is filled in as you verify and submitted with Milestone 2. ## Mathematical Background: from eigenvalues to singular values **Idea.** This is the eigenvalue and diagonalization machinery from our linear-algebra thread, extended to non-square matrices. A [symmetric matrix](https://en.wikipedia.org/wiki/Symmetric_matrix) can be written as a weighted sum of rank-one pieces built from its eigenvectors (its [spectral decomposition](https://en.wikipedia.org/wiki/Spectral_theorem)). The SVD is the same idea made to work for any rectangular matrix, and an image is a rectangular matrix. **Recall.** If $A$ is a square matrix, an [**eigenvector**](https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors) $\mathbf{x}$ and its corresponding **eigenvalue** $\lambda$ satisfy: $A\mathbf{x} = \lambda\mathbf{x}$ This means that multiplying $A$ by $\mathbf{x}$ stretches or compresses $\mathbf{x}$ by a factor of $\lambda$. The directions of these special vectors $\mathbf{x}$ remain unchanged by the transformation. For a 2×2 matrix, eigenvalues are found by solving the [characteristic equation](https://en.wikipedia.org/wiki/Characteristic_polynomial): $\det(A - \lambda I) = 0$ Then for each $\lambda$, you solve $(A - \lambda I)\mathbf{x} = \mathbf{0}$ to find the corresponding eigenvector $\mathbf{x}$. Theory tells us that when there are as many eigenvectors as there are columns of $A$, then it is possible to use them as the columns of a matrix $P$ such that $ A = P D P^{-1}$ where $D$ is a diagonal matrix whose elements are the eigenvalues of $A$ (such an $A$ is called [diagonalizable](https://en.wikipedia.org/wiki/Diagonalizable_matrix)). If $A$ has the additional property of being a symmetric matrix, $A=A^{T}$, then the [spectral theorem](https://en.wikipedia.org/wiki/Spectral_theorem) gives [orthonormal](https://en.wikipedia.org/wiki/Orthonormality) eigenvectors $\mathbf{x}_i$ and real eigenvalues $\lambda_i$ with a simplified diagonal decomposition and summation formula, $A = Q D Q^{T} = \sum_{i} \lambda_i\, \mathbf{x}_i \mathbf{x}_i^{T},$ which is a sum of rank-one pieces, each a direction $\mathbf{x}_i$ scaled by $\lambda_i$. **Worked example (symmetric case - see course notes for non-symmetric case).** For $A=\begin{bmatrix}0&1\\1&0\end{bmatrix}$ the eigenpairs are $\lambda_1=1,\ \mathbf{x}_1=\tfrac{1}{\sqrt2}(1,1)$ and $\lambda_2=-1,\ \mathbf{x}_2=\tfrac{1}{\sqrt2}(1,-1)$. The spectral sum reproduces $A$ exactly, $1\cdot\tfrac12\begin{bmatrix}1&1\\1&1\end{bmatrix}+(-1)\cdot\tfrac12\begin{bmatrix}1&-1\\-1&1\end{bmatrix}=\begin{bmatrix}0&1\\1&0\end{bmatrix},$ which is the first check you will confirm in code (Step 1 below). Note the geometry: the eigenvectors are perpendicular, and $A$ (a reflection across the line $y=x$) acts by stretching one diagonal direction by $+1$ and the other by $-1$. **The singular value decomposition.** An image matrix is rectangular, so it has no eigendecomposition. Even so, every matrix carries a spectral-like structure: every real matrix $A\in\mathbb{R}^{m\times n}$ admits a **singular value decomposition** (SVD), $A = U\Sigma V^{T} = \sum_{i=1}^{r}\sigma_i\, \mathbf{u}_i \mathbf{v}_i^{T},$ where $U$ and $V$ have orthonormal columns (the [left and right singular vectors](https://en.wikipedia.org/wiki/Singular_value_decomposition)), and $\Sigma$ is diagonal with the [singular values](https://en.wikipedia.org/wiki/Singular_value) $\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0$ in decreasing order. Each term $\sigma_i\mathbf{u}_i\mathbf{v}_i^{T}$ is a rank-one layer (an [outer product](https://en.wikipedia.org/wiki/Outer_product) scaled by $\sigma_i$), and the ordering means the first layers carry the most of the image. Geometrically the three factors act as [rotate](https://en.wikipedia.org/wiki/Rotation_matrix) ($V^{T}$), stretch along axes ($\Sigma$), rotate ($U$): the SVD says every matrix, however lopsided, is a rotation-stretch-rotation. > [!note] Where this comes from > The SVD is the spectral theorem applied to $A^{T}A$. That matrix is symmetric and [positive semidefinite](https://en.wikipedia.org/wiki/Definite_matrix), so it has orthonormal eigenvectors (the columns of $V$) and nonnegative eigenvalues, and the singular values are $\sigma_i=\sqrt{\lambda_i(A^{T}A)}$. Nothing new is introduced (the SVD is the symmetric-eigenvalue story of $A^{T}A$, repackaged for a rectangular $A$). **Worked example (rectangular case, run to completion).** Let $A=\begin{bmatrix}2&0\\0&2\\-2&0\end{bmatrix}$. Then $A^{T}A=\begin{bmatrix}8&0\\0&4\end{bmatrix}$ and the eigenvalues of $A^T A$ are 8 and 4. So, the singular values of $A$ are $ \sigma_1 = \sqrt{8} = 2\sqrt{2}, \quad \sigma_2 = \sqrt{4} = 2 $ Because $A^T A$ is diagonal, the eigenvectors (right singular vectors) are simply the standard basis vectors: $ V = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} $ Now compute the left singular vectors using $\mathbf{u}_i = \frac{1}{\sigma_i} A \mathbf{v}_i$: $ \mathbf{u}_1 = \frac{1}{2\sqrt{2}} A \begin{bmatrix} 1 \\ 0 \end{bmatrix} = \frac{1}{2\sqrt{2}} \begin{bmatrix} 2 \\ 0 \\ -2 \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 0 \\ -1 \end{bmatrix} $ $ \mathbf{u}_2 = \frac{1}{2} A \begin{bmatrix} 0 \\ 1 \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 0 \\ 2 \\ 0 \end{bmatrix} = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} $ So the left singular vectors form the matrix $ U = \begin{bmatrix} \frac{1}{\sqrt{2}} & 0 \\ 0 & 1 \\ -\frac{1}{\sqrt{2}} & 0 \end{bmatrix} $ The diagonal matrix of singular values is $ \Sigma = \begin{bmatrix} 2\sqrt{2} & 0 \\ 0 & 2 \\ 0 & 0 \end{bmatrix} $ Putting it all together, $ A = U \Sigma V^T = \begin{bmatrix} \frac{1}{\sqrt{2}} & 0 \\ 0 & 1 \\ -\frac{1}{\sqrt{2}} & 0 \end{bmatrix} \begin{bmatrix} 2\sqrt{2} & 0 \\ 0 & 2 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} $ This decomposition reveals the dominant direction (along the x-axis) with a strong singular value $2\sqrt{2}$, and a secondary vertical direction with singular value 2. The third row of $A$ mirrors the first, and the structure of $U$ records it: check that the columns of $U$ are orthonormal (each has unit length, and their dot product is zero). **Continuing the example: layers, and what dropping one costs.** The same decomposition written as rank-one layers is $A = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T + \sigma_2 \mathbf{u}_2 \mathbf{v}_2^T$. Compute each term explicitly: $ \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T = 2\sqrt{2} \cdot \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 0 \\ -1 \end{bmatrix} \begin{bmatrix} 1 & 0 \end{bmatrix} = 2 \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ -1 & 0 \end{bmatrix}, \qquad \sigma_2 \mathbf{u}_2 \mathbf{v}_2^T = 2 \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix} \begin{bmatrix} 0 & 1 \end{bmatrix} = 2 \begin{bmatrix} 0 & 0 \\ 0 & 1 \\ 0 & 0 \end{bmatrix} $ Their sum is $A$ exactly: keeping $100\%$ of the layers reconstructs the matrix. Now keep only the dominant layer, $ A_1 = \sigma_1 \mathbf{u}_1 \mathbf{v}_1^T = \begin{bmatrix} 2 & 0 \\ 0 & 0 \\ -2 & 0 \end{bmatrix}, \qquad A - A_1 = \begin{bmatrix} 0 & 0 \\ 0 & 2 \\ 0 & 0 \end{bmatrix}, $ and the error is exactly the discarded layer: $A - A_1 = \sigma_2\mathbf{u}_2\mathbf{v}_2^{T}$, entry for entry. Nothing about the kept layer appears in the error, and its size is set entirely by $\sigma_2$. This is the identity your warm-up code must reproduce, and it is the whole compression story in miniature: what a truncation costs is exactly the layers you threw away. ### Truncation: the discarded layers are the error Keeping the first $k$ layers defines the rank-$k$ truncation $\boxed{\,A_k=\sum_{i=1}^{k}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^{T},\,}$ and, exactly as in the hand example, subtraction leaves the discarded layers and nothing else: $A-A_k=\sum_{i=k+1}^{r}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^{T}.$ So the cost of truncation is set by the singular values you dropped: when $\sigma_{k+1},\dots,\sigma_r$ are small, so is the error. How much was lost gets reported two ways (the note below), and $100\%$ is the built-in check: keeping all $r$ layers must reproduce $A$ to [machine precision](https://en.wikipedia.org/wiki/Machine_epsilon). > [!note] Two percentages (do not conflate them) > The task dial in this project is the **percentage of singular values kept**, $100\,k/r$. A more telling measure is the **energy captured**, $100\sum_{i\le k}\sigma_i^2 / \sum_i\sigma_i^2$, the share of the total $\sigma^2$ content carried by the kept layers. The two diverge dramatically: on our sprite, keeping $0.2\%$ of the singular values ($k=1$) already captures $79.1\%$ of the energy. Report both, along with the **average pixel error** (mean absolute difference in gray levels), which is the number your eyes can sanity-check against the panels. > [!warning] Ordering matters (a common trap) > Singular values are indexed in *decreasing* order, $\sigma_1\ge\sigma_2\ge\cdots$, and a good approximation keeps the *first* (largest) $k$ terms. Keeping the smallest terms instead would discard the image and retain the noise. Confirm your language returns the singular values largest-first (all four below do). ## Implementation Guidelines The workflow is the same in every language, so we state it once and then implement it step by step. > [!abstract] The workflow (state it once, run it in any language) > Given an image and target percentages $p_1 < p_2 < p_3 < 100$ of singular values to keep: > 1. **Warm up**: define the verification helper, then reproduce the worked examples numerically (the $2\times2$ eigendecomposition, the $3\times2$ SVD, and $A - A_1 = \sigma_2\mathbf{u}_2\mathbf{v}_2^{T}$). > 2. **Load** the image as a matrix and cast to floating point. > 3. **Decompose**: take the economy SVD $A = U\Sigma V^{T}$ and inspect the decay of $\sigma_1 \ge \sigma_2 \ge \cdots$. > 4. **Truncate and check**: for each target percentage, $k = \operatorname{round}(p\,r/100)$ and $A_k = U_k \Sigma_k V_k^{T}$; confirm $100\%$ reproduces $A$ to machine precision; report the energy captured $\sum_{i\le k}\sigma_i^2 / \sum_i \sigma_i^2$, the average pixel error, and the storage $k(m+n+1)/(mn)$. > 5. **Render**: the reconstruction panels and the energy curve. > 6. **Generalize to color**: an RGB image is three matrices; run 2 through 5 per channel and restack (its own section after the verification warning below). Pick your language below; each callout contains all five steps in order, run top to bottom. Code is four-language and portable (Python is the run-verified reference; Mathematica is numerical only, i.e., no symbolic linear algebra). Download the two reference images from the course website (links on the Assignments page) and put them next to your script; [`link_gray.png`](data/link_gray.png) is a single-channel 8-bit PNG, so every language reads the identical pixel matrix. > [!example]- Python (reference) > > **Step 1: define the helper, then warm up on the hand examples.** The helper compares two matrices and reports the largest entry-by-entry difference; every check in this project runs through it. > ```python > import numpy as np > def same(name, X, Y, tol=1e-9): # verification helper > d = np.max(np.abs(np.asarray(X) - np.asarray(Y))) > print(f"{name}: max difference {d:.1e}", "OK" if d < tol else "CHECK THIS") > A2 = np.array([[0., 1.], [1., 0.]]) > w, Q = np.linalg.eigh(A2) # ascending: [-1, 1] > print("eigenvalues:", w) > same("spectral rebuild", Q @ np.diag(w) @ Q.T, A2) > A = np.array([[2., 0.], [0., 2.], [-2., 0.]]) > U3, s3, Vt3 = np.linalg.svd(A, full_matrices=False) > print("singular values:", s3) # [2.828, 2.] = [2*sqrt(2), 2] > A1 = s3[0] * np.outer(U3[:, 0], Vt3[0]) # the dominant layer > layer2 = s3[1] * np.outer(U3[:, 1], Vt3[1]) # the discarded layer > same("A - A_1 equals the discarded layer", A - A1, layer2) > ``` > *Look for:* eigenvalues $\pm1$ (ascending order in `eigh`), singular values $2\sqrt2, 2$, and both helper checks printing OK with differences at roundoff ($\sim10^{-16}$). > > **Step 2: load the image as a matrix.** PIL reads the PNG as 8-bit integers; cast to float before any arithmetic. > ```python > import matplotlib.pyplot as plt > from PIL import Image > M = np.asarray(Image.open("link_gray.png")).astype(float) # 512x512, values 0..255 > print("shape:", M.shape) # (512, 512) > plt.imshow(M, cmap="gray", vmin=0, vmax=255); plt.axis("off"); plt.show() > ``` > *Look for:* shape `(512, 512)` and values spanning 0..255 (floats, not uint8). > > **Step 3: compute the SVD and inspect the decay.** `np.linalg.svd` returns `Vt`, which is already $V^{T}$. > ```python > U, s, Vt = np.linalg.svd(M, full_matrices=False) # s descending > print("sigma_1..6:", np.round(s[:6], 3)) # [97205.807 25235.034 21359.274 20216.085 15705.16 13168.314] > plt.semilogy(np.arange(1, s.size+1), np.maximum(s, 1e-14), ); plt.xlabel("i"); plt.ylabel("sigma_i"); plt.show() > ``` > *Look for:* $\sigma_1\approx97205.8$ towering over the rest, then a cliff at $i=83$: the sprite is *exactly* rank $83$ ($\sigma_{84}\approx10^{-10}$ is pure roundoff). Pixel art is genuinely low rank; photographs only approximately so. The `np.maximum(..., 1e-14)` keeps the log plot from choking on the zero tail. > > **Step 4: truncate at the target percentages and check.** $k = 1, 5, 20, 50$ is $0.2\%, 1\%, 4\%, 10\%$ of $r=512$; the average pixel error says how far off a typical pixel is (out of 255), and $100\%$ must reproduce the image to machine precision. > ```python > rank_k = lambda k: (U[:, :k] * s[:k]) @ Vt[:k, :] # U_k diag(s_k) V_k^T > tot = np.sum(s**2) > for k in (1, 5, 20, 50): > Ak = rank_k(k) > energy = 100*np.sum(s[:k]**2)/tot > avg = np.mean(np.abs(M - Ak)) # average pixel error, gray levels > print(f"k={k:3d} ({100*k/s.size:5.2f}%) energy {energy:5.2f}% avg pixel error {avg:5.2f}") > same("100% rebuilds the image", rank_k(s.size), M, tol=1e-8) > ``` > *Look for:* energy $79.08\%, 93.72\%, 99.67\%, 99.97\%$; average pixel errors $85.21, 38.44, 5.04, 1.10$ gray levels; the $100\%$ helper check OK at $\sim10^{-10}$. > > **Step 5: render the panels and the energy curve.** Clip reconstructions to $[0,255]$ (low-rank output can overshoot the pixel range); the energy curve is pure singular-value arithmetic, no reconstructions needed. > ```python > ks, sig = [1, 5, 20, 50], ["0.2%", "1%", "4%", "10%"] > cum = np.cumsum(s**2) # energy, for the dual labels > fig, ax = plt.subplots(1, 5, figsize=(15, 3.6)) > ax[0].imshow(M, cmap="gray", vmin=0, vmax=255); ax[0].set_title("original (100%)"); ax[0].axis("off") > for a, k, sg in zip(ax[1:], ks, sig): > a.imshow(np.clip(rank_k(k), 0, 255), cmap="gray", vmin=0, vmax=255) > a.set_title(f"k = {k} ({sg} of sigma, {100*cum[k-1]/tot:.1f}% energy)"); a.axis("off") > plt.tight_layout(); plt.show() > plt.plot(np.arange(1, s.size+1), 100*cum/tot); plt.axhline(95, ls="--", c="gray") > plt.xlabel("k"); plt.ylabel("energy captured (%)"); plt.show() > ``` > *Look for:* rank 20 already recognizable; the energy curve crosses $95\%$ already at $k=6$. > [!example]- MATLAB > > **Step 1: define the helper, then warm up on the hand examples.** The helper reports the largest entry-by-entry difference between two matrices; every check runs through it. > ```matlab > same = @(name, X, Y) fprintf('%s: max difference %.1e\n', name, max(abs(X - Y), [], 'all')); > A2 = [0 1; 1 0]; > [Q, D] = eig(A2); disp(diag(D).') % ascending: -1 1 > same('spectral rebuild', Q*D*Q', A2) > A = [2 0; 0 2; -2 0]; > [U3, S3, V3] = svd(A, 'econ'); s3 = diag(S3); > disp(s3.') % 2.8284 2.0000 > A1 = s3(1) * U3(:,1) * V3(:,1)'; % the dominant layer > same('A - A1 equals the discarded layer', A - A1, s3(2) * U3(:,2) * V3(:,2)') > ``` > *Look for:* eigenvalues $\pm1$, singular values $2\sqrt2, 2$, and both helper differences at roundoff ($\sim10^{-16}$). > > **Step 2: load the image as a matrix.** `imread` gives `uint8`; cast with `double` before any arithmetic. > ```matlab > M = double(imread('link_gray.png')); % uint8 -> double, values 0..255 > disp(size(M)) % [512 512] > imshow(uint8(M)) % display expects uint8 or [] scaling > ``` > *Look for:* size `[512 512]`; forgetting `double()` makes later arithmetic silently integer-saturate. > > **Step 3: compute the SVD and inspect the decay.** MATLAB returns $V$ (not $V^{T}$); the transpose happens at reconstruction. > ```matlab > [U, S, V] = svd(M, 'econ'); s = diag(S); % economy SVD; s descending > disp(s(1:6).') % 97205.8 25235.0 21359.3 20216.1 15705.2 13168.3 > semilogy(1:numel(s), max(s, 1e-14)); xlabel('i'); ylabel('\sigma_i') > ``` > *Look for:* the same six leading values as the Python reference, and the cliff at $i=83$ (the sprite is exactly rank $83$). > > **Step 4: truncate at the target percentages and check.** Note the explicit transpose `V(:,1:k)'`. > ```matlab > r = numel(s); tot = sum(s.^2); > for k = [1 5 20 50] > Ak = U(:,1:k) * diag(s(1:k)) * V(:,1:k)'; % note V, transposed here > fprintf('k=%3d (%5.2f%%) energy %5.2f%% avg pixel error %5.2f\n', ... > k, 100*k/r, 100*sum(s(1:k).^2)/tot, mean(abs(M - Ak), 'all')); > end > same('100% rebuilds the image', U * diag(s) * V', M) > ``` > *Look for:* energy $79.08, 93.72, 99.67, 99.97\%$; average pixel errors $85.21, 38.44, 5.04, 1.10$; the $100\%$ helper difference at $\sim10^{-10}$. > > **Step 5: render the panels and the energy curve.** `uint8()` saturates (clamps) automatically; export with `exportgraphics`, not `saveas`. > ```matlab > ks = [1 5 20 50]; sig = {'0.2%', '1%', '4%', '10%'}; > cum = cumsum(s.^2); % energy, for the dual labels > figure('Position', [100 100 1500 360]); > subplot(1,5,1); imshow(uint8(M)); title('original (100%)') > for i = 1:4 > k = ks(i); Ak = U(:,1:k) * diag(s(1:k)) * V(:,1:k)'; > subplot(1,5,i+1); imshow(uint8(Ak)) > title(sprintf('k = %d (%s of \\sigma, %.1f%% energy)', k, sig{i}, 100*cum(k)/tot)) > end > figure; plot(100*cum/tot); yline(95, '--'); xlabel('k'); ylabel('energy captured (%)') > ``` > *Look for:* five panels, original first; the energy curve crosses $95\%$ already at $k=6$. > [!example]- R > > **Step 1: define the helper, then warm up on the hand examples.** The helper reports the largest entry-by-entry difference between two matrices; every check runs through it. > ```r > same <- function(name, X, Y) cat(sprintf("%s: max difference %.1e\n", name, max(abs(X - Y)))) > A2 <- matrix(c(0, 1, 1, 0), 2, 2) > e <- eigen(A2) # e$values DECREASING: 1, -1 > print(e$values) > same("spectral rebuild", e$vectors %*% diag(e$values) %*% t(e$vectors), A2) > A <- matrix(c(2, 0, -2, 0, 2, 0), 3, 2) # column-major -> [2 0; 0 2; -2 0] > sv3 <- svd(A) > print(sv3$d) # 2.828427 2.000000 > A1 <- sv3$d[1] * sv3$u[, 1] %*% t(sv3$v[, 1]) # the dominant layer > same("A - A1 equals the discarded layer", A - A1, sv3$d[2] * sv3$u[, 2] %*% t(sv3$v[, 2])) > ``` > *Look for:* R orders eigenvalues *decreasing* ($1, -1$), the reverse of Python and MATLAB; singular values $2\sqrt2, 2$; both helper differences at roundoff ($\sim10^{-16}$). > > **Step 2: load the image as a matrix.** `png::readPNG` returns values in $[0,1]$, so scale by 255; base `image()` draws matrices rotated, hence the flip. > ```r > library(png) # install.packages("png") if needed > M <- png::readPNG("link_gray.png") * 255 # to the common 0..255 range > cat(dim(M), "\n") # 512 512 > image(t(M[nrow(M):1, ]), col = gray.colors(256), asp = 1) # base R image() flips rows > ``` > *Look for:* `dim` 512 512; if the sprite is sideways or mirrored, the flip is missing. > > **Step 3: compute the SVD and inspect the decay.** `svd()` returns `sv$d` (descending), `sv$u`, `sv$v`. > ```r > sv <- svd(M) > print(round(sv$d[1:6], 3)) # 97205.807 25235.034 21359.274 20216.085 15705.160 13168.314 > plot(pmax(sv$d, 1e-14), log = "y", type = "l", xlab = "i", ylab = expression(sigma[i])) > ``` > *Look for:* the same six leading values as the Python reference, and the cliff at $i=83$ (the sprite is exactly rank $83$); `pmax` keeps the log axis happy on the zero tail. > > **Step 4: truncate at the target percentages and check.** `diag(x, k, k)` keeps a $1\times1$ matrix from collapsing to a scalar when $k=1$. > ```r > r <- length(sv$d); tot <- sum(sv$d^2) > rank_k <- function(k) sv$u[, 1:k] %*% diag(sv$d[1:k], k, k) %*% t(sv$v[, 1:k]) > for (k in c(1, 5, 20, 50)) > cat(sprintf("k=%3d (%5.2f%%) energy %5.2f%% avg pixel error %5.2f\n", > k, 100*k/r, 100*sum(sv$d[1:k]^2)/tot, mean(abs(M - rank_k(k))))) > same("100% rebuilds the image", rank_k(r), M) > ``` > *Look for:* energy $79.08, 93.72, 99.67, 99.97\%$; average pixel errors $85.21, 38.44, 5.04, 1.10$; dropping the `diag(., k, k)` guard breaks only $k=1$, a classic silent bug. > > **Step 5: render the panels and the energy curve.** Clamp with `pmin/pmax`; use the `draw()` closure plus `png(); draw(); dev.off()` if you also want the figure saved under headless `Rscript`. > ```r > ks <- c(1, 5, 20, 50); sig <- c("0.2%", "1%", "4%", "10%") > cum <- cumsum(sv$d^2) # energy, for the dual labels > par(mfrow = c(1, 5), mar = c(0, 0, 2, 0)) > show <- function(A, t) image(t(A[nrow(A):1, ])/255, col = gray.colors(256), > axes = FALSE, asp = 1, zlim = c(0, 1), main = t) > show(M, "original (100%)") > for (i in 1:4) show(pmin(pmax(rank_k(ks[i]), 0), 255), > sprintf("k = %d (%s of sigma, %.1f%% energy)", ks[i], sig[i], 100*cum[ks[i]]/tot)) > par(mfrow = c(1, 1), mar = c(4, 4, 1, 1)) > plot(100*cum/tot, type = "l", xlab = "k", ylab = "energy captured (%)") > abline(h = 95, lty = 2) > ``` > *Look for:* five panels, original first; the energy curve crosses $95\%$ already at $k=6$. > [!example]- Mathematica (numerical only) > > **Step 1: define the helper, then warm up on the hand examples.** The helper reports the largest entry-by-entry difference between two matrices; every check runs through it. (`E`, `N`, `D` are protected names; use your own.) > ```wolfram > same[name_, x_, y_] := Print[name, ": max difference ", Max@Abs[x - y]]; > A2 = {{0., 1.}, {1., 0.}}; > {vals, vecs} = Eigensystem[A2]; (* {1., -1.} with row eigenvectors *) > Print[vals] > same["spectral rebuild", Transpose[vecs] . DiagonalMatrix[vals] . vecs, A2] > A = {{2., 0.}, {0., 2.}, {-2., 0.}}; > {u3, sig3, v3} = SingularValueDecomposition[A]; > Print[Diagonal[sig3]] (* {2.82843, 2.} *) > a1 = sig3[[1, 1]] Outer[Times, u3[[All, 1]], v3[[All, 1]]]; > same["A - A1 equals the discarded layer", A - a1, sig3[[2, 2]] Outer[Times, u3[[All, 2]], v3[[All, 2]]]] > ``` > *Look for:* eigenvalues as $\{1., -1.\}$ (Mathematica returns eigenvectors as *rows*); singular values $2\sqrt2, 2$; both helper differences at roundoff ($\sim10^{-16}$). > > **Step 2: load the image as a matrix.** `ImageData` returns values in $[0,1]$, so scale by 255 to the common range; keep everything numeric. > ```wolfram > img = Import["link_gray.png"]; (* numerical only *) > M = N[ImageData[img]]*255.; (* to the common 0..255 range *) > Print[Dimensions[M]] (* {512, 512} *) > Image[M/255.] > ``` > *Look for:* `{512, 512}`; the `N[...]` keeps the pipeline in machine numbers. > > **Step 3: compute the SVD and inspect the decay.** `SingularValueDecomposition` returns $\Sigma$ as a matrix; take `Diagonal`. > ```wolfram > {u, sig, v} = SingularValueDecomposition[M]; (* numerical; sig is a diagonal matrix *) > s = Diagonal[sig]; (* descending singular values *) > Print[N@Take[s, 6]] (* {97205.8, 25235., 21359.3, 20216.1, 15705.2, 13168.3} *) > ListLogPlot[Clip[s, {10.^-14, Infinity}], Joined -> True, AxesLabel -> {"i", "sigma"}] > ``` > *Look for:* the same six leading values as the Python reference, and the cliff at $i=83$ (the sprite is exactly rank $83$); the `Clip` keeps the log axis happy on the zero tail. > > **Step 4: truncate at the target percentages and check.** Slices use `1;;k`; `v` holds columns, so `Transpose` at reconstruction. > ```wolfram > r = Length[s]; tot = Total[s^2]; > rankK[k_] := u[[All, 1;;k]] . sig[[1;;k, 1;;k]] . Transpose[v[[All, 1;;k]]]; > Do[Module[{Ak = rankK[k]}, > Print[{k, 100. k/r, 100. Total[s[[1;;k]]^2]/tot, Mean[Flatten[Abs[M - Ak]]]}]], > {k, {1, 5, 20, 50}}] > same["100% rebuilds the image", rankK[r], M] > ``` > *Look for:* energy $79.08, 93.72, 99.67, 99.97$; average pixel errors $85.21, 38.44, 5.04, 1.10$; the $100\%$ helper difference at $\sim10^{-10}$. > > **Step 5: render the panels and the energy curve.** Clamp with `Clip`; `Export` saves files but does not display, so evaluate the graphics on their own lines to see them. > ```wolfram > ks = {1, 5, 20, 50}; sig = {"0.2%", "1%", "4%", "10%"}; > cum = Accumulate[s^2]; (* energy, for the dual labels *) > GraphicsRow[Prepend[ > Table[Labeled[Image[Clip[rankK[ks[[i]]], {0., 255.}]/255.], > "k = " <> ToString[ks[[i]]] <> " (" <> sig[[i]] <> " of sigma, " <> > ToString[NumberForm[100. cum[[ks[[i]]]]/tot, {4, 1}]] <> "% energy)"], {i, 4}], > Labeled[Image[M/255.], "original (100%)"]], ImageSize -> 900] > ListLinePlot[100. cum/tot, AxesLabel -> {"k", "energy captured (%)"}, > GridLines -> {None, {95}}] > ``` > *Look for:* five panels, original first; the energy curve crosses $95\%$ already at $k=6$. What the guided run should produce: ![Reconstruction panels; each title carries both dials, the share of singular values kept and the energy captured](Media/lowrank_panels.png) ![Singular-value decay (cliff at the exact rank 83), and energy captured versus k (95% at k = 6)](Media/lowrank_spectrum.png) > [!warning] Verify against ground truth (required) > Record these checked numbers; "it rendered" is not verification. > - Warm-up: eigenvalues $\pm1$; singular values $2\sqrt2\approx2.828,\ 2$; both helper checks (spectral rebuild, $A-A_1$ equals the discarded layer) at roundoff, $\sim10^{-16}$. > - Sprite ($512\times512$): $\sigma_{1..6} \approx 97205.8,\ 25235.0,\ 21359.3,\ 20216.1,\ 15705.2,\ 13168.3$. > - The decay cliff: $\sigma_{83}\approx39.0$ but $\sigma_{84}\approx5\times10^{-10}$; the sprite is exactly rank $83$. > - Energy captured: $79.08\%$ ($k=1$), $93.72\%$ ($k=5$), $99.67\%$ ($k=20$), $99.97\%$ ($k=50$); $95\%$ first reached at $k=6$. > - Average pixel error (gray levels of 255): $85.21,\ 38.44,\ 5.04,\ 1.10$ at $k=1,5,20,50$. > - $100\%$ reconstruction: helper difference $\approx 4\times10^{-10}$ (machine precision). > - Storage at rank $k$ is $k(m+n+1)$ versus $mn=262{,}144$: $0.39\%,\ 1.96\%,\ 7.82\%,\ 19.55\%$ at $k=1,5,20,50$; the exact rank $83$ costs $32.5\%$. ### Generalization: color is three matrices (Step 6) So far the image was one matrix because each pixel held one number, a gray level. A color pixel holds three: the brightnesses of a red, a green, and a blue subpixel, 8 bits each (this is [RGB](https://en.wikipedia.org/wiki/RGB_color_model), and "24-bit color" is exactly these three bytes). A color image is therefore an $H\times W\times 3$ array, i.e., three grayscale matrices photographing the same scene: a red matrix, a green matrix, and a blue matrix. The SVD is a matrix tool, so the generalization asks for nothing new: decompose each channel separately, truncate each at the same $k$, and restack the three reconstructions into a color image. The shipped photograph is [`link_plush.png`](data/link_plush.png), Link again, this time as a $768\times512$ photo of the plush (each channel is a $768\times512$ matrix, so $r=512$ as before). Two things to watch. First, the spectra: the sprite's singular values fell off a cliff at its exact [rank](https://en.wikipedia.org/wiki/Rank_%28linear_algebra%29) $83$, while a photograph decays smoothly with no cliff (effective rank $494$ of $512$ here), and the three channels decay at slightly different rates (blue is slowest: the tunic carries the detail). Second, at very small $k$ the three channels disagree about where the edges are, which shows up as color fringing; it fades as $k$ grows. > [!example]- Step 6 code, Python (reference) > > **Step 6: decompose each channel, truncate, restack.** The third array index selects the channel; everything else is Steps 3 and 4 run three times. > ```python > rgb = np.asarray(Image.open("link_plush.png")).astype(float) # 768 x 512 x 3 > print("shape:", rgb.shape) > def channel_k(chan, k): > Uc, sc, Vtc = np.linalg.svd(chan, full_matrices=False) > return (Uc[:, :k] * sc[:k]) @ Vtc[:k, :] > for k in (5, 20, 50): > out = np.stack([channel_k(rgb[:, :, c], k) for c in range(3)], axis=2) > print(k, round(np.mean(np.abs(rgb - out)), 2)) # 24.72 10.9 5.44 > show = np.clip(np.stack([channel_k(rgb[:, :, c], 50) for c in range(3)], axis=2), 0, 255) > plt.imshow(show.astype(np.uint8)); plt.axis("off"); plt.show() > for c, col in zip(range(3), ("crimson", "green", "royalblue")): # channel spectra > sc = np.linalg.svd(rgb[:, :, c], compute_uv=False) > plt.semilogy(sc/sc[0], color=col) > plt.xlabel("i"); plt.ylabel("sigma_i / sigma_1"); plt.show() > ``` > *Look for:* average pixel errors $24.72, 10.90, 5.44$ at $k=5,20,50$; pink/cyan fringing at $k=5$ where the channels disagree; three smooth spectra with no cliff (put the sprite's on the same axes and the contrast is unmissable). > [!example]- Step 6 code, MATLAB > > **Step 6: decompose each channel, truncate, restack.** The channel is the third index; write each reconstructed channel into a preallocated array. > ```matlab > rgb = double(imread('link_plush.png')); % 768 x 512 x 3 > disp(size(rgb)) > out = zeros(size(rgb)); > for k = [5 20 50] > for c = 1:3 > [Uc, Sc, Vc] = svd(rgb(:,:,c), 'econ'); > out(:,:,c) = Uc(:,1:k) * Sc(1:k,1:k) * Vc(:,1:k)'; > end > fprintf('k=%d avg pixel error %.2f\n', k, mean(abs(rgb - out), 'all')); > end > figure; imshow(uint8(out)) % the k = 50 reconstruction > figure; hold on > for c = 1:3 > sc = svd(rgb(:,:,c)); plot(sc/sc(1)); > end > set(gca, 'YScale', 'log'); xlabel('i'); ylabel('\sigma_i/\sigma_1'); legend('R','G','B') > ``` > *Look for:* the same three error numbers as the Python reference; `uint8` saturation handles the clipping on display. > [!example]- Step 6 code, R > > **Step 6: decompose each channel, truncate, restack.** The channel is the third index of the array. > ```r > rgb <- png::readPNG("link_plush.png") * 255 # 768 x 512 x 3 > cat(dim(rgb), "\n") > channel_k <- function(chan, k) { > svc <- svd(chan) > svc$u[, 1:k] %*% diag(svc$d[1:k], k, k) %*% t(svc$v[, 1:k]) > } > out <- array(0, dim(rgb)) > for (k in c(5, 20, 50)) { > for (c in 1:3) out[, , c] <- channel_k(rgb[, , c], k) > cat(sprintf("k=%d avg pixel error %.2f\n", k, mean(abs(rgb - out)))) > } > plot(as.raster(pmin(pmax(out, 0), 255)/255)) # the k = 50 reconstruction > matplot(sapply(1:3, function(c) { d <- svd(rgb[, , c])$d; d/d[1] }), > type = "l", log = "y", lty = 1, col = c("red", "green3", "blue"), > xlab = "i", ylab = "sigma_i / sigma_1") > ``` > *Look for:* the same error numbers; `as.raster` wants values in $[0,1]$, hence the divide by 255. > [!example]- Step 6 code, Mathematica (numerical only) > > **Step 6: decompose each channel, truncate, restack.** `SingularValueList` gives just the spectrum; `C` is a protected name, hence `chan`. > ```wolfram > rgb = N[ImageData[Import["link_plush.png"]]]*255.; (* 768 x 512 x 3 *) > Print[Dimensions[rgb]] > channelK[chan_, k_] := Module[{uc, sc, vc}, > {uc, sc, vc} = SingularValueDecomposition[chan]; > uc[[All, 1;;k]] . sc[[1;;k, 1;;k]] . Transpose[vc[[All, 1;;k]]]]; > restack[k_] := Transpose[Table[channelK[rgb[[All, All, c]], k], {c, 3}], {3, 1, 2}]; > Do[Print[{k, Mean[Flatten[Abs[rgb - restack[k]]]]}], {k, {5, 20, 50}}] > Image[Clip[restack[50], {0., 255.}]/255.] > ListLogPlot[Table[Module[{d = SingularValueList[rgb[[All, All, c]]]}, d/First[d]], {c, 3}], > Joined -> True, PlotLegends -> {"R", "G", "B"}, AxesLabel -> {"i", "sigma ratio"}] > ``` > *Look for:* the same error numbers; `Transpose[..., {3,1,2}]` turns the list of three channel matrices back into an $H\times W\times 3$ array; `C` is a protected name, hence `chan`. What the color run should produce: ![Plush reconstructions at 1%, 4%, 10% of the singular values per channel; the energy percentages in the titles combine all three channels](Media/lowrank_rgb_panels.png) ![Channel spectra of the plush photo against the sprite: smooth decay, no cliff](Media/lowrank_rgb_spectra.png) > [!warning] Verify the color run (required) > - Channel $\sigma_1$: $124138.3$ (R), $119875.1$ (G), $109302.3$ (B); effective rank $494$ in every channel (no cliff). > - Average pixel error over all three channels: $24.72,\ 10.90,\ 5.44$ at $k=5,20,50$; restacking all $512$ layers per channel rebuilds the photo to $\sim10^{-11}$. > - Energy at $k=5$: $97.83\%$ (R), $97.43\%$ (G), $95.27\%$ (B), yet the $k=5$ panel is visibly blurry with color fringing. The large white background flatters the energy count; the eye lives in the discarded tail. Say so in your write-up. ### Your images: three spectral personalities Now run the study on three images you choose, picked in advance to have three different spectral personalities. The order of operations matters every time: decompose first, read the spectrum, and only then decide what to display. 1. **Image A, concentrated content.** Choose something you predict holds *little latent detail* and should compress beautifully: large smooth regions, strong geometry, little texture (a logo, a cartoon still, a flag, clean architecture). Prediction to test: fast decay, high energy at tiny $k$, a small $k$ already acceptable to the eye. 2. **Image B, channel personality.** Choose something whose red, green, and blue data should *genuinely differ*: dominated by saturated, distinct colors (a sunset, neon signs, stained glass, a bold poster). Report the three channel spectra separately: where does each channel reach $95\%$ energy, which channel is the expensive one, and at what $k$ does the color fringing die out? 3. **Image C, distributed content.** Choose something you predict carries *a lot of latent detail* and should resist compression: fine texture everywhere (foliage, gravel, a crowd, fur, fabric). Prediction to test: slow decay, heavy tail, many terms before the eye is satisfied. For each image, present the spectrum (per channel if RGB) and comment on it: decay shape, cliff or slide, effective rank, where it sits between our sprite and our plush, and whether a smooth background flatters the energy count. Then define your own interesting percentages *for that image*, with stated reasons (e.g., cheapest recognizable, the $95\%$ energy point, visually indistinguishable), and display them alongside $100\%$ with energy captured, average pixel error, and storage cost for each. Close with a cross-image comparison: rank the three by compressibility and defend the ranking from the spectra, not from the looks. > [!note] "Latent" is the right word > The layers $\sigma_i\mathbf{u}_i\mathbf{v}_i^{T}$ are the image's latent structure: content carried by whole rank-one patterns rather than any single pixel. An image with little latent detail concentrates its content in a few layers; one with a lot spreads it across hundreds. This is the same latent-factor language used by recommender systems (see Real-World Context). Some conversion helpers: > [!example]- Getting a grayscale matrix from your picture > ```python > from PIL import Image > Image.open("my_photo.jpg").convert("L").save("my_gray.png") # "L" = 8-bit grayscale > ``` > ```matlab > imwrite(rgb2gray(imread('my_photo.jpg')), 'my_gray.png') > ``` > ```r > library(png); library(jpeg) > rgb <- jpeg::readJPEG("my_photo.jpg") > gray <- 0.299*rgb[,,1] + 0.587*rgb[,,2] + 0.114*rgb[,,3] # standard luma weights > png::writePNG(gray, "my_gray.png") > ``` > ```wolfram > Export["my_gray.png", ColorConvert[Import["my_photo.jpg"], "Grayscale"]] > ``` > *Look for:* a single-channel matrix after loading (two dimensions, not three). #### Student task loop for implementation, analysis and reflection 1. **Predict** what the reconstruction should look like and roughly what energy it captures, before running. 2. **Implement** the truncation at that percentage. 3. **Compare** the energy captured and average pixel error against your prediction, and the visual quality against the panels. 4. **Interpret** what the result says about your image's redundancy and the compression tradeoff. ## Reflection Framework ### Geometry and structure - In the warm-up examples, verify the singular vectors are orthonormal. What do the three SVD factors (rotate, stretch, rotate) each do? - How fast do your image's singular values decay, and what does that decay say about its redundancy? ### Compression accounting - Why can $0.2\%$ of the singular values carry $79.1\%$ of the energy? Reconcile the two percentages, and explain the decay cliff at $k=83$ from the sprite's blocky geometry. - The plush photo reaches $95\%$ energy by $k=2$ in two of its three channels, yet the $k=5$ reconstruction is visibly blurry. What does a large smooth background do to the energy count, and why does the eye disagree? - What rank gives *your* image acceptable quality, and what storage ratio does that represent? ### Mathematical insights - **Energy captured:** relate the energy curve to the panels and the average pixel errors: where does additional energy stop buying visible quality, and why? - **[Condition number](https://en.wikipedia.org/wiki/Condition_number):** $\sigma_1/\sigma_r$ measures how close the matrix is to lower rank, and hence how sensitive it is. For an exactly rank-deficient matrix like our sprite, use the smallest *nonzero* singular value ($\sigma_1/\sigma_{83}$ here); dividing by a roundoff-level $\sigma_{512}$ says nothing. ## (Optional) Mathematical Extensions These go beyond the project's stated level; attempt them if interested. - **Optimality ([Eckart-Young-Mirsky](https://en.wikipedia.org/wiki/Low-rank_approximation)):** the truncated SVD is not just a good rank-$k$ approximation, it is provably the best possible one. Formalizing "best" requires matrix norms; read or reconstruct the argument. - **PCA connection:** center the pixel columns and relate the singular vectors to principal components. - **Randomized SVD:** implement a sketch-based top-$k$ SVD and compare accuracy and speed to the full decomposition. - **Transform-coding comparison:** contrast the SVD's data-dependent basis with the fixed [discrete cosine transform](https://en.wikipedia.org/wiki/Discrete_cosine_transform) basis used by JPEG. - **[Matrix completion](https://en.wikipedia.org/wiki/Matrix_completion):** recover missing entries by seeking a low-rank fit (the [nuclear norm](https://en.wikipedia.org/wiki/Matrix_norm#Schatten_norms) relaxation), which links this project to recommender systems. ## (Optional) Real-World Context Low-rank approximation via the SVD is one of the most widely used ideas in applied linear algebra. ### Applications - **[Image and data compression](https://en.wikipedia.org/wiki/Low-rank_approximation):** keep the dominant layers, discard the rest (the idea behind this project, and a cousin of transform coders like [JPEG](https://en.wikipedia.org/wiki/JPEG)). - **[Principal component analysis](https://en.wikipedia.org/wiki/Principal_component_analysis):** PCA is the SVD of centered data, used for dimensionality reduction across the sciences. - **[Latent semantic analysis](https://en.wikipedia.org/wiki/Latent_semantic_analysis):** truncated SVD of a term-document matrix uncovers topics in text. - **[Recommender systems](https://en.wikipedia.org/wiki/Collaborative_filtering):** the same latent-factor idea drives collaborative filtering (see the companion project on latent factors). ### Technical challenges - **[Computational cost](https://en.wikipedia.org/wiki/Singular_value_decomposition#Numerical_approach):** a full SVD costs about $O(mn\min(m,n))$, which is heavy for large images or matrices. - **[Randomized and truncated SVD](https://en.wikipedia.org/wiki/Randomized_algorithm):** when only the top $k$ layers are needed, randomized methods approximate them far more cheaply. - **[Numerical stability](https://en.wikipedia.org/wiki/Numerical_stability):** the SVD is computed stably, but forming $A^{T}A$ explicitly (the textbook route) squares the condition number and should be avoided in practice. ### Why it matters For the shipped $512\times512$ sprite, keeping $4\%$ of the singular values costs under $8\%$ of the storage yet captures $99.7\%$ of the energy, and what was lost is not a guess: it is exactly the sum of the discarded layers. The same arithmetic, at scale, is how latent-factor models compress user-item tables with millions of entries.