Singular Value Decomposition and Image Compression

Eigensystems as dominant “stretching” directions¶

The Spectral theorem states that if a matrix is Hermitian ($A^\dagger = A$, including real-symmetric), then: (1) The eigenvectors are an orthogonal basis, and (2) the eigenvalues are real. This is consistent with our intuition from the dynamical systems view: A “blob” of initial conditions gets anisotropically stretched at different rates along different directions. Eigenvalues are thus invariant to rotations (ie, multiplication by an orthonormal matrix)

Why are Hermitian matrices so important?¶

Many physical systems arise from the minimization of a scalar potential, such as a free energy or loss function. The dynamics of such systems are described by the gradient of the potential,

The Jacobian matrix of partial derivatives determines the local geometry and concavity of this scalar field, and thus stability and acceleration along various directions

Dynamical systems (continuous $\dot{\mathbf{x}} = \mathbf{f}(\mathbf{x})$ or discrete $\mathbf{x}_{t} = \mathbf{f}(\mathbf{x}_{t-1})$) always have a Jacobian, but they do not necessarily have a corresponding scalar potential, which requires that the Jacobian be Hermitian.

The loss function in machine learning. Often in machine learning, we need to compare to datapoints, or to compare a high-dimensional model prediction to some ground truth value. If we have two vectors $\hat{\mathbf{y}}, \mathbf{y}$, a loss function or distance function assigns a scalar value to the pair,

The gradient of the loss function with respect to vector $\hat{\mathbf{y}}$ describes how much the loss function changes as we make small changes to the prediction along each dimension.

For example, we often seek to minimize the mean squared error or Euclidean distance between a vector of predictions $\hat{\mathbf{y}}$ and a vector of ground truth values $\mathbf{y}$,

The gradient of the loss function with respect to $\hat{\mathbf{y}}$ is then the signed element-wise difference between the two vectors,

Now let’s repeat that iteration many times. We will apply our random Hermitian matrix $A$ many times to a cloud of points, and plot the result each time.

We can now plot the eigenvectors of $A$, which are the directions along which the cloud is stretched.

We can see that the eigenvectors of a Hermitian matrix indeed determine the dominant stretching directions of the cloud. What about the eigenvalues? To understand their role, we can project the point cloud onto the unit vectors $\mathbf{\hat{u}}_i$ associated with the eigenvectors of $A$. Given our point cloud $X_t$ at a given time $t$, $X_t \in \mathbb{R}^{N \times 2}$, we can project it onto the unit vector $\mathbf{\hat{u}}_i$ by taking the dot product $X_t \cdot \mathbf{\hat{u}}_i$:

where $i,j \in \{1, 2\}$ are the row and column indices of $X$. We can plot the maximum projection of the point cloud onto each eigenvector as a function of time, in order to measure the span of the point cloud along each eigenvector direction. We can compare this sequence to the geometric series of eigenvalues associated with each eigenvector, $\lambda_i^t$.

We can see that the eigenspectrum of Hermitian matrices has a straightforward interpretation, which is partially a consequence of the spectral theorem. The eigenvectors correspond to dominant stretching directions of a cloud of points under the action of a matrix, and the eigenvalues correspond to rates of compression or stretching per iteration. In the case of non-Hermitian matrices, we will observe rotation in addition to stretching, resulting in complex-valued eigenvalues and eigenvectors.

What about non-square matrices? Singular Value Decomposition¶

Suppose we have a matrix $A \in \mathbb R^{N \times M}$. For example, if we have a video of a fluid flow or sensor readings, we might have $N$ frames from $M$ sensors/pixels in each frame. If $A$ describes an algebraic system, we might have $N$ equations and $M$ unknowns, which can be overdetermined ($N > M$) or underdetermined ($N < M$).

However, we can instead calculate the eigenspectrum of the square matrices $A^\dagger A \in \mathbb R^{M \times M}$ or $A A^\dagger \in \mathbb R^{N \times N}$. We refer to the eigenvectors as the “right” or “left” singular vectors, respectively, and their associated eigenvalues are called the singular values. For an $N \times M$ matrix, we have $M$ right singular vectors and $N$ left singular vectors.

Note that, when $A$ is square, $A^\dagger A$ and $A A^\dagger$ share the same eigenvalues. To see this, let $\mathbf{x}$ be an eigenvector of $A^\dagger A$ with corresponding eigenvalue $a$. $A^\dagger A \mathbf{x} = a \mathbf{x} \implies A A^\dagger A \mathbf{x} = a \mathbf{x}$. Consequently, $a$ is also an eigenvalue of $A A^\dagger$ with corresponding eigenvector $A \mathbf{x}$.

For a non-square $A$, we can factorize $A$ as

where the columns of $U \in \mathbb R^{N \times N}$ are the right singular vectors, the columns of $V \in \mathbb R^{M \times M}$ are the left singular vectors, and $\Sigma \in \mathbb R^{N \times M}$ is a square diagonal matrix. The non-zero entries of $\Sigma$ are the square roots of the non-zero eigenvalues of $A^\dagger A$ or $A A^\dagger$. Since the eigenvalues and eigenvectors can be listed in any order, by convention we normally list the elements in descending order of eigenvalue magnitude, $|\sigma_i| \geq |\sigma_{i+1}|$.

The singular vectors decompose a dataset into characteristic directions in rowspace (left eigenvectors) and characteristic directions in columnspace (right eigenvectors).Importantly, unlike traditional eigenvalue decomposition, SVD applies to non-square matrices. Depending on the context, the singular vectors have different interpretations, including (1) Geometric, (2) Dynamical, (3) Data compression, pattern discovery

• See image below from Chris Rycroft’s AM205 class

The singular values $\sigma_i$ are the square roots of the eigenvalues of either $A^\dagger A$ or $A A^\dagger$. Because $A^\dagger A$ or $A A^\dagger$ are both symmetric matrices, the spectral theorem ensures that the singular values are always real and non-negative, and so are usually indexed so that they are listed in descending order, $\sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_k$.

The left singular vectors $U = \{\mathbf{u}_i\}$ are the eigenvectors of $A^\dagger A$. They are orthogonal unit vectors. The right singular vectors $V = \{\mathbf{v}_i\}$ are the eigenvectors of $A A^\dagger$. They are preimages of the $U$ vectors, so that $A \mathbf{v}_i = \sigma_i \mathbf{u}_i$. They are also orthogonal unit vectors.

Truncated versus full SVD¶

In full SVD, we calculate all $N$ singular values and vectors, and the matrices $U$ and $V$ are square with dimensions $N \times N$ and $M \times M$, respectively, while the matrix $\Sigma$ is rectangular with dimensions $N \times M$. In truncated SVD, we calculate only the first $k$ singular values and vectors associated with the largest singular values $\sigma_i$,

We define the truncated matrices $\tilde U$ of size $N \times k$, $\tilde \Sigma$ of size $k \times k$. $V$ still has size $M \times M$. The truncated SVD often approximates the original matrix $A$ well, making it useful for data compression and dimensionality reduction.

Image from from Chris Rycroft’s AM205 class

Dynamical systems interpretation of SVD¶

In dynamical systems, we might have skew-product structure, where the dynamics of some variables depend only on a subset of the others

Examples: $f(z_t) = z_t + c$ is a a time like variable, $f(z_t) = \eta$ is a random number generator (like we saw for the Ornstein-Uhlenbeck simulation); $f(z_t) = \sin(\omega t)$ is a periodic driver.

Can write this system has rectangular structure $$

=

Can't use function '$' in math mode at position 146: …he sequence of $̲z_t$ values. No…

Our analogy to square matrix dynamics is now a bit more involved, because we can't iterate the system repeatedly without knowing the sequence of $z_t$ values. Notice that if the coupling to the $z$ driver system is zero, then our square "response" subsystem represents an autonomous linear system with Hermitian structure: $x$ affects $y$ as much as $y$ affects $x$,

=

$$

When $z \ne 0$, we cannot do any spectral analysis anymore, because our coefficient matrix $A \in \mathbb{R}^{2 \times 3}$. However, we can instead analyze the two square matrices $AA^\dagger \in \mathbb{R}^{2 \times 2}$ and $A^\dagger A \in \mathbb{R}^{3 \times 3}$. This is like analyzing two iterations of a square dynamical system.

