## Introduction A digital image is a grid of pixels, where each pixel represents color or brightness at a specific location. In grayscale images, each pixel is stored as a single number—typically between 0 (black) and 255 (white)—forming a matrix of intensity values. For color images, three separate matrices represent the red, green, and blue (RGB) channels. This matrix structure allows us to apply linear algebra operations, like transformations and decompositions, directly to image data. Additionally, what linear algebra tells us about the data's internal structure can give us a means to form approximations to the data with less memory cost. ![[Teaching/MATH307/Projects/Media/link.gif]] ## Project Description Read the following mathematics on the fundamental theorem of linear algebra, and using it and a suitable computational environment, create low-rank approximations to a greyscale image of your choice. (1 point) Extend this result to RGB image with an animation showing the effect of adding additional terms to the approximation. (2 points) **Key Note (Deliverable)**: I am currently working on a format for your submission that is efficient for both you/us and also demonstrates that something was learned about the code, context, and mathematics surrounding this problem. ## Eigenvalues/Eigenvectors and Diagonalization If $A$ is a square matrix, an **eigenvector** $\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: $\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$. For example, $A = \begin{bmatrix}0 & 1 \\ 1 & 0 \end{bmatrix}\implies \lambda_{1} = 1, \mathbf{x}_{1}= \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \, \, \lambda_{2} = -1, \mathbf{x}_{1}= \begin{bmatrix} 1 \\ -1 \end{bmatrix}$ and one can check that if $ P=\begin{bmatrix}1 & 1 \\[6pt] 1 & -1\end{bmatrix},\quad D=\begin{bmatrix}1 & 0 \\[3pt] 0 & -1\end{bmatrix}, \quad P^{-1}=\tfrac12\begin{bmatrix}1 & 1 \\[3pt] 1 & -1\end{bmatrix},$ then $P\,D\,P^{-1}= \begin{bmatrix}1 & 1 \\[6pt] 1 & -1\end{bmatrix} \begin{bmatrix}1 & 0 \\[3pt] 0 & -1\end{bmatrix} \tfrac12\begin{bmatrix}1 & 1 \\[3pt] 1 & -1\end{bmatrix} = \begin{bmatrix}0 & 1 \\[3pt] 1 & 0\end{bmatrix} = A.$ ## Orthogonal Diagonalization of Symmetric Matrices The **matrix transpose** is formed by flipping rows and columns. If $ A = \begin{bmatrix} a & b \\ c & d \end{bmatrix}, \quad \text{then} \quad A^T = \begin{bmatrix} a & c \\ b & d \end{bmatrix} $ If the transpose of a matrix is equal to itself, then we say that matrix is **symmetric**. Linear algebra tells us that a symmetric matrix **always** has enough eigenvectors to perform a diagonal decomposition of a square matrix and that the eigenvectors form an **orthonormal set**, which allows us to write the inverse of the eigenvector matrix through transposition. For example, $ A = \begin{bmatrix}0 & 1 \\ 1 & 0 \end{bmatrix} \quad\Rightarrow\quad \lambda_{1} = 1, \ \mathbf{x}_{1}= \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad \lambda_{2} = -1, \ \mathbf{x}_{2}= \begin{bmatrix} 1 \\ -1 \end{bmatrix} $ We can check that these eigenvectors are orthogonal: $ \mathbf{x}_1 \cdot \mathbf{x}_2 = \begin{bmatrix} 1 \\ 1 \end{bmatrix} \cdot \begin{bmatrix} 1 \\ -1 \end{bmatrix} = 1 \cdot 1 + 1 \cdot (-1) = 0 $ To make them orthonormal, normalize each vector: $ \hat{\mathbf{x}}_1 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad \hat{\mathbf{x}}_2 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \end{bmatrix} $ Now define $ Q = \begin{bmatrix} \hat{\mathbf{x}}_1 & \hat{\mathbf{x}}_2 \end{bmatrix} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}, \quad D = \begin{bmatrix} 1 & 0 \\ 0 & -1 \end{bmatrix} $ Since $Q$ is orthogonal (its columns are orthonormal), we have $ A = Q D Q^T $ This is the **orthogonal diagonalization** of the symmetric matrix $A$. ## Spectral representation of a symmetric matrix The **spectral representation** expresses a symmetric matrix as a weighted sum of outer products, i.e., column-vector left multiplied onto a row-vector, of its orthonormal eigenvectors. Specifically, if $A$ is symmetric with eigenvalues $\lambda_i$ and orthonormal eigenvectors $\mathbf{x}_i$, then $ A = \sum_{i} \lambda_i \, \mathbf{x}_i \mathbf{x}_i^T $ This form shows how $A$ acts by scaling along its principal directions. Applying this to $ A = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} $ with eigenvalues $\lambda_1 = 1$, $\lambda_2 = -1$, and orthonormal eigenvectors: $ \mathbf{x}_1 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \end{bmatrix}, \quad\mathbf{x}_2 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \end{bmatrix} $ Compute each outer product: $ \mathbf{x}_1 \mathbf{x}_1^T = \frac{1}{2} \begin{bmatrix} 1 \\ 1 \end{bmatrix} \begin{bmatrix} 1 & 1 \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} $ $ \mathbf{x}_2 \mathbf{x}_2^T = \frac{1}{2} \begin{bmatrix} 1 \\ -1 \end{bmatrix} \begin{bmatrix} 1 & -1 \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} $ Now assemble the spectral representation: $ A = \lambda_1 \mathbf{x}_1 \mathbf{x}_1^T + \lambda_2 \mathbf{x}_2 \mathbf{x}_2^T = \left(1\right) \cdot \frac{1}{2} \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix}+(-1) \cdot \frac{1}{2} \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} $ $ = \frac{1}{2} \left( \begin{bmatrix} 1 & 1 \\ 1 & 1 \end{bmatrix} - \begin{bmatrix} 1 & -1 \\ -1 & 1 \end{bmatrix} \right)= \frac{1}{2} \begin{bmatrix} 0 & 2 \\ 2 & 0 \end{bmatrix} = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} $ This confirms that the spectral representation correctly reconstructs $A$. ## The Singular Value Decomposition The fundamental theorem of linear algebra is that every matrix is endowed with a spectral decomposition. Specifically, even if a matrix is not square, $A \in \mathbb{R}^{m \times n}$, then it admits a **Singular Value Decomposition** (SVD): $ A = U \Sigma V^T $ $U \in \mathbb{R}^{m \times m}$ is an orthogonal matrix whose columns are the **left singular vectors** of $A$, $V \in \mathbb{R}^{n \times n}$ is an orthogonal matrix whose columns are the **right singular vectors** of $A$, $\Sigma \in \mathbb{R}^{m \times n}$ is a diagonal matrix (possibly rectangular) containing the **singular values** of $A$, which are non-negative and ordered from largest to smallest. The columns of $U$ and $V$ provide orthonormal bases for the **range** and **domain** of $A$, respectively, and the singular values tell us how much $A$ stretches along each of those directions. This decomposition works for **any** real matrix—square or rectangular, symmetric or not—and forms the basis for many practical computations in scientific computing, such as low-rank approximation, data compression, and solving ill-posed systems. Let $A \in \mathbb{R}^{3 \times 2}$ be: $ A = \begin{bmatrix} 2 & 0 \\ 0 & 2 \\ -2 & 0 \end{bmatrix} $ Each row is a point in 2D: $(2, 0)$, $(0, 2)$, and $(-2, 0)$. We think of this as a data matrix of three 2D observations. First, compute $A^T A$: $ A^T A = \begin{bmatrix} 2 & 0 & -2 \\ 0 & 2 & 0 \end{bmatrix} \begin{bmatrix} 2 & 0 \\ 0 & 2 \\ -2 & 0 \end{bmatrix} = \begin{bmatrix} 8 & 0 \\ 0 & 4 \end{bmatrix} $ 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 is a reflection of the first, which is reflected in the structure of $U$. ## Spectral-Like Decomposition We can also express $A$ directly in a spectral-like form using its singular value decomposition. Since $A$ is not square, it does not have a traditional spectral decomposition, but the SVD provides a closely related structure: $ A = \sum_{i=1}^{r} \sigma_i \, \mathbf{u}_i \mathbf{v}_i^T $ where $\sigma_i$ are the singular values, $\mathbf{u}_i$ are the left singular vectors, and $\mathbf{v}_i$ are the right singular vectors. For the matrix $ A = \begin{bmatrix} 2 & 0 \\ 0 & 2 \\ -2 & 0 \end{bmatrix} $ we already found: * $\sigma_1 = 2\sqrt{2}, \quad \sigma_2 = 2$ * $\mathbf{u}_1 = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 0 \\ -1 \end{bmatrix}, \quad \mathbf{u}_2 = \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}$ * $\mathbf{v}_1 = \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad \mathbf{v}_2 = \begin{bmatrix} 0 \\ 1 \end{bmatrix}$ Then the decomposition becomes: $ 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} $ $ \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} $ So: $ A = \begin{bmatrix} 1 & 0 \\ 0 & 0 \\ -1 & 0 \end{bmatrix} \cdot 2 + \begin{bmatrix} 0 & 0 \\ 0 & 1 \\ 0 & 0 \end{bmatrix} \cdot 2 = \begin{bmatrix} 2 & 0 \\ 0 & 2 \\ -2 & 0 \end{bmatrix} $ This decomposition reveals $A$ as the sum of two rank-one matrices, each describing a directional stretching determined by the singular values and corresponding singular vectors. It’s the natural analog of spectral decomposition for arbitrary (non-square) matrices. ## Key Conclusion **Low-Rank Approximation**: Given $A\in\mathbb{R}^{m \times n}$, and whatever context its data represents, we always have access to the expression $A=\sum_{i=1}^{r} \sigma_i \, \mathbf{u}_i \mathbf{v}_i^T$ and if $\sigma_{1}\leq \sigma_{2}\leq\sigma_{3}\ll \sigma_{4}\leq \dots \leq \sigma_{r}$, then we could decide to only retain the first three terms of the sum to approximate $A$. Doing so would give us an approximation of the data in $A$ for the cost of six vectors and three scalars! ## Deliverable Checklist - [ ] Code - [ ] Computer codes that check the associated results in the mathematical discussion portion of this project - [ ] takes in image data and verifies - [ ] converts the image data to matrix data - [ ] acts on the matrix data with singular value decomposition libraries - [ ] Renders the data of the first few low-rank approximations - [ ] Communication - [ ] One-page overview of your solution - [ ] A five to seven minute recording of your walkthrough of this one-page overview and working code base