Python

How to Calculate QR Decomposition in Python

Jake @Scicoding

May 17, 2023 • 7 min read

QR decomposition is one of the essential matrix factorization techniques. It simplifies numerous computational tasks making it a necessary ingredient in many engineering and scientific algorithms and applications.

One of the key advantages of QR decomposition is its role in solving linear least-squares problems. It transforms the original problem into a more manageable form and enables more accurate and efficient solutions. These are vital in data fitting and optimization. Moreover, QR decomposition plays an important role in eigenvalue and eigenvector computations.

QR decomposition holds great importance in science and engineering applications due to its versatility in simplifying complex problems, enhancing computational efficiency, and enabling deeper insights into various systems.

In this post, we will dive into using QR decomposition with Python. You'll learn how to perform QR factorization using the most important Python packages for mathematical computation: NumPy and SciPy. In addition, you learn two methods to implement a custom algorithm for QR decomposition without having to use anything else but the Python core modules.

Applications of QR decomposition

QR factorization is an inherent component of many algorithms. It is applied in engineering and science, showcasing its versatility and importance across a wide range of disciplines. Its ability to simplify complex problems and enhance computational efficiency makes it an invaluable tool in our journey toward innovation and discovery!

Eigenvalue decomposition: the QR algorithm or QR iteration is a practical algorithm for calculating the eigenvalue decomposition. That is, it provides both the eigenvalues and eigenvectors.
Singular value decomposition: Similar to the eigenvalue decomposition, a variant of the QR iteration can be used to calculate the singular value decomposition (SVD).
Data compression: QR decomposition can be used for efficient rank-deduction, and, as such, is useful for data compression applications.
Signal processing: QR factorization can be used for various signal processing applications, such as linear prediction, adaptive filtering, and spectral analysis.
Principle component analysis (PCA): Similarly, to rank reduction and SVD, QR decomposition is used within many PCA algorithms.
Linear systems and least squares problems: QR factorization provides a convenient way of solving linear systems. Using the QR decomposition, consider a linear system defined by matrix \(\mathbf{A} = \mathbf{QR}\). Then, we have \(\mathbf{Ax} = \mathbf{b} \Leftrightarrow\) \(\mathbf{QRx} = \mathbf{b} \Leftrightarrow\) \(\mathbf{Rx} = \mathbf{Q}^\mathrm{H}\mathbf{b}\). Since \(\mathbf{R}\) is triangular, we can solve \(\mathbf{x}\) by simple backwards substitution.

Test data

Before getting into the actual methods of calculating QR decomposition in Python, we can create common test data. This way we can verify that our different schemes provide similar results.

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

Test input matrix for QR decomposition

QR Decomposition using NumPy

NumPy has readily available support for QR decomposition. The interface for the matrix factorization can be accessed using the linear algebra module np.linalg. The QR decomposition can be calculated using np.linalg.qr.

Using np.linalg.qr is very straightforward. There are only two parameters. The first one is the matrix to be factorized. The second one indicates mode that is used to control the factorization. By default mode='reduced'. This means that np.linalg.qr returns matrix \(\mathbf{Q} \in \mathcal{C}^{m \times k}, \mathbf{R} \in \mathcal{C}^{k \times n},\) where \((m ,n)\) are the dimensions of the input matrix and \(k = \min(m, k)\).

mode	Q dimensions	R dimensions
reduced	m x k	k x n
complete	m x m	m x n
r	-	k x n

In addition, there is mode=raw, which return h and tau parameters from the internally used Householder reflections algorithm.

Now, let's have a look at how to use np.linalg.qr in practice.

import numpy as np

Q, R = np.linalg.qr(input_matrix)

Calculating QR Decomposition with NumPy

Next, we'll check that the results look as expected. We should be able to reconstruct the input matrix by calculating \[\mathbf{Q}\mathbf{R}.\]

# Reconstructing the input matrix from the QR decomposition
Q @ R

Reconstruct the input matrix

The reconstructed matrix matches with the input matrix as expected.

array([[1., 2., 3.],
       [4., 5., 6.],
       [7., 8., 9.]])

Reconstructed input matrix

We can also verify that \(\mathbf{Q}\) is orthonormal by calculating \(\mathbf{Q}^H \mathbf{Q}\)

# Check orthonormality
Q.T.conjugate() @ Q

Check orthonormality

This returns the identity matrix as expected

array([[ 1.00000000e+00, -2.48218629e-16, -2.62094065e-17],
       [-2.48218629e-16,  1.00000000e+00,  2.82872043e-17],
       [-2.62094065e-17,  2.82872043e-17,  1.00000000e+00]])

The identity matrix returned by the orthonormality check

Finally, we can check that \(\mathbf{R}\) is actually upper triangular

array([[-8.12403840e+00, -9.60113630e+00, -1.10782342e+01],
       [ 0.00000000e+00,  9.04534034e-01,  1.80906807e+00],
       [ 0.00000000e+00,  0.00000000e+00, -8.88178420e-16]])

Upper triangular matrix R

QR Decomposition using SciPy

The second most popular numerical computation library scipy also provides an interface to QR decomposition. It is included in the linear algebra package scipy.linalg. The actual method is scipy.linalg.qr. With respect to NumPy, SciPy provides some additional flexibility.

It works pretty much the same as np.linalg.qr. You just provide it with an input matrix and it will return the \(\mathbf{Q}\) and \(\mathbf{R}\) matrices.

SciPy uses a similar mode parameter to adjust the dimensions of the returned matrices. In addition, there is mode=raw, which returns the tau matrix for insight into the interworkings of the underlying algorithm.

mode	Q dimensions	R dimensions
full	m x m	m x k
economic	m x k	k x n
r	-	k x n

Using scipy.linalg.qr is as easy as the corresponding NumPy method.

import scipy

Q, R = scipy.linalg.qr(input_matrix)

Calculating QR Decomposition with SciPy

Simple! Now, we can validate the results

# Reconstructing the input matrix from the QR decomposition
Q @ R

Reconstruct the input matrix from SciPy QR decomposition

This of course returns the same input matrix as we had defined for our test data.

Pure Python algorithms for QR decomposition

We have now gone through two essential numerical computation packages and seen how to perform QR decomposition on those. If you happen to be in a situation, where external libraries are not an option. Then, implementing QR decomposition from scratch might be the only option. Luckily, the algorithms are not that complex and these do not require advanced linear algebra preliminaries matrix multiplication and transpose are enough.

Gram-Schmidt process

The Gram-Schmidt process is an algorithm for orthogonalizing a set of vectors in an inner product space, most commonly the Euclidean space \( \mathbb{R}^n \). Given a set of vectors, the process iteratively orthonormalizes the vectors to create an orthogonal (or orthonormal) basis for the span of the original vectors.

Let's break down the steps of the Gram-Schmidt process for QR decomposition:

Given a matrix \( \mathbf{A}) where \( a_i \) are the column vectors of the matrix, we start by normalizing the first vector:

\[ q_1 = \frac{a_1}{|a_1|}. \]

For each subsequent vector \( a_j \), we subtract the projection of \( a_j \) onto each of the \( q_i \) for \( i \lt j \), and then normalize:

\[ v_j = a_j - \sum_{i=1}^{j-1} (q_i^T a_j) q_i \]

\[q_j = \frac{v_j}{|v_j|} \]

Continue this process for all column vectors in \( \mathbf{A} \).

The \( \mathbf{Q} \) matrix is then formed by the orthonormal vectors \( q_i \) and the \( \mathbf{R} \) matrix is formed by the coefficients of the projections. Specifically, \( \mathbf{R} \) is an upper triangular matrix where each element \( r_{ij} \) is given by the dot product \( q_i \cdot a_j \) for \( i \leq j \) and \( r_{ij} = 0 \) for \( i > j \).

Here is a Python function for performing the QR decomposition using the Gram-Schmidt process:

def qr_decomposition(A):
    A_ = transpose(A)
    n = len(A)
    Q = [[0] * n for _ in range(n)]
    R = [[0] * n for _ in range(n)]

    for i in range(n):
        # Begin the Gram-Schmidt process
        v = A_[i]
        
        for j in range(i):
            R[j][i] = dot_product(Q1[j], A_[i])
            v = vector_subtract(v, scalar_multiply(Q[j], R[j][i]))
        
        R[i][i] = vector_norm(v)
        Q[i] = scalar_multiply(v, 1 / R[i][i])

    return transpose(Q1), R

def transpose(U):
    return [[row[i] for row in U] for i in range(len(U[0]))]

def dot_product(u, v):
    return sum(ui * vi for ui, vi in zip(u, v))

def vector_subtract(u, v):
    return [ui - vi for ui, vi in zip(u, v)]

def scalar_multiply(u, c):
    return [ui * c for ui in u]

def vector_norm(v):
    return sum(vi ** 2 for vi in v) ** 0.5

Q, R = qr_decomposition(A)

In this code, qr_decomposition is the main function, which takes a matrix A as input and returns the matrices Q and R.

These helper functions perform the following operations:

dot_product(u, v): This function calculates the dot product of two vectors u and v.
vector_subtract(u, v): This function subtracts vector v from vector u.
scalar_multiply(u, c): This function multiplies vector u by a scalar c.
vector_norm(v): This function calculates the norm (length) of vector v.

Inside the qr_decomposition function, we initialize the Q and R matrices with zero values. Then, for each column vector in A, we perform the Gram-Schmidt process. The R[j][i] value is calculated as the dot product of Q[j] and A[i], and subtracted from A[i] to get the next orthogonal vector v. We then normalize v to get Q[i], and the norm of v becomes the R[i][i] value.

Please note that Python indexing starts at 0, so Q[0] and R[0][0] in the code correspond to \(q_1\) and \(r_{11}\) in the mathematical description, and so on.

Summary

In the article, various methods for performing QR decomposition in Python were explored, leveraging both built-in scientific computing libraries and custom implementation. The NumPy and SciPy libraries, known for their efficient and robust numerical routines, are employed to perform QR decomposition, demonstrating how these libraries facilitate complex mathematical operations with straightforward, readable code.

Moreover, we delved into the intricate details of the Gram-Schmidt process, a foundational algorithm for QR decomposition. A custom Python function implementing this algorithm is provided, offering insight into the underlying mechanics of QR decomposition. By comparing these approaches, we presented a comprehensive understanding of QR decomposition, from high-level library functions to the nitty-gritty of algorithmic implementation.