# How to Calculate Singular Value Decomposition (SVD) in Python

Singular value decomposition (SVD) is one of the most significant matrix decomposition techniques. It should be in the toolbox of anyone working on linear algebra.

Python has many great packages that provide support for basic linear algebra. We will go through some of these, and also see how to create a naive python-only SVD implementation without any external packages:

After reading this post, you can perform SVD using the following packages

1. NumPy: This is our go-to package for any numerical computations using Python. It provides great Python support for linear algebra, particularly to perform SVD, we can use  function numpy.linalg.svd.
2. SciPy: Whenever something is not directly provided by NumPy, the next best bet is to check SciPy. It provides a wide range of functions and algorithms for scientific computing in Python. It has also it's one linear algebra module from which we can find a function scipy.linalg.svd that computes the SVD.
3. scikit-learn: This is one of the more popular machine-learning packages for Python. Machine learning is a linear algebra-heavy and computationally intensive field, so we can also find operations for some popular linear algebra operations within the machine learning packages. scikit-learn has a class TruncatedSVD, which can use to perform a special type of SVD, where only the most significant singular values are considered and the corresponding matrices can be provided in a condensed form.
4. PyTorch: PyTorch is another machine learning package for Python, which is more steered towards artificial intelligence and deep learning. Like scikit-learn, PyTorch also provides interfaces to some linear algebra operations. For SVD, we can utilize a function torch.svd.

In addition, we provide a pure Python implementation using the power iteration method for calculating dominant singular values and vectors.

## A short introduction to Singular Value Decomposition

Singular value decomposition or SVD, in short, is a matrix decomposition method, which can be used to decompose a matrix into a product of three matrices with very specific structures.

Consider for instance a matrix $$\mathbf{A}$$ with dimension $$n \times m$$. Using SVD, we can decompose this matrix into $\mathbf{A} = \mathbf{U}\mathbf{D}\mathbf{V}^\text{H},$
where

• $$\mathbf{U}$$ is the left singular matrix, which has orthonormal column vectors and dimensions $$n \times n$$. The column vectors are called left singular vectors.
• $$\mathbf{D}$$ is the diagonal singular value matrix with dimensions $$n \times m$$. The diagonal values in $$\mathbf{D}$$ are the singular values of $$\mathbf{A}$$. These are non-negative numbers that are usually presented in decreasing order.
• $$\mathbf{V}$$ is the right singular matrix, which has orthonormal column vectors and dimensions $$m \times m$$. The column vectors are called right singular vectors.

Here, orthonormal means that the corresponding column vectors have unit-norm and are orthogonal to each other. This essentially means that the following holds $\mathbf{U}\mathbf{U}^\text{H} = \mathbf{I},$ where $$\mathbf{I}$$ is identity matrix.

### Applications

It's quite safe to say that SVD is one of the most important matrix factorization techniques. Thus, there are numerous applications for SVD. It is used in

• Efficient matrix inversion: singular value decomposition can greatly simplify matrix inversion for some specific matrix structures. It is beneficial in algorithms, where repetitive inversion is required with only minor changes in a matrix structure.
• Data compression: singular value decomposition can reveal the most significant linear components (directions, subspaces, etc...) that hold the most information. This has obvious applications in data compression. See for instance our recent post on using PCA in image compression. PCA uses SVD internally to calculate components holding the most variance.
• Data denoising: similar to data compression, noise removal singular value decomposition can be used to remove noise in the data.
• Any other linear algebra-heavy application: SVD is applied in many linear algebra-heavy processes such as AI, machine learning, latent semantic indexing, collaborative filtering, and natural language processing.

## Calculating SVD in Python

As mentioned before there are various packages that provide an interface to the SVD. We will shortly go through a few of the more popular ones.

First, we need to create a test data set to have a common point of reference in the comparison of the results.

