3
$\begingroup$

I am not very sophisticated in my linear algebra, so please excuse any messiness in my terminology, etc.

I am attempting to reduce the dimensionality of a dataset using Singular Value Decomposition, and I am having a little trouble figuring out how this should be done. I have found a lot of material about reducing the rank of a matrix, but not reducing its dimensions.

For instance, if I decompose using SVD: $A = USV^T$, I can reduce the rank of the matrix by eliminating singular values below a certain threshold and their corresponding vectors. However, doing this returns a matrix of the same dimensions of $A$, albeit of a lower rank.

What I actually want is to be able to express all of the rows of the matrix in terms of the top principal components (so an original 100x80 matrix becomes a 100x5 matrix, for example). This way, when I calculate distance measures between rows (cosine similarity, Euclidean distance), the distances will be in this reduced dimension space.

My initial take is to multiply the original data by the singular vectors: $AV_k$. Since $V$ represents the row space of $A$, I interpret this as projecting the original data into a subspace of the first $k$ singular vectors of the SVD, which I believe is what I want.

Am I off base here? Any suggestions on how to approach this problem differently?

$\endgroup$
2
  • $\begingroup$ Are you looking for a matrix $\tilde{A}$ of reduced dimension such that $\mathcal{R}(A^T) \approx \mathcal{R}(\tilde{A}^T)$ in some sense? $\endgroup$
    – copper.hat
    Jun 26, 2012 at 17:55
  • $\begingroup$ I'm not completely familiar with your notation, but I think that's right. I want to approximate the larger matrix with a matrix of reduced dimension while maintaining as much of the original variance as possible (hence the use of the SVD). $\endgroup$
    – Matt
    Jun 26, 2012 at 18:42

2 Answers 2

3
$\begingroup$

If you want to do the $100\times80$ to $100\times5$ conversion, you can just multiply $U$ with the reduced $S$ (after eliminating low singular values). What you will be left with is a $100\times80$ matrix, but the last $75$ columns are $0$ (provided your singular value threshold left you with only $5$ values). You can just eliminate the columns of $0$ and you will be left with $100\times5$ representation.

The above $100\times5$ matrix can be multiplied with the $5\times80$ matrix obtained by removing the last $75$ rows of $V$ transpose, this results in the approximation of $A$ that you are effectively using.

$\endgroup$
1
  • $\begingroup$ Thank you! This is exactly what I wanted to know. $\endgroup$
    – Matt
    Oct 1, 2012 at 14:37
1
$\begingroup$

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

$\endgroup$
1
  • $\begingroup$ Wow Lucy, thanks for all the great detail and the code sample! This makes a wonderful reference on this topic, and I still do a fair amount of SVD so I'll be coming back to this regularly. $\endgroup$
    – Matt
    Jan 20 at 13:38

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .