Show code cell source

Hide code cell source

import matplotlib.pyplot as plt
%matplotlib inline
import matplotlib_inline
matplotlib_inline.backend_inline.set_matplotlib_formats('svg')
import seaborn as sns

Singular Value Decomposition

Singular value decomposition (SVD) is a factorization of a matrix into three matrices. It is used in many applications such as data compression, denoising, and solving linear systems of equations. In scientific machine learning, it is also used in principal component analysis (PCA) and the numerical construction of Karhunen–Loève expansions.

A broader treatment appears in Brunton and Kutz (2019).

Let \(\mathbf{X}\) be an \(n \times m\) matrix. Think of \(\mathbf{X}\) as the matrix you obtain when doing \(n\) experiments and measuring \(m\) different things. The SVD of \(\mathbf{X}\) is given by

\[ \mathbf{X} = \mathbf{U} \mathbf{\Sigma} \mathbf{V}^T \]

where \(\mathbf{U}\) is an \(n \times n\) orthogonal matrix, \(\mathbf{\Sigma}\) is an \(n \times m\) matrix with non-negative real numbers on the diagonal and zeros elsewhere, and \(\mathbf{V}\) is an \(m \times m\) orthogonal matrix.

The columns of \(\mathbf{U}\) are called the left singular vectors of \(\mathbf{X}\), the columns of \(\mathbf{V}\) are called the right singular vectors of \(\mathbf{X}\), and the diagonal elements of \(\mathbf{\Sigma}\) are called the singular values of \(\mathbf{X}\).

Economy-size SVD

Assume that \(n \geq m\). Then, \(\mathbf{\Sigma}\) has the form:

\[\begin{split} \mathbf{\Sigma} = \begin{bmatrix} \hat{\mathbf{\Sigma}} \\ \mathbf{0} \end{bmatrix}, \end{split}\]

where \(\hat{\mathbf{\Sigma}}\) is an \(m \times m\) diagonal matrix with non-negative real numbers on the diagonal and zeros elsewhere, and \(\mathbf{0}\) is an \((n-m) \times m\) matrix with zeros. Now, only the first \(m\) columns of \(\mathbf{U}\) are needed to represent \(\mathbf{X}\). We write the economy-size SVD as:

\[ \mathbf{X} = \mathbf{U}_m \hat{\mathbf{\Sigma}} \mathbf{V}^T, \]

where \(\mathbf{U}_m\) is an \(n \times m\) matrix with the first \(m\) columns of \(\mathbf{U}\).

Truncated SVD

The truncated SVD is a low-rank approximation of \(\mathbf{X}\). It is given by:

\[ \mathbf{X} \approx \mathbf{U}_k \hat{\mathbf{\Sigma}}_k \mathbf{V}_k^T, \]

where \(\mathbf{U}_k\) is an \(n \times k\) matrix with the first \(k\) columns of \(\mathbf{U}\), \(\hat{\mathbf{\Sigma}}_k\) is a \(k \times k\) diagonal matrix with the first \(k\) singular values of \(\mathbf{\Sigma}\), and \(\mathbf{V}_k\) is an \(m \times k\) matrix with the first \(k\) columns of \(\mathbf{V}\).

We can also write:

\[ \mathbf{X} \approx \mathbf{X}_k = \sum_{i=1}^k \sigma_i \mathbf{u}_i \mathbf{v}_i^T, \]

where \(\sigma_i\) is the \(i\)-th singular value, and \(\mathbf{u}_i\) and \(\mathbf{v}_i\) are the \(i\)-th columns of \(\mathbf{U}\) and \(\mathbf{V}\), respectively.

One can show that the matrix \(\mathbf{X}_k\) is the best rank-\(k\) approximation of \(\mathbf{X}\) in the Frobenius norm. The Frobenius norm is the square root of the sum of the squared matrix entries, \(\|\mathbf{A}\|_F=\sqrt{\sum_{i,j}A_{ij}^2}\), so it measures the total entrywise approximation error.

Demonstration - SVD for image compression

We compress the grayscale image below using the truncated SVD.

from PIL import Image
with open("neom-DMGDdksVoWI-unsplash.jpg", "rb") as f:
    img = Image.open(f)
    img = img.convert("L")  # convert to grayscale
img

Extract the matrix that represents the image:

import numpy as np
X = np.array(img)
X.shape

(1798, 2400)

Do SVD of the matrix:

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

Let’s look at the singular values as a function of the index:

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(s)
ax.set(yscale='log', xlabel="Singular value index", ylabel="Singular value")
finalize_axes(keep_box=False)

array([<Axes: xlabel='Singular value index', ylabel='Singular value'>],
      dtype=object)

Typically, we pick the number of singular values such that the sum of the squares of the singular values is a certain percentage of the sum of the squares of all the singular values. For example, we can pick \(k\) so that:

\[ \frac{\sum_{i=1}^k \sigma_i^2}{\sum_{i=1}^{\min(n,m)} \sigma_i^2} \geq 0.998. \]

To facilitate our choice we can do the following plot:

Show code cell source

Hide code cell source

cum_var = np.cumsum(s**2)/np.sum(s**2)
best_k = np.argmax(cum_var > 0.998)

fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_standard"])
ax.plot(cum_var)
ax.axhline(0.998, color='red', linestyle='--')
ax.axvline(best_k, color='red', linestyle='--')
ax.text(best_k + 10, 0.99, f'99.8% at k={best_k}', verticalalignment='center', color='red')
ax.set(yscale='log', xlabel="Singular value index", ylabel="Cumulative fraction of total variance")
finalize_axes(keep_box=False)

array([<Axes: xlabel='Singular value index', ylabel='Cumulative fraction of total variance'>],
      dtype=object)

Let’s plot the compressed images for some choices of \(k\):

ks = [1, 2, 4, 8, 16, 32, 64, 128, 256, 682]

for k in ks:
    X_k = U[:, :k] @ np.diag(s[:k]) @ Vt[:k, :]
    fig, ax = plt.subplots(figsize=FIGURE_SIZES["half_tall"])
    ax.imshow(X_k, cmap='gray')
    ax.axis('off')