We define our input matrix as $A = \begin{bmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{bmatrix}$

### NumPy

NumPy provides a simple interface for the SVD via np.linalg.svd. All we need to do is pass the input matrix, and we get the left/right singular matrix and singular values as an output. Simple!

Let's then verify that everything looks as assumed.

The result seems to be an identity matrix as expected. Thus, all good.

Next, we can do a similar check for the right singular vectors.

All good with the right singular vectors also.

The last item to check is that we can reproduce the singular values by multiplying the input matrix from the left and right with the corresponding singular matrices.

This gives us a matrix with the assumed singular values on the diagonal.

These match the singular values, which we received from the SVD: print(D).

### SciPy

SciPy is our favorite scientific computation package. The go-to place to check for functionality whenever it cannot be found from NumPy. SciPy provides a very similar interface for SVD as NumPy. It can be accessed from scipy.linalg.svd.

Let's compare the results to the NumPy.

print(D)

[16.848  1.068  0.   ]

Pretty much as expected, the results are identical to the ones provided by NumPy.

### scikit-learn

Moving to the machine learning libraries. Scikit-learn provides an interface to truncated SVD via sklearn.decomposition.TruncatedSVD, which can be used to calculate a subset of singular values for a given matrix (n_components). Unfortunately, we can only receive the right singular vectors. However, if you are only interested in a few dominant singular values, this may be a suitable solution

This should print out the singular values as

### PyTorch

Another popular machined learning package we are considering is PyTorch. It provides access to the singular value decomposition via torch.linalg.svd, which behaves similarly to NumPy and SciPy with the exception that the input type needs to be a tensor

The resulting singular values are

### Python-only singular value computation

Sometimes you can not handle the overhead of NumPy or SciPy and you need to work with vanilla Python. Fortunately, there are algorithms for SVD that are fairly straightforward to implement and provide reasonable performance for modestly sized matrices. To this end, we will apply the Power iteration algorithm for SVD.

In short, power iteration is a method that finds the largest singular vector for a given matrix. The process is such that, we find dominant vectors and values one-by-one and iterative remove them in the process.

#### Power iteration method basics

The power iteration method finds the dominant eigenvalue and vector for matrix $$\mathbf{A}\mathbf{A}^\text{T}$$. Noting from the SVD properties, we have $$\mathbf{A} = \mathbf{U}\mathbf{D}\mathbf{V}^\text{T}$$. Thus, we get $\mathbf{A}\mathbf{A}^\text{T} = \mathbf{U}\mathbf{D}\mathbf{V}^\text{T}\mathbf{V}\mathbf{D}^\text{T}\mathbf{U}^\text{T} = \mathbf{U}\mathbf{D}^2\mathbf{U}^\text{T},$ where the last steps comes from $$\mathbf{V}^\text{T}\mathbf{V} = \mathbf{I}$$ and $$\mathbf{D}$$ being diagonal matrix with singular values on the diagonal.

Following this, what the power iteration method provides us is iteratively the matrix $$\mathbf{U}$$ and $$\mathbf{D}^2$$.

1. Let $$n = 1$$, inialize search space to $$\mathbf{B}_{base}^{n} = \mathbf{A}\mathbf{A}^\text{T}$$, and randomly initialize vector $$\mathbf{b}^n_1$$.
2. Iterative until convergence or a fixed number of iterations $$k = 1,\ldots$$
$\mathbf{b}_{k + 1} = {\mathbf{B}_{base}^{n}\mathbf{b}_{k} \over \|\mathbf{B}_{base}^{n}\mathbf{b}_{k} \|}.$
3. Remove the last singular vector direction from the search base: $\mathbf{B}_{base}^{n +1} = \mathbf{B}_{base}^{n} - d^n \mathbf{b}^n_k (\mathbf{b}^n_k)^\text{T}.$
4. Get the singular values from the Rayleigh quotient
$d_n = \sqrt{{(\mathbf{b}^n_k)^\text{T} \mathbf{B}_{base}^n (\mathbf{b}^n_k) \over \|(\mathbf{b}^n_k)\|}}$
5. The $$n^\text{th}$$ singular vector is $$\mathbf{b}_k^n.$$

Now, we have $$\mathbf{U}$$ and $$\mathbf{D}$$. Finally, we can obtain, the right singular matrix by $\mathbf{V}^\text{T} = \mathbf{D}^{-1}\mathbf{U}^T \mathbf{A}$.

That is it. This is simple enough to implement in vanilla Python. You can find the implementation steps below.

#### Preliminaries

Before getting into the actual power iteration algorithm, we need to introduce a few simple utility functions. We will focus on code simplicity rather than optimal performance.

One thing that we need is matrix multiplication by a vector and vector multiplication by another vector. We call this dot, short for dot-product

Another one, which we will need is a matrix transpose. This will simply flip the order of rows & columns for a given input array

### Power iteration algorithm

Now, we are ready to move on to implementing the power iteration algorithm. Our implementation pretty much follows the steps listed above and should be quite self-explanatory.

This should print out similar results to the ones that, we got from the external packages

This works, but there is a lot of room for optimization. However, for small matrices, this works fine. Also, extending this to complex numbers is quite straightforward.

## Summary

We went through 5 ways of computing singular value decomposition in Python. As usual, you can use the preexisting packages like NumPy or SciPy. They utilize optimized algorithms for computing SVDs, which will really pay off when the size of the matrices increases. Although, if you are in a situation, where external packages are not feasible, you can use the power iteration method as a starting point for pure Python implementation.