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}\]

# Input test matrix for SVD
input_matrix = [
    [1.0, 2.0, 3.0],
    [4.0, 5.0, 6.0],
    [7.0, 8.0, 9.0]
]

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!

import numpy as np

U, D, V = np.linalg.svd(input_matrix)

Let's then verify that everything looks as assumed.

# Check that left singular matrix is orthonornal (unitary)
np.set_printoptions(precision=3, suppress=True)
print(U @ U.transpose())

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

array([[ 1.,  0., -0.],
       [ 0.,  1., -0.],
       [-0., -0.,  1.]])

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

# Check that left singular matrix is orthonornal (unitary)
np.set_printoptions(precision=3, suppress=True)
print(V @ V.transpose())

All good with the right singular vectors also.

array([[ 1.,  0., -0.],
       [ 0.,  1., -0.],
       [-0., -0.,  1.]])

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.

# Recover the singular values
U.transpose() @ input_matrix @ V.transpose()

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

array([[16.848,  0.   , -0.   ],
       [ 0.   ,  1.068,  0.   ],
       [-0.   , -0.   , -0.   ]])

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

[16.848  1.068  0.   ]

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.

import scipy.linalg

U, D, V = scipy.linalg.svd(input_matrix)

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

import sklearn.decomposition

svd = sklearn.decomposition.TruncatedSVD(n_components=3)

svd.fit(input_matrix)

print(svd.singular_values_)

This should print out the singular values as

[16.848  1.068  0.   ]

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

import torch.linalg

tensor_input = torch.tensor(input_matrix)

U, D, V = torch.linalg.svd(tensor_input)

print(D)

The resulting singular values are

tensor([1.6848e+01, 1.0684e+00, 1.0313e-07])

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

def dot(A, B):
    # Multiply matrix by a matrix
    if type(A[0]) != list:
        A = [A]

    if type(B[0]) != list:
        B = [[b] for b in B]

    ret = [[0] * len(B[0]) for i in range(len(A))]
    
    for row in range(len(ret)):
        for col in range(len(ret[0])):
            ret[row][col] = 0
            
            for i in range(len(B)):
                ret[row][col] += A[row][i] * B[i][col]
            
    # Return scalar if dimensions are 1x1
    if len(ret) == 1 and len(ret[0]) == 1:
        return ret[0][0]
    elif len(ret[0]) == 1:
        # Return vector if dimensions are 1xn
        return [r[0] for r in ret]

    # Return matrix
    return ret

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

def transpose(A):
    # Calculate matrix transpose
    if type(A[0]) != list:
        A = [A]

    rows = len(A)
    cols = len(A[0])

    B = [[0] * rows for i in range(cols)]
    
    for row in range(rows):
        for col in range(cols):
            B[col][row] = A[row][col]

    return B

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.

import random

# Squared input matrix (A*A^T)
input_squared = dot(input_matrix, transpose(input_matrix))

# Number of iterations
iterations = 100

# Number of SVDs to recover
N = min(len(input_squared), len(input_squared[0]))

# Return values 
# Left signular vectors
U = [[0] * len(input_squared[0]) for i in range(N)]

# Singular values
D = [0] * N

for n in range(N):
    # Randomly initialize search vector
    b = [random.random() for i in range(len(input_squared[0]))]

    dominant_svd = None
    for k in range(iterations):
        # Input matrix multiplied by b_k
        projection = dot(input_squared, b) 
                
        # Norm of input matrix multiplied by b_k
        norm = dot(projection, projection)**0.5

        # Calculate b_{k+1}
        b_next = [d / norm for d in projection]

        # Calculate the Rayleight quotient 
        dominant_svd = dot(b, projection) / dot(b, b)

        b = b_next

    D[n] = dominant_svd**0.5

    for i in range(len(b)):
        U[i][n] = b[i]


    # Remove current dimensions from the input before moving on 
    # the next singular value
    outer_product = [[0] * len(b) for j in range(len(b))]
    for i in range(len(b)):
        for j in range(len(b)):
            outer_product[i][j] = dominant_svd*b[i]*b[j]

    for i in range(len(input_squared)):
        for j in range(len(input_squared[0])):
            input_squared[i][j] -= outer_product[i][j]

Dinv = [[0] * N for i in range(N)]
for i in range(N):
    Dinv[i][i] = 1/D[i]

V = dot(Dinv, dot(transpose(U), input_matrix))

import pprint
print("Left singular vectors")
pprint.pprint(U)
print("Right singular vectors")
pprint.pprint(V)
print("Singular values")
pprint.pprint(D)

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

Left singular vectors
[[0.3863177031186115, 0.9223657800770559],
 [0.9223657800770583, -0.386317703118617]]
Right singular vectors
[[0.4286671335486261, 0.5663069188480351, 0.7039467041474441],
 [-0.8059639085893282, -0.11238241409663538, 0.581199080396058]]
Singular values
[9.508032000695724, 0.7728696356734851]

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.

Further resources

• NumPy SVD algorithm reference manual.
• SciPy SVD algorithm reference manual.
• Scikit-learn truncated SVD algorithm reference manual.
• PyTorch SVD algorithm reference manual.
• Singular value decomposition Wikipedia page.
• Power iteration method Wikipedia page.
• Eigenvalue decomposition in Python.