# The Singular Value Decomposition

## 1. Using a matrix to store our data
We imagine that the data we are interested in are shaped into a $n\times m$ matrix, $\mathbf{X}$. Most of our data will be real numbers only, so $\mathbf{X}\in\mathbb{R}^{n\times m}$.

The data in $\mathbf{X}$ might be spatially 2-D (a single image, for example) or it might be composed such that each column is data from a set of spatial postions at different instants in time,

\begin{equation}
\mathbf{X} = \begin{bmatrix} 
\mid & \mid &  & \mid  \\ 
\mathbf{x_1} & \mathbf{x_2} & \cdots  & \mathbf{x_m} \\
\mid & \mid &  & \mid \end{bmatrix}\quad\text{,}
\end{equation}

where each column $\mathbf{x_k}$ is data from a set of $n$ sensors at a specific instant in time.

## 2. Definition of the Singular Value Decomposition

Any matrix $\mathbf{X}$ can be written as the product of 3 special matrices, $\mathbf{U}$, $\mathbf{\Sigma}$ and $\mathbf{V^T}$ (where $\mathbf{V^T}$ is the transpose of $\mathbf{V}$),

\begin{equation}
\mathbf{X} = \mathbf{U}\,\mathbf{\Sigma}\,\mathbf{V^T}\quad\text{.}
\end{equation}

$\mathbf{U}$ and $\mathbf{V}$ are both *orthogonal* matrices because,

\begin{align}
\mathbf{U}\mathbf{U^T} &=\mathbf{U^T}\mathbf{U}=\mathbf{I}\quad\text{and} \\
\mathbf{V}\mathbf{V^T} &=\mathbf{V^T}\mathbf{V}=\mathbf{I}\quad\text{.}
\end{align}

The $\mathbf{\Sigma}$ matrix only has non-zero entries on the diagonal. These are the **singular values** and are arranged in descending order.

(If our data contains complex numbers, $\mathbf{X}\in\mathbb{C}^{n\times m}$, then the SVD still works but now $\mathbf{U}$ and $\mathbf{V}$ are *unitary* matrices so that $\mathbf{U}\mathbf{U^*}=\mathbf{I}$ where $^*$ denotes conjugate transponse.)

## 3. Shapes of the matrices in the SVD

If $\mathbf{X}\in\mathbb{R}^{n\times m}$ then $\mathbf{U}\in\mathbb{R}^{n\times n}$, $\mathbf{\Sigma}\in\mathbb{R}^{n\times m}$ and $\mathbf{V^T}\in\mathbb{R}^{m\times m}$.



Let's see how we can evaluate the SVD using Python:


```python
import numpy as np 
import matplotlib.pyplot as plt
%matplotlib inline
```


```python
m = 3
n = 5
X = np.random.rand(n,m)     # Create a matrix of random numbers (between 0. and 1.) with n rows and m columns
```


```python
X                           # In a notebook, this writes X to the screen - similar to print(X)
```


```python
U,S,VT = np.linalg.svd(X)   # Use numpy's linalg.svd to evaluate the SVD
```


```python
U.shape, S.shape, VT.shape  # Find the shapes of U, S and VT 
```

We expected $\mathbf{U}\in\mathbb{R}^{n\times n}$, $\mathbf{\Sigma}\in\mathbb{R}^{n\times m}$ and $\mathbf{V^T}\in\mathbb{R}^{m\times m}$, but `numpy.linalg.svd` has returned $\mathbf{\Sigma}$ as a vector of length $n$. `numpy` has returned a compact description of $\mathbf{\Sigma}$, noting that it is a diagonal matrix and, since $n>m$ in our case, the bottom $n-m$ rows of $\mathbf{\Sigma}$ are zeros.

If we need to, we can reconstruct the full $\mathbf{\Sigma}$ using:


```python
Sfull = np.zeros([n,m])     # Make a matrix (array) that is all zeros and of the correct shape
Sfull[:m,:m] = np.diag(S)   # Insert the diagonal matrix made of the vector S. 
                            # The [:m,:m] notation means the first m rows and first m columns
```


```python
Sfull                       # display the full Sigma matrix (note the order of the singular values)
```

Now we can check the SVD by reconstructing X:


```python
Xrecon = U @ Sfull @ VT     # @ is matrix multiplication in Python 3
Xrecon 
```

## 4. Properties of U and V

We can also demonstrate some properties of $\mathbf{U}$ (and $\mathbf{V}$). The vectors that form the columns of $\mathbf{U}$ (and the columns of $\mathbf{V}$) have unit magnitude and are orthogonal to the other columns of the matrix, i.e. they are orthonomal.


```python
U[:,0] @ U[:,0].T           # Check that the magnitude of the first column of U = 1.
```


```python
U[:,0] @ U[:,1].T           # Check that the first column is orthogonal to the second column
```

---

Key intepretation of $\mathbf{U}$ and $\mathbf{V}$ (from Brunton and Kutz):

* Columns of $\mathbf{U}$ are an orthonormal basis for the column space of $\mathbf{X}$ 
* Columns of $\mathbf{V}$ (not $\mathbf{V^T}$ are an orthonormal basis for the row space of $\mathbf{X}$ 

If each column of $\mathbf{X}$ contains spatial data at a different time instant, then $\mathbf{U}$ relates to spatial patterns, and $\mathbf{V}$ to tempoeral patterns.

---


## 5. Matrix appoximation using the SVD

The SVD identifies the *rank* of a matrix by the number of singular values. The rank is the maximum number of linearly independent columns (or rows) of the matrix. 


```python
m = 3
n = 5
X = np.random.rand(n,m)      # Make a new (n,m) matrix of random numbers
U,S,VT = np.linalg.svd(X)    # Compute the SVD
S                            # Show the singular values (since X is made of random numbers, we expect all columns (and rows)
                             # to be independent and so X will be of rank m
```


```python
X[:,2] = X[:,1]*2.           # Change the third column (indices start from 0 in Python) to be twice the second column
U,S,VT = np.linalg.svd(X)    # Now when we compute the SVD, we expect two singular values because the third column is not idependent of the second
S
```

We can approximate $\mathbf{X}$ in a lower rank using the SVD. The optimal rank $r$ approximation is given by the largest $r$ singular values. A nice way to show this is by using an image as the data matrix $\mathbf{X}$.


```python
img = np.load('data/img.npy')
plt.imshow(img, cmap = 'gray', vmin = 0, vmax = 1)
```


```python
U,S,VT = np.linalg.svd(img)     # Perform the SVD
```


```python
S.size                          # Tells us the rank of the img matrix
```

We now approximate the `img` matrix using a new matrix of rank $r$.



The SVD now contains only the leading $r$ terms of $\mathbf{\Sigma}$ (the largest $r$ singular values), $\mathbf{\hat{U}}$ is the first $r$ columns of $\mathbf{U}$ and $\mathbf{\hat{V}^T}$ is the first $r$ rows of $\mathbf{V^T}$.


```python
r = 50                                                 # Rank of new matrix
imgRecon = U[:,:r] @ np.diag(S[:r]) @ VT[:r,:]         # Use first r columns of U, first r values in S, and first r rows in VT
plt.imshow(imgRecon, cmap='gray', vmin = 0, vmax = 1)
plt.title("Approximation; Rank = " + str(r) ) 
```

This represents a signficant reduction in the amount of data needed to describe the image:


```python
size_orig = img.flatten().size
size_approx = U[:,:r].flatten().size + S[:r].flatten().size + VT[:r,:].flatten().size
size_approx / size_orig
```

Even a rank 10 approximation is surprisingly clear: 


```python
r = 10 
imgRecon = U[:,:r] @ np.diag(S[:r]) @ VT[:r,:]
plt.imshow(imgRecon, cmap='gray')
plt.title("Approximation; Rank = " + str(r) ) 
```

The matrix formed when a single column of $\mathbf{U}$ and a single row of $\mathbf{V^T}$ are multiplied together is of rank 1. Our appoximation is the superposition of $r$ of these rank 1 matrices, each multiplied by a their respective singular value. 

We can inspect the individual rank 1 matrices:


```python
f=plt.figure(figsize=(10,5))
for k in range(10):
    plt.subplot(2,5,k+1)
    plt.title('k = ' + str(k))
    rank1matrix = np.outer(U[:,k] , VT[k,:])   # Form the rank 1 matrix
    plt.imshow(rank1matrix, cmap='gray')
    plt.axis('off')
```

The approximation of a matrix $\mathbf{X}$ up to rank $r$ is often given as an example of data compression. In engineering, we take a slightly different perspective. If a reduced rank matrix is a good approximation to the original data, then this tells us about the underlying low rank structure of our data. It is often this low rank structure engineers can apply from one example of the problem to the next; it is both interpretable and transferrable.

The SVD is useful in a wide range of data science applications. In the next notebook, we will look at one of the most common, Principal Component Analysis.


```python

```

134c79b5b71f1318e86670c86c9da6fbc3a81365
ipynb
Jupyter Notebook
01 - Singular Value Decomposition.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
01 - Singular Value Decomposition.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
01 - Singular Value Decomposition.ipynb
grahampullan/datascienceBB
87c1127d9deddd3e4e11f1a785b4f9e93a74fb19
MIT
Qwen/Qwen-72B
1. YES
2. YES