Singular Value Decomposition and Image Compression

Eigensystems as dominant “stretching” directions¶

• Spectral theorem: 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

• 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 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¶

What am I leaving out?¶

• Non-Hermitian matrices rotate their input space, in addition to anisotopically squeezing it
• • A hot topic of research in stat mech and active matter theory; these describe non-reciprocal interactions (broken detailed balance if the matrix describes kinetics, “odd” elasticity, etc)

What about non-square matrices? Singular Value Decomposition¶

• $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

• $A$ 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.

  • Note: $A^\dagger A$ and $A A^\dagger$ share the same eigenvalues. 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}$$

    (6)
  • Consequently, $a$ is also an eigenvalue of $A A^\dagger$ with corresponding eigenvector $A \mathbf{x}$.

• We can now factorize $A$ as:

  $$A = U \Sigma V^\dagger$$

  (7)

  where the columns of $U$ are the right singular vectors, the columns of $V$ are the left singular vectors, and $\Sigma$ 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.

• Decompose a dataset into characteristic directions in rowspace (left eigenvectors) and characteristic directions in columnspace (right eigenvectors)

• Unlike traditional eigenvalue decomposition, SVD applies to non-square matrices. We can also apply SVD to complex matrices!

• Different interpretations: (1) Geometric, (2) Dynamical, (3) Data compression, pattern discovery

• See image below from Chris Rycroft’s AM205 class

From this diagram, we can define several quantities

• The singular values $\sigma_i$ are the square roots of the eigenvalues of either $A^\dagger A$ or $A A^\dagger$. They are always real and non-negative, and are usually indexed so that they are listed in descending order.

• The left singular vectors $U = \{\mathbf{u}_i\}$ are the eigenvectors of $A^\dagger A$. They are orthogonal unit vectors that point in the directions of the principal axes of $AS$

• 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$.

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.

• See image below for a visual representation of the truncated SVD 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 $$

=

$$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$/

We can’t 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 as propagation of an impulse in $z$ across the response subsystem (coupling of dynamical variables to the driver. 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 approach might remind you of density matrices

Data representation and compression views¶

• $A \in \mathbb R^{N \times M}$ is a data matrix, where we have $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”
• SVD can be used to reduce the dimensionality of a problem by truncating relevant eigenvalues; a basic compression method for high-dimensional data
• The eigenvalue spectrum provides a clue into the dimensionality of the underlying data generating process

Application: Image compression¶

Important: SVD applies to any rectangular matrix. We can reconstruct individual images by treating them as collections of columns, and reconstructing by adding together subsets of characteristic columns

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

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

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|$, 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 $N > M$, then we have more equations than unknowns, and the system is overdetermined.

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^+$, as the matrix with the same dimensions as $\Sigma$ that has the reciprocal of each non-zero element along the diagonal, and zeros everywhere else. If $\Sigma$ is square, then the pseudoinverse $\Sigma^+$ is the true inverse.

We can then multiply both sides 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

What does the pseudoinverse mean?¶

The pseudoinverse is a generalization of the inverse to non-square matrices. In cases where no $\mathbf{x}$ exists that satisfies $A \mathbf{x} = \mathbf{b}$, the pseudoinverse gives us the “best” solution in the sense that it minimizes the squared error $||A \mathbf{x} - \mathbf{b}||^2$.

Appendix¶