5  Linear Algebra of Images

Linear algebra is the first common language of modern computer vision. A grayscale image can be viewed as a matrix, a color image as a tensor, a collection of images as a data matrix, and a learned representation as a vector in a high-dimensional feature space. This chapter develops the clean mathematical foundations behind vectorization, inner products, norms, linear operators, projections, SVD, PCA, low-rank approximation, and feature embeddings. The audience is MA Applied Mathematics and MS Statistics students, so the emphasis is on precise objects, dimensions, assumptions, and statistical interpretation.

5.1 Learning Goals

After reading this chapter, students should be able to:

1. represent images as vectors, matrices, tensors, and points in high-dimensional Euclidean spaces;
2. define useful inner products, norms, distances, similarities, and projections for image data;
3. write image collections as data matrices and explain centering, scaling, covariance, and mean images;
4. interpret image processing operations as linear maps or affine maps;
5. use least squares and orthogonal projections to approximate images by lower-dimensional subspaces;
6. state and apply the singular value decomposition and the Eckart–Young theorem;
7. derive principal component analysis from variance maximization and reconstruction-error minimization;
8. explain why high-dimensional image data require care with conditioning, sample size, and statistical generalization;
9. connect classical linear representations with modern AI components such as patch embeddings, linear layers, low-rank adaptation, and linear probing.

5.2 4.1 Why Linear Algebra Comes Early in Computer Vision

The previous chapters described images as arrays, coordinate-based objects, and sampled functions. Linear algebra gives a different but equally important viewpoint:

\[ \text{image} \quad \longleftrightarrow \quad \text{vector in a high-dimensional space}. \]

If a grayscale image has height \(H\) and width \(W\), then it has \(p=HW\) pixels. After choosing an ordering of pixels, the image becomes a vector

\[ x = \operatorname{vec}(I) \in \mathbb{R}^{p}. \]

If the image has \(C\) channels, then \(p=HWC\), and

\[ x = \operatorname{vec}(I) \in \mathbb{R}^{HWC}. \]

This simple transformation is the beginning of many important ideas:

• a dataset of images becomes a point cloud in \(\mathbb{R}^p\);
• distances between images become Euclidean or cosine distances;
• denoising becomes projection onto a lower-dimensional subspace;
• compression becomes low-rank approximation;
• recognition becomes classification in a feature space;
• deep learning layers become compositions of affine maps and nonlinearities;
• vision transformers begin by mapping image patches into vectors.

The central idea of this chapter is this:

\[ \text{computer vision representation} = \text{choice of vector space} + \text{choice of coordinates} + \text{choice of metric}. \]

The same image can have very different mathematical behavior depending on how it is embedded into a vector space and how similarity is measured.

Linear algebra is not only a computational convenience. It determines what a model considers close, smooth, low-dimensional, redundant, stable, or informative.

5.3 4.2 Images as Vectors

Let

\[ I \in \mathbb{R}^{H \times W} \]

be a grayscale image. A vectorization map is a bijection between pixel indices \((i,j)\) and vector indices \(k\). One common row-major convention is

\[ k = iW + j, \qquad 0 \le i < H, \quad 0 \le j < W. \]

Then

\[ x_k = I[i,j] \quad \text{where} \quad k=iW+j. \]

Equivalently,

\[ x = \operatorname{vec}(I) = \begin{bmatrix} I[0,0] \\ I[0,1] \\ \vdots \\ I[0,W-1] \\ I[1,0] \\ \vdots \\ I[H-1,W-1] \end{bmatrix} \in \mathbb{R}^{HW}. \]

For color images,

\[ I \in \mathbb{R}^{H \times W \times C}. \]

A vectorization convention must specify whether the channel index changes fastest or slowest. For example,

\[ k = (iW+j)C+c \]

gives

\[ x_k = I[i,j,c]. \]

Vectorization is mathematically simple but convention-dependent. If two software libraries use different memory layouts, the same formula \(x=\operatorname{vec}(I)\) may correspond to different coordinate orderings. Always state the convention when dimensions matter.

5.3.1 Interactive demo: vectorization order

The figure below shows a small image and how its pixels are ordered into a vector. The important point is not the particular ordering, but that an ordering is a coordinate choice.

5.4 4.3 Inner Products and Norms for Images

Once an image is treated as a vector, we can use the Euclidean inner product. For two grayscale images

\[ I,J \in \mathbb{R}^{H \times W}, \]

the Frobenius inner product is

\[ \langle I,J \rangle_F = \sum_{i=0}^{H-1}\sum_{j=0}^{W-1} I[i,j]J[i,j]. \]

