SVD can be used for the best matrix approximation of matrix $X$ in terms of the Frobenius norm. For instance, images can be compressed, or in the recommender system, we can approximate recommendations through lower-rank matrices. Suppose we define a new matrix $U_k \in \mathbb{R}^{s \times k}$ as the first $k$ columns of $U$. Then we observe that our k-rank approximated matrix $X_k$
$$X_k=X U_k U_k^T=U_k U_k \text{U$\Sigma $V}^T=U_k \left(I_{k\times k} \,\,\,0_{k\times (s-k)}\right)\text{$\Sigma $V}^T =\text{U$\Sigma $}_k V^T=X_k$$
where $\Sigma _k\in \mathbb{R}^{s\times n}$ as
$$\Sigma _k=\left(
\begin{array}{ccccccc}
\sigma _1 & 0 & \ldots & 0 & 0 & \cdots & 0 \\
0 & \sigma _2 & \ldots & 0 & 0 & \cdots & 0 \\
0 & 0 & \ddots & 0 & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & \sigma _k & 0 & \cdots & 0 \\
0 & 0 & \ldots & 0 & 0 & \cdots & 0 \\
\vdots & \vdots & \ddots & 0 & \vdots & \ddots & \vdots \\
0 & 0 & \ldots & 0 & 0 & \ldots & 0 \\
\end{array}
\right)$$
The above fact is expressed as Eckart-Young-Mirsky theorem that states that for any data matrix $X$ and any matrix $\hat{X}$ with $rank(\hat{X})=k$, we have the error of approximation as follows
$$\left\| X-\hat{X}\right\| _{\text{Fro}}^2\geq \left\| X-X_k\right\| _{\text{Fro}}^2=\left\| X-X U_k U_k^T\right\| _{\text{Fro}}^2=\sum _{j\geq k+1}^{\min (s,n)} \sigma _j^2$$
Thus by computing SVD of the matrix $X$ and taking k-first columns of $U$ we obtain the best k-rank matrix approximation in terms of the Frobenius norm. The last term in the above formula is the sum of singular values beyond $X_k$. So alternatively, we can use the k-singular values of diagonal matrix $\Sigma$, i.e. to truncate all the eigenvalues after k.
Let's compute approximations of the original image by setting all singular values after the given index to zero and print the ratios of the Frobenius norm of approximations to the original matrix. As you will see from the plotted pictures, the closer the ratio is to one, the better an approximation. We will work on the image as below.
import numpy as np
#!pip install scikit-image
from skimage import data
image1 = data.astronaut()
print(image1.shape)
plt.imshow(image1)
plt.show()
matrix = image1.reshape(-1, 1536) # -> (512, 1536)
The first step is to perform SVD decomposition.
U, sing_vals, V_transpose = np.linalg.svd(matrix)
To approximate our image from SVD matrices, we form a diagonal matrix $\Sigma_k$ based on the thresholded singular values, i.e., we zero out all the singular values after the threshold. To recover the approximation matrix, we multiply matrices $X_{k}= U \Sigma_k V^T$.
The approximation results are shown below.
#To perform PCA center the data
matrix = matrix - np.mean(matrix, axis=1).reshape(matrix.shape[0], 1)
#In practice the full standardisation is often applied
#matrix = matrix/np.std(matrix, axis=1).reshape(matrix.shape[0], 1)
def plot_approx_img(indexes, U, sing_vals, V_transpose ):
_, axs = plt.subplots(nrows=len(indexes)//2, ncols=2, figsize=(14,26))
axs = axs.flatten()
for index, ax in zip(indexes, axs):
sing_vals_thresholded = np.copy(sing_vals)
sing_vals_thresholded[index:] = 0 # zero out all the singular values after threshold
ratio = np.sum(sing_vals_thresholded**2)/np.sum(sing_vals**2) # compute the ratio
Sigma = np.zeros(matrix.shape)
#form zero-truncated matrix of singular values
Sigma[:matrix.shape[0], :matrix.shape[0]] = np.diag(sing_vals_thresholded)
matrix_approx = U@(Sigma@V_transpose) # approximated matrix
#image reshaping and clipping
image_approx = matrix_approx.reshape(image1.shape)
img = np.rint(image_approx).astype(int).clip(0, 255)
ax.imshow(img)
ax.set_title(['Threshold is '+str(index), 'Explained in ' + str(round(ratio, 4))])
thresholds = np.array([2] + [4*2**k for k in range(7)])
plot_approx_img(thresholds, U, sing_vals, V_transpose)
thresholds = np.array([2] + [4*2**k for k in range(7)])
plot_approx_img(thresholds, U, sing_vals, V_transpose)
So as many singular values, you use so much our lower rank matrix (image) is approximated. Look at the picture below
Look at approximations from SVD here
For more details, look at https://lucynowacki.github.io/blog/svd/index.html