Loosely, the left matrix represents the propagation of an impulse in $z$ across the response subsystem (coupling of dynamical variables to the driver. This tells us about the dependence of response variables on earlier values

Loosely, the right matrix tells us about the driver system. If the driver couplings have the form $b = [3, 2]$, then the right matrix is perturbed by elementwise-adding the structured matrix $b b^T$ to the right matrix. The diagonal elements tell us something about how much the different dynamical variables are affected. This interpretation is reminiscent of the density matrix of mixed quantum systems.

Data representation and compression views¶

We consider the case where $A \in \mathbb R^{N \times M}$ is a data matrix, describing $N$ observations of an $M$ dimensional process. Examples include a video ($N$ is number of frames, $M$ is number of pixels), or a collection of $N$ points in $3D$ space ($M = 3$). The left eigenvectors give characteristic “modes” in the data; they have the same dimensionality. The right eigenvectors give us the characteristic “loadings” onto those modes. SVD can thusbe used to reduce the dimensionality of a problem by truncating relevant eigenvalues; a basic compression method for high-dimensional data. The decay of the singular values provides a clue into the intrinsic dimensionality of the underlying data generating process.

Application: Image compression¶

The SVD of an image can be used to compress the image by truncating the sum over singular values and vectors. Recall that the SVD of a matrix $A$ is given by

We can verify that this decomposition is valid by reconstructing the original matrix $A$ from its SVD components.

Numpy’s SVD solver returns the matrices $U, S, V^\dagger$. We use a keyword argument full_matrices=False to request reduced SVD, which returns $U$ with shape $N \times \min(N, M)$, $S$ with shape $\min(N, M)$, and $V^\dagger$ with shape $\min(N, M) \times M$. We can convert the vector $S$ to a diagonal matrix $\Sigma$ using the np.diag function.

Note that Numpy returns $V^\dagger$ directly, and so we do not need to transpose it when we compute the reconstruction.

SVD applies to any rectangular matrix. If $A$ corresponds to an $N \times M$ image, we can reconstruct individual images by treating them as collections of columns, and reconstructing by adding together subsets of characteristic columns,

We can also truncate this sum at some $k < \text{min}(N, M)$in order to approximate the matrix $A$ instead.

If the singular values decay rapidly, then we can approximate the matrix by truncating this sum to the first $k$ terms

We can now try varying the number of terms $k$ that we retain in the sum, and see how well we can reconstruct the original image.

We can directly look at the singular value distribution, to see if it matches our expectation based on the images.

Density matrices and product states¶

From our image compression example, we can see that we can approximate a matrix by a sum over rank-1 (outer product) matrices

where the sum is over all singular values and vectors.

For a Hermitian matrix, $A = A^\dagger$, we can instead write this as a sum over rank-1 matrices with the same vector in both slots

where the $\mathbf{u}_i$ are orthonormal eigenvectors and $\lambda_i$ are the corresponding eigenvalues.

Recall that, in quantum mechanics, that if we have a system with “pure” eigenstates $\Phi^\dagger \equiv \{\phi_i\}_N$, then we can define the density matrix for a given mixed state as

We can see truncated SVD is equivalent to approximating $\rho$ by only including the most-probable terms with largest $|s_i|$ (up to a renormalization), thus removing rare states.

The Moore-Penrose pseudoinverse¶

Suppose we have an overdetermined system of equations

where $A \in \mathbb{R}^{N \times M}$, $\mathbf{x} \in \mathbb{R}^M$, and $\mathbf{b} \in \mathbb{R}^N$. If the number of rows is greater than the number of columns ($N > M$), then we have more equations than unknowns, and the system is overdetermined. In this case, it is likely that no $\mathbf{x}$ exists that satisfies $A \mathbf{x} = \mathbf{b}$ exactly. However, we can instead seek the “best” solution in the sense of minimizing the squared error $||A \mathbf{x} - \mathbf{b}||^2$.

Deriving the pseudoinverse¶

We write the rectangular matrix $A$ in terms of its singular value decomposition

where, $U \in \mathbb{R}^{N \times N}$, $\Sigma \in \mathbb{R}^{N \times M}$, and $V \in \mathbb{R}^{M \times M}$. We multiply both sides by $U^\dagger \in \mathbb{R}^{N \times N}$ and note that $U^\dagger U = I$, producing the equation

We next define the pseudoinverse of $\Sigma$, denoted by $\Sigma^+ \in \mathbb{R}^{M\times N}$, as the matrix with the same dimensions as $\Sigma$ that has the reciprocal of each non-zero element along the primary diagonal, and zeros everywhere else. Note that this pseudoinverse has a shape that is transposed relative to $\Sigma$. This pseudoinverse has the property that $\Sigma^+ \Sigma = I \in \mathbb{R}^{M \times M}$. However, $\Sigma \Sigma^+ \neq I \in \mathbb{R}^{N \times N}$.

Our previous expression was

We multiply both sides of our expression above by $\Sigma^+$ to get

We can then multiply both sides by $V$ to get

We can then define the pseudoinverse $A^+$ of a rectangular matrix $A$ as

To understand the pseudoinverse, we consider the case where $N$ is greater than $M$, so that $A$ is overdetermined.

Within the $min(M,N)$ subspace associated with the $M$ columns of $A$, the pseudotinverse is self-adjoint, $A^+ = A^\dagger$ and so $A^+ A = I$.