If \(x=\operatorname{vec}(I)\) and \(z=\operatorname{vec}(J)\), then

\[ \langle I,J \rangle_F = x^Tz. \]

The associated Frobenius norm is

\[ \|I\|_F = \sqrt{\langle I,I\rangle_F} = \left(\sum_{i,j} I[i,j]^2\right)^{1/2}. \]

For color images,

\[ \langle I,J\rangle_F = \sum_{i,j,c} I[i,j,c]J[i,j,c], \qquad \|I\|_F^2 = \sum_{i,j,c} I[i,j,c]^2. \]

5.4.1 Distances between images

The Euclidean image distance is

\[ d_2(I,J) = \|I-J\|_F. \]

The mean squared error is

\[ \operatorname{MSE}(I,J) = \frac{1}{HWC}\|I-J\|_F^2. \]

The peak signal-to-noise ratio, often used in image reconstruction, is

\[ \operatorname{PSNR}(I,J) = 10\log_{10}\left(\frac{M^2}{\operatorname{MSE}(I,J)}\right), \]

where \(M\) is the maximum possible intensity value, such as \(1\) for normalized images or \(255\) for eight-bit images.

5.4.2 Cosine similarity

In representation learning, magnitude may be less important than direction. The cosine similarity between nonzero vectors \(x,z \in \mathbb{R}^p\) is

\[ \operatorname{cos}(x,z) = \frac{x^Tz}{\|x\|_2\|z\|_2}. \]

If \(x\) and \(z\) are feature vectors from a neural network, cosine similarity measures angular closeness in feature space. This becomes important for contrastive learning and vision-language models later in the book.

From a statistical viewpoint, a distance is an assumption about noise and variation. Squared Euclidean distance corresponds naturally to independent Gaussian noise with equal variance in each coordinate. If image coordinates have different noise levels or correlations, then a Mahalanobis distance may be more appropriate.

The Mahalanobis distance associated with a positive definite covariance matrix \(\Sigma\) is

\[ d_{\Sigma}(x,z)^2 = (x-z)^T\Sigma^{-1}(x-z). \]

This distance downweights directions of high variance and upweights directions of low variance. It is mathematically natural but statistically difficult in images because \(p=HWC\) can be extremely large.

5.5 4.4 Image Datasets as Data Matrices

Suppose we have \(n\) images, each vectorized into \(x_i \in \mathbb{R}^{p}\). The image dataset can be written as a data matrix

\[ X = \begin{bmatrix} - & x_1^T & - \\ - & x_2^T & - \\ & \vdots & \\ - & x_n^T & - \end{bmatrix} \in \mathbb{R}^{n \times p}. \]

Each row is one image. Each column is one pixel or channel coordinate.

The sample mean image is

\[ \bar{x} = \frac{1}{n}\sum_{i=1}^n x_i \in \mathbb{R}^{p}. \]

The centered data matrix is

\[ X_c = X - \mathbf{1}\bar{x}^T, \]

where \(\mathbf{1}\in \mathbb{R}^n\) is the vector of all ones.

The empirical covariance matrix is

\[ S = \frac{1}{n-1}X_c^TX_c \in \mathbb{R}^{p \times p}. \]

This object is central in PCA, denoising, and statistical modeling.

5.5.1 High-dimensional issue

For modern images, \(p\) is often much larger than \(n\). For example, \(224\times 224\times 3\) images have

\[ p = 224\cdot 224\cdot 3 = 150{,}528 \]

coordinates. If \(n=5000\), then the covariance matrix has size

\[ 150{,}528 \times 150{,}528, \]

which is far too large to store directly in ordinary memory.

In image statistics, \(p \gg n\) is common. The sample covariance matrix is then singular because

\[ \operatorname{rank}(S) \le \min(n-1,p). \]

This is one reason why dimensionality reduction, regularization, and learned representations are essential in computer vision.

5.6 4.5 Linear Maps on Images

A map \(T:\mathbb{R}^p\to\mathbb{R}^q\) is linear if

\[ T(ax+bz)=aT(x)+bT(z) \]

for all vectors \(x,z\) and scalars \(a,b\). Every finite-dimensional linear map can be represented by a matrix

\[ T(x)=Ax, \qquad A\in\mathbb{R}^{q\times p}. \]

Many image operations are linear or approximately linear:

Operation	Linear?	Mathematical form
brightness scaling	yes	\(x \mapsto ax\)
brightness shift	affine, not linear	\(x \mapsto x+b\mathbf{1}\)
grayscale conversion	yes	channel matrix multiplication
blur/filtering	yes	sparse structured matrix
downsampling after averaging	yes	\(x \mapsto Ax\)
rotation with interpolation	usually linear in intensities	\(x \mapsto Ax\)
gamma correction	no	\(x \mapsto x^\gamma\) coordinatewise
ReLU	no	\(x \mapsto \max(x,0)\)
softmax	no	normalized exponential map

5.6.1 Color conversion as a linear map

A simple grayscale conversion has the form

\[ Y = 0.299R + 0.587G + 0.114B. \]

For each pixel,

\[ \begin{bmatrix}Y\end{bmatrix} = \begin{bmatrix} 0.299 & 0.587 & 0.114 \end{bmatrix} \begin{bmatrix} R\\G\\B \end{bmatrix}. \]

For an entire image, the same \(1\times 3\) matrix is applied at every pixel.

5.6.2 Filtering as a matrix

A one-dimensional smoothing filter with kernel

\[ h = \begin{bmatrix}1/4 & 1/2 & 1/4\end{bmatrix} \]

maps a signal \(x\in\mathbb{R}^m\) to \(y\in\mathbb{R}^m\) by

\[ y_i = \frac14x_{i-1}+\frac12x_i+\frac14x_{i+1}, \]

with a boundary convention. This can be written as

\[ y = Ax, \]

where \(A\) is a banded matrix. In two dimensions, blurs and finite-difference filters can also be written as large sparse matrices, although in practice they are implemented without explicitly forming the matrix.

5.7 4.6 Kronecker Products and Separable Image Operators

Many image operators have tensor-product structure. If

\[ I \in \mathbb{R}^{H\times W}, \]

then a row operation can be written as

\[ I \mapsto IA^T, \qquad A\in\mathbb{R}^{W'\times W}, \]

and a column operation as

\[ I \mapsto BI, \qquad B\in\mathbb{R}^{H'\times H}. \]

Applying both gives

\[ I \mapsto BIA^T. \]

Under vectorization,

\[ \operatorname{vec}(BIA^T) = (A\otimes B)\operatorname{vec}(I), \]

where \(\otimes\) denotes the Kronecker product.

The identity

\[ \operatorname{vec}(BIA^T) = (A\otimes B)\operatorname{vec}(I) \]

is one of the most useful linear algebra identities for image processing. It explains why many huge image operators can be represented compactly by smaller one-dimensional operators.

A two-dimensional Gaussian blur is often approximately separable: first blur horizontally, then vertically. Instead of storing a huge \(HW\times HW\) matrix, one stores two much smaller one-dimensional kernels.

5.8 4.7 Orthogonal Projection and Least Squares

A central problem in representation is approximating an image vector \(x\in\mathbb{R}^p\) using a smaller set of basis images.

Let

\[ U = \begin{bmatrix}u_1 & u_2 & \cdots & u_k\end{bmatrix} \in \mathbb{R}^{p\times k} \]

have orthonormal columns:

\[ U^TU=I_k. \]

The orthogonal projection of \(x\) onto the subspace

\[ \mathcal{U}=\operatorname{span}\{u_1,\ldots,u_k\} \]

is

\[ \widehat{x}=UU^Tx. \]

The coefficient vector is

\[ a = U^Tx \in \mathbb{R}^k, \]

and the reconstruction is

\[ \widehat{x}=Ua. \]

The projection solves the least-squares problem

\[ \widehat{x} = \arg\min_{z\in\mathcal{U}} \|x-z\|_2^2. \]

Equivalently,

\[ a^* = \arg\min_{a\in\mathbb{R}^k}\|x-Ua\|_2^2. \]

If \(U\) has orthonormal columns, then \(a^*=U^Tx\).

If \(B\in\mathbb{R}^{p\times k}\) is a general basis matrix with full column rank, the least-squares coefficients are

\[ a^* = (B^TB)^{-1}B^Tx. \]

The projection matrix is

\[ P_B = B(B^TB)^{-1}B^T. \]

Projection is both a geometric and statistical idea. Geometrically, it gives the closest point in a subspace. Statistically, it can remove noise if signal lies mostly in a low-dimensional subspace and noise spreads across many directions.

5.9 4.8 Singular Value Decomposition of an Image Matrix

For a single grayscale image

\[ I\in\mathbb{R}^{H\times W}, \]

the singular value decomposition is

\[ I = U\Sigma V^T, \]

where

\[ U\in\mathbb{R}^{H\times H}, \qquad V\in\mathbb{R}^{W\times W}, \]

are orthogonal matrices and \(\Sigma\) contains nonnegative singular values

\[ \sigma_1\ge \sigma_2\ge \cdots \ge 0. \]

The rank-\(r\) truncated SVD is

\[ I_r = \sum_{\ell=1}^r \sigma_\ell u_\ell v_\ell^T. \]

This gives a low-rank approximation of the image.

5.9.1 Eckart–Young theorem

The rank-\(r\) truncated SVD is optimal in Frobenius norm:

\[ I_r = \arg\min_{\operatorname{rank}(B)\le r}\|I-B\|_F. \]

Moreover,

\[ \|I-I_r\|_F^2 = \sum_{\ell>r}\sigma_\ell^2. \]

This theorem explains why SVD is the mathematical foundation for one classical form of image compression.

Low rank does not mean semantic understanding. A low-rank approximation may preserve broad smooth structure but lose small objects, edges, text, and fine textures. Modern vision models often need features that are not captured by simple pixel-space low-rank structure.

5.9.2 Interactive demo: low-rank image approximation

The next demo shows how increasing the rank improves reconstruction. Smooth structure often appears early; fine details require more singular components.

5.10 4.9 SVD of a Dataset and Principal Component Analysis

For a dataset matrix

\[ X_c\in\mathbb{R}^{n\times p}, \]

the singular value decomposition is

\[ X_c = U\Sigma V^T. \]

The right singular vectors

\[ v_1,\ldots,v_r \in \mathbb{R}^p \]

are principal directions in image space. The principal component scores are

\[ Z = X_cV_k \in \mathbb{R}^{n\times k}, \]

where

\[ V_k = \begin{bmatrix}v_1 & \cdots & v_k\end{bmatrix}. \]

Each row of \(Z\) gives the low-dimensional coordinates of an image.

The empirical covariance matrix is

\[ S = \frac{1}{n-1}X_c^TX_c. \]

Since

\[ X_c^TX_c = V\Sigma^2V^T, \]

the principal directions are eigenvectors of \(S\), and the eigenvalues are

\[ \lambda_j = \frac{\sigma_j^2}{n-1}. \]

5.10.1 PCA as variance maximization

The first principal direction solves

\[ v_1 = \arg\max_{\|v\|_2=1} v^TSv. \]

The quantity

\[ v^TSv \]

is the sample variance of the projected data \(X_cv\).

The next directions solve the same problem subject to orthogonality constraints:

\[ v_j = \arg\max_{\|v\|_2=1,\ v\perp v_1,\ldots,v_{j-1}} v^TSv. \]

5.10.2 PCA as reconstruction-error minimization

PCA also solves a reconstruction problem. Among all \(k\)-dimensional linear subspaces, PCA finds the subspace minimizing total squared reconstruction error:

\[ \min_{\dim(\mathcal{U})=k}\sum_{i=1}^n \left\|x_i-\bar{x}-P_{\mathcal{U}}(x_i-\bar{x})\right\|_2^2. \]

The solution is

\[ \mathcal{U}=\operatorname{span}\{v_1,\ldots,v_k\}. \]

PCA has two equivalent interpretations:

\[ \text{maximize variance of projections} \quad \Longleftrightarrow \quad \text{minimize squared reconstruction error}. \]

For image data, these interpretations explain why PCA can be used both for visualization and compression.

5.10.3 Interactive demo: PCA of synthetic images

The next demo creates a small synthetic image dataset with two sources of variation: horizontal location and vertical location of a bright blob. PCA identifies the dominant directions of variation.

5.11 4.10 Dual PCA When Pixels Outnumber Images

When \(p\) is very large, directly forming

\[ S = \frac{1}{n-1}X_c^TX_c\in\mathbb{R}^{p\times p} \]

may be impossible. If \(n\ll p\), one can instead work with the smaller Gram matrix

\[ G = X_cX_c^T \in \mathbb{R}^{n\times n}. \]

Suppose

\[ G u_j = \sigma_j^2 u_j. \]

Then the corresponding right singular vector is

\[ v_j = \frac{1}{\sigma_j}X_c^T u_j. \]

This is the algebraic idea behind classical eigenimage and eigenface methods. Instead of diagonalizing a huge pixel covariance matrix, one diagonalizes an \(n\times n\) matrix.

Dual PCA saves computation but does not solve the statistical problem automatically. If the dataset is small, the leading components may describe lighting, alignment, background, or acquisition artifacts rather than semantic content.

5.12 4.11 Low-Dimensional Subspaces, Manifolds, and Image Variation

PCA assumes that important variation is approximately linear. This is useful but limited.

Consider a simple image containing one object translated horizontally. The set of all translated images is generally not a linear subspace. If \(x(t)\) is the vectorized image at horizontal shift \(t\), then

\[ \mathcal{M} = \{x(t): t\in\mathbb{R}\} \]

is closer to a one-dimensional curved manifold in \(\mathbb{R}^p\) than to a one-dimensional linear subspace.

For small shifts, however, a first-order approximation gives

\[ x(t+\Delta t) \approx x(t) + \Delta t\,x'(t), \]

so locally the manifold is approximately linear. This explains why tangent spaces and local linearization appear repeatedly in vision geometry, optical flow, and representation learning.

Linear methods are often good local approximations. Nonlinear vision models, including deep networks, can be viewed as learning representations in which complicated visual manifolds become easier to separate or approximate linearly.

5.13 4.12 Tensors, Channels, and Feature Maps

A modern vision model rarely operates only on vectors. It usually operates on tensors.

A batch of images may be stored as

\[ X\in\mathbb{R}^{N\times H\times W\times C} \]

or

\[ X\in\mathbb{R}^{N\times C\times H\times W}, \]

depending on the software convention.

A convolutional layer transforms an image tensor into a feature tensor

\[ F \in \mathbb{R}^{H'\times W'\times D}, \]

where \(D\) is the number of learned channels. At each spatial location, the feature vector

\[ F[i,j,:]\in\mathbb{R}^D \]

summarizes local visual information.

5.13.1 Tensor flattening for neural networks

Many neural networks eventually flatten or pool feature maps into a vector

\[ z\in\mathbb{R}^d. \]

A classifier then applies an affine map

\[ s = Wz+b, \qquad W\in\mathbb{R}^{K\times d}, \quad b\in\mathbb{R}^{K}, \]

where \(K\) is the number of classes. The output \(s\) contains logits.

A softmax layer transforms logits into class probabilities:

\[ \widehat{p}_k = \frac{e^{s_k}}{\sum_{\ell=1}^K e^{s_\ell}}. \]

The affine part is linear algebra. The full classifier is nonlinear because of softmax and because \(z\) is usually produced by many nonlinear layers.

5.14 4.13 Patch Embeddings and Vision Transformers

Vision transformers begin by dividing an image into patches. Suppose the image has shape

\[ H\times W\times C, \]

and each patch has size

\[ P\times P\times C. \]

Each patch is vectorized into

\[ p_i \in \mathbb{R}^{P^2C}. \]

A learned patch embedding is an affine map

\[ z_i = Ep_i + b, \qquad E\in\mathbb{R}^{d\times P^2C}, \quad b\in\mathbb{R}^{d}. \]

Thus, before attention appears, a vision transformer already begins with a linear algebra step:

\[ \text{patch} \quad \mapsto \quad \text{token vector}. \]

The number of patches is

\[ N_p = \frac{H}{P}\cdot \frac{W}{P} \]

when \(P\) divides both \(H\) and \(W\).

Example. A \(224\times224\) RGB image with patch size \(16\times16\) has

\[ N_p = 14\cdot14 = 196 \]

patches. Each patch has dimension

\[ 16\cdot16\cdot3 = 768. \]

A patch embedding maps each vector in \(\mathbb{R}^{768}\) to a token in \(\mathbb{R}^d\).

This connects classical vectorization with modern foundation models. The difference is that the vectorization is now local, patchwise, and followed by learned attention-based processing.

5.15 4.14 Linear Probing and Feature Embeddings

A trained vision model often produces an embedding

\[ \phi_\theta(I) = z \in \mathbb{R}^d. \]

A linear probe freezes \(\phi_\theta\) and trains only a linear classifier:

\[ \widehat{y} = \arg\max_{k} (Wz+b)_k. \]

The purpose is to test whether the representation already makes the classes approximately linearly separable.

A good representation is not merely a smaller vector. It is a coordinate system in which important tasks become simpler. Linear probing asks whether a difficult visual classification problem has become approximately linear after representation learning.

This idea will return in later chapters on representation learning, self-supervised learning, CLIP-style models, and foundation models.

5.16 4.15 Statistical Cautions in High-Dimensional Image Spaces

Linear algebra gives powerful tools, but raw pixel space is statistically challenging.

5.16.1 Curse of dimensionality

If \(x\in\mathbb{R}^p\) and \(p\) is large, distances can become less informative. For many high-dimensional distributions, pairwise distances concentrate. This can make nearest-neighbor reasoning in raw pixel space unreliable.

5.16.2 Nuisance variation

Large pixel-space variation may be caused by nuisance factors:

• lighting;
• camera exposure;
• background;
• pose;
• small translations;
• sensor noise;
• compression artifacts.

PCA may faithfully identify high-variance directions, but high variance is not the same as semantic importance.

5.16.3 Rank and sample size

If \(p\gg n\), the data matrix \(X_c\in\mathbb{R}^{n\times p}\) has rank at most \(n-1\). Therefore, any covariance estimate based on raw pixels is low rank and unstable without additional assumptions.

5.16.4 Bias from preprocessing

Cropping, resizing, normalization, augmentation, and alignment all change the geometry of the data. These steps are not neutral. They change the vector space in which linear algebra is being applied.

5.17 4.16 Python Example: Vectorization, Inner Products, and MSE

import numpy as np

H, W = 32, 32
rng = np.random.default_rng(7243)

# Synthetic grayscale image: smooth blob plus noise
x_grid = np.linspace(-1, 1, W)
y_grid = np.linspace(-1, 1, H)
Xg, Yg = np.meshgrid(x_grid, y_grid)
I = np.exp(-6 * (Xg**2 + Yg**2))
J = I + 0.08 * rng.normal(size=(H, W))

# Vectorization
x = I.reshape(-1)
z = J.reshape(-1)

fro_inner = np.sum(I * J)
vec_inner = x @ z

mse = np.mean((I - J)**2)
fro_error = np.linalg.norm(I - J, ord="fro")

print("Frobenius inner product:", fro_inner)
print("Vector inner product:", vec_inner)
print("MSE:", mse)
print("Frobenius reconstruction error:", fro_error)

This code illustrates the identity

\[ \langle I,J\rangle_F = \operatorname{vec}(I)^T\operatorname{vec}(J). \]

5.18 4.17 Python Example: SVD Image Compression

import numpy as np
import matplotlib.pyplot as plt

H, W = 96, 96
x = np.linspace(-1, 1, W)
y = np.linspace(-1, 1, H)
X, Y = np.meshgrid(x, y)

# Synthetic image with smooth and sharp structure
I = (
    np.exp(-8 * ((X + 0.35)**2 + (Y - 0.2)**2))
    + 0.8 * np.exp(-14 * ((X - 0.35)**2 + (Y + 0.25)**2))
    + 0.25 * (X > 0).astype(float)
)
I = (I - I.min()) / (I.max() - I.min())

U, s, Vt = np.linalg.svd(I, full_matrices=False)

ranks = [1, 3, 5, 10, 20, 40]
plt.figure(figsize=(10, 6))
for idx, r in enumerate(ranks, start=1):
    Ir = (U[:, :r] * s[:r]) @ Vt[:r, :]
    error = np.linalg.norm(I - Ir, ord="fro") / np.linalg.norm(I, ord="fro")
    plt.subplot(2, 3, idx)
    plt.imshow(Ir, cmap="gray", vmin=0, vmax=1)
    plt.title(f"rank {r}, rel err={error:.3f}")
    plt.axis("off")
plt.tight_layout()
plt.show()

energy = np.cumsum(s**2) / np.sum(s**2)
plt.figure(figsize=(6, 4))
plt.plot(np.arange(1, len(s) + 1), energy)
plt.xlabel("rank r")
plt.ylabel("cumulative singular-value energy")
plt.title("Energy captured by rank-r approximation")
plt.grid(True)
plt.show()

5.19 4.18 Python Example: PCA on an Image Dataset

import numpy as np
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

rng = np.random.default_rng(7243)
H, W = 32, 32
n = 300

xs = np.linspace(-1, 1, W)
ys = np.linspace(-1, 1, H)
Xg, Yg = np.meshgrid(xs, ys)

images = []
labels = []
for i in range(n):
    cx = rng.uniform(-0.5, 0.5)
    cy = rng.uniform(-0.5, 0.5)
    sigma = rng.uniform(0.12, 0.20)
    img = np.exp(-((Xg - cx)**2 + (Yg - cy)**2) / (2 * sigma**2))
    img += 0.05 * rng.normal(size=(H, W))
    images.append(img)
    labels.append([cx, cy])

X_data = np.array([img.reshape(-1) for img in images])
labels = np.array(labels)

pca = PCA(n_components=10)
Z = pca.fit_transform(X_data)

print("Data matrix shape:", X_data.shape)
print("Scores shape:", Z.shape)
print("Explained variance ratio:", pca.explained_variance_ratio_[:5])

plt.figure(figsize=(6, 5))
plt.scatter(Z[:, 0], Z[:, 1], c=labels[:, 0], s=20)
plt.xlabel("PC 1 score")
plt.ylabel("PC 2 score")
plt.title("PCA scores colored by true horizontal position")
plt.colorbar(label="horizontal center")
plt.show()

fig, axes = plt.subplots(2, 5, figsize=(10, 4))
for ax, comp in zip(axes.ravel(), pca.components_[:10]):
    ax.imshow(comp.reshape(H, W), cmap="gray")
    ax.axis("off")
plt.suptitle("First 10 principal directions as images")
plt.tight_layout()
plt.show()

5.20 4.19 TensorFlow Example: Patch Embedding as a Linear Layer

import tensorflow as tf

batch_size = 8
H, W, C = 224, 224, 3
patch_size = 16
embed_dim = 128

images = tf.random.normal((batch_size, H, W, C))

# Extract non-overlapping patches.
patches = tf.image.extract_patches(
    images=images,
    sizes=[1, patch_size, patch_size, 1],
    strides=[1, patch_size, patch_size, 1],
    rates=[1, 1, 1, 1],
    padding="VALID",
)

# Shape: batch_size x 14 x 14 x (16*16*3)
print("patch tensor shape:", patches.shape)

num_patches = (H // patch_size) * (W // patch_size)
patch_dim = patch_size * patch_size * C

patch_vectors = tf.reshape(patches, (batch_size, num_patches, patch_dim))
print("patch vectors shape:", patch_vectors.shape)

# A learned linear map from patch_dim to embed_dim.
projection = tf.keras.layers.Dense(embed_dim)
tokens = projection(patch_vectors)
print("token shape:", tokens.shape)

Mathematically, this implements

\[ z_i = Ep_i+b \]

for every image patch \(p_i\in\mathbb{R}^{P^2C}\).

5.21 4.20 AI-Assisted Study Components

Modern AI tools can help students study linear algebra in computer vision, but they should be used with mathematical checking. The main danger is that an AI assistant may give a correct-looking formula with mismatched dimensions.

5.21.1 Dimension-checking prompt

Use the following prompt after deriving a formula:

I am studying linear algebra for computer vision. Check only the dimensions in the following derivation. Let \(X\in\mathbb{R}^{n\times p}\) be a centered image data matrix. Let \(V_k\in\mathbb{R}^{p\times k}\) contain the first \(k\) principal directions. I wrote \(Z=V_kX\). Is this dimensionally correct? If not, explain the correction without changing the mathematical goal.

Expected correction:

\[ Z = XV_k \in \mathbb{R}^{n\times k}. \]

5.21.2 Concept-checking prompt

Explain the difference between low-rank approximation of one image matrix \(I\in\mathbb{R}^{H\times W}\) and PCA of a dataset matrix \(X\in\mathbb{R}^{n\times p}\). Use formulas and state what the singular vectors represent in each case.

A good answer should say that single-image SVD decomposes spatial row/column structure inside one matrix, while dataset PCA decomposes variation across examples in a high-dimensional pixel space.

5.21.3 Modeling prompt

Suppose a PCA model on face images mainly separates images by brightness instead of identity. Give a mathematical explanation using variance maximization, and propose two preprocessing steps that could reduce this effect.

This encourages students to connect the PCA objective

\[ \max_{\|v\|=1}v^TSv \]

with nuisance variation.

5.22 4.21 Common Mistakes

1. Confusing image SVD with dataset PCA.
   SVD of one image \(I\in\mathbb{R}^{H\times W}\) studies row-column structure. PCA of \(X\in\mathbb{R}^{n\times p}\) studies variation across images.

2. Forgetting centering in PCA.
   PCA should usually be applied to \(X_c=X-\mathbf{1}\bar{x}^T\), not to raw \(X\).

3. Assuming Euclidean distance matches perception.
   Small pixel shifts can create large Euclidean distances even when images look similar to humans.

4. Ignoring dimensions.
   Many formulas fail because \(X\), \(X^T\), \(V\), and \(U\) are multiplied in the wrong order.

5. Interpreting high variance as semantic importance.
   PCA finds directions of large variance, not necessarily directions relevant for classification or understanding.

6. Forming huge matrices explicitly.
   Image operators are often too large to store. Use sparse, convolutional, separable, or implicit representations.

5.23 4.22 Summary

This chapter developed the linear algebra of images. The main ideas are:

• images can be vectorized into high-dimensional vectors;
• image datasets can be represented as data matrices;
• inner products, norms, and distances define a geometry on image space;
• many image operations are linear maps or affine maps;
• projection and least squares provide a foundation for approximation and denoising;
• SVD gives optimal low-rank approximation for a single image matrix;
• PCA gives optimal linear dimension reduction for centered datasets;
• high-dimensional image statistics require care because \(p\) is often much larger than \(n\);
• modern AI vision models still rely heavily on linear algebra through patch embeddings, affine layers, feature vectors, and linear probes.

The next chapter moves from finite-dimensional linear algebra to frequency analysis. Fourier methods reveal how image structure can be described by oscillations, filters, and spectral content.

5.24 Exercises

5.24.1 Conceptual exercises

1. Explain why an image can be represented both as a matrix and as a vector. What information is lost if the vectorization convention is not recorded?

2. Give two examples of image transformations that are linear and two that are nonlinear.

3. Explain why Euclidean distance between raw pixels may fail to represent perceptual similarity.

4. What is the difference between a basis image and a feature vector?

5. Why is centering important before PCA?

6. Explain why PCA may learn lighting variation instead of object identity.

5.24.2 Mathematical exercises

1. Let \(I,J\in\mathbb{R}^{H\times W}\). Prove that

   \[ \langle I,J\rangle_F = \operatorname{vec}(I)^T\operatorname{vec}(J). \]

2. Let \(U\in\mathbb{R}^{p\times k}\) have orthonormal columns. Prove that \(P=UU^T\) is an orthogonal projection matrix by showing that

   \[ P^T=P \qquad \text{and} \qquad P^2=P. \]

3. Let \(B\in\mathbb{R}^{p\times k}\) have full column rank. Derive the least-squares solution

   \[ a^*=(B^TB)^{-1}B^Tx \]

   for minimizing \(\|x-Ba\|_2^2\).

4. Suppose \(X_c=U\Sigma V^T\). Show that the covariance matrix

   \[ S=\frac{1}{n-1}X_c^TX_c \]

   has eigenvectors \(v_j\) and eigenvalues \(\sigma_j^2/(n-1)\).

5. Prove that the reconstruction error for PCA using \(k\) components is

   \[ \sum_{j>k}\sigma_j^2 \]

   for the centered data matrix, up to the normalization convention.

6. Verify the Kronecker identity

   \[ \operatorname{vec}(BIA^T)=(A\otimes B)\operatorname{vec}(I) \]

   for a \(2\times2\) image using explicit symbols.

5.24.3 Python exercises

1. Generate a synthetic image and compute its vectorization. Verify numerically that the Frobenius inner product equals the vector dot product.

2. Implement rank-\(r\) SVD reconstructions for a grayscale image. Plot reconstruction error versus \(r\).

3. Create a dataset of translated Gaussian blobs. Apply PCA and interpret the first two principal directions.

4. Compare Euclidean distances between an image and shifted copies of the same image. What happens as the shift increases?

5. Implement a simple linear classifier on PCA features for a synthetic two-class image dataset.

6. Use TensorFlow or PyTorch to extract image patches and apply a learned linear embedding to each patch.

5.24.4 AI-assisted exercises

1. Ask an AI assistant to explain the difference between PCA and SVD. Then identify one missing assumption or ambiguity in the answer.

2. Give an AI assistant a deliberately incorrect PCA formula, such as \(Z=V_kX\), and ask it to check dimensions.

3. Ask an AI assistant to design a preprocessing pipeline before PCA for face images. Evaluate whether the steps reduce nuisance variation.

4. Ask for a proof of the projection formula \(UU^Tx\). Check whether the proof uses the orthonormality assumption correctly.

5. Ask an AI assistant whether low-rank approximation is enough for object recognition. Write a short critique using examples from this chapter.

5.25 Notes and References

The linear algebra viewpoint is foundational throughout computer vision and statistical learning. Classical image-processing and computer-vision references include Gonzalez and Woods (2018) and Szeliski (2022). PCA, covariance methods, and statistical pattern recognition connect to Bishop (2006). Deep-learning uses of affine maps, learned representations, and neural network layers connect to Goodfellow et al. (2016). Later chapters return to patch embeddings and transformers, connecting this chapter with Dosovitskiy et al. (2020).

Bishop, Christopher M. 2006. Pattern Recognition and Machine Learning. Springer.

Dosovitskiy, Alexey, Lucas Beyer, Alexander Kolesnikov, et al. 2020. “An Image Is Worth 16x16 Words: Transformers for Image Recognition at Scale.” arXiv Preprint arXiv:2010.11929. https://arxiv.org/abs/2010.11929.

Gonzalez, Rafael C., and Richard E. Woods. 2018. Digital Image Processing. 4th ed. Pearson.

Goodfellow, Ian, Yoshua Bengio, and Aaron Courville. 2016. Deep Learning. MIT Press. https://www.deeplearningbook.org/.

Szeliski, Richard. 2022. Computer Vision: Algorithms and Applications. 2nd ed. Springer. https://szeliski.org/Book/.