Additional links:
- Original class material
<< Previous lecture | Next lecture >>
todo clean
Singular Value Decomposition of a Matrix
In Data Analysis, it is often required to compress the data, either to make it more manageable or to be able to visualize the most important features, reducing redundancy. The two tasks, usually named data compression and dimensionality reduction, are mathematically equivalent and closely related to the concept of Singular Value Decomposition (SVD) of a matrix.
Recall that:
An invertible matrix is said to be orthogonal if or, equivalently, .
Now, consider a given matrix . It is known that it can be factorized into the product of three matrices,
where and are orthogonal matrices, while is a rectangular matrix which is non-zero only on the diagonal. Such decomposition is named Singular Value Decomposition (SVD) of .
Of particular interest in our analysis are the values on the diagonal of , named singular values of , and usually denoted as . In particular, it is known that the singular values:
- are always greater or equal to 0, i.e. , ;
- are ordered in descending order, i.e. $\sigma_1 \geq \sigma_2 \geq \dots \geq 0;
- can be used to determine the rank of , since it is equal to the number of singular values strictly greater than zero, i.e. if and for some index , then .
A useful properties of the SVD of is that it can be used to compress the informations contained in itself. Indeed, note that the SVD decomposition allows to rewrite as the sum of simple matrices, i.e.
where and are the columns of and , respectively. Each term is a rank-1 matrix named dyad, and the -th singular value represent the importance of the -th dyad in the construction of . In particular, the SVD decomposition allows to deconstruct into the sum of matrices with decreasing information content.
The SVD decomposition can be used to compress the matrix by considering its -rank approximation , defined as
It has been already showed that the -rank approximation of is the -rank matrix that minimizes the distance (expressed in Frobenius norm) from , i.e.
Implementation: computing the SVD of a Matrix
The functions required to compute SVD decomposition of a matrix in Python are contained into the numpy
package. In the following, we will consider the example matrix reported into the code snippet below.
Then, we consider the shape of and we save that value into the variables and , to follow the mathematical notation above:
Now we are ready to compute the SVD of . This can be done by running the function np.linalg.svd()
, that takes as input a matrix and returns a triplet U, s, VT
, representing the matrices , and a vectorized version of that only contains the diagonal (to save memory!).
Note that, to do some computation, it can be useful to compute the full matrix .
A first, sanity check is to verify that the algorithm works as expected. To do that, we should verify that (left as an exercise to the students)
Application 1: Compute the numerical rank of
As a first exercise, you are asked to compute the rank of by using the equation
Note that the result has to be compared with the output of the built-in function in numpy
:
Visualize the solution
<pre>
Check the rank of A by s
r = np.sum(s > 1e-10) print(fβThe computed rank of A is: {r}.β)
Application 2: The -rank approximation
You are now asked to compute the -rank approximation of . In particular, you are asked to:
- Given an integer , compute the -rank approximation .
- Compute the approximation error in Frobenius norm, .
Visualize the solution
<pre>
Given k, compute the k-rank approximation A_k of A
k = 3 A_k = U[:, :k] @ S[:k, :k] @ VT[:k, :] print(fβShape of A_k: {A_k.shape}. Rank of A_k: {np.linalg.matrix_rank(A_k)}.β)
Compute the error in Frobenius norm due to the approximation A_k
print(fβ||A - A_k||_2 = {np.linalg.norm(A - A_k, 2)}.β)
Application 3: SVD for Image Compression
From a computational point of view, a grey-scale image is a matrix with shape , such that the element in position contains the intensity of the pixel in the corresponding position. An RGB image is a triplet of matrices such that in position , each of the three matrices represents the amount of Red, Green and Blue in the corresponding pixel.
For example, consider an image from the skimage.data
submodule:
To visualize a matrix as an image, it can be used the plt.imshow()
function from the matplotlib.pyplot
module. If the image is a grey-scale image, it is required to set the cmap='gray'
option to visualize it as a grey-scale image.
Besides the visualization, the image is still a matrix and all the techniques developed for matrices can be used on images. In particular, the -rank approximation, by truncating the least important informations of the image, can be used to compress it with minimal information loss.
Note that, the -rank approximation of the image can be written as
where each is a scalar number, is an -dimensional vector, while is an -dimensional vector. As a consequence, the number of values required to memorize is , while the number of values required to memorize the whole image is . As a consequence, the compression factor of can be computed as
As an exercise, compute the -rank approximation of the cameraman image for different values of and observe the behavior.