Linear equations with Python: the QR decomposition

See the first article in this series Solving linear equations using matrices and Python.

In this second article on methods for solving systems of linear equations using Python, we will see the QR Decomposition method.

This method is very similar to the LU decomposition. The equation to be solved is of the form Ax = B. In this particular case, the matrix A = QR, where Q is an orthogonal matrix and R is an upper triangular matrix.

$Ax=B$

$A=QR$

What are the properties of an orthogonal matrix?

It is a square matrix;
Multiplying Q for its transpose, we obtain the identity matrix;
$QQ^{T}=Q^{T}Q=I$
The inverse of an orthogonal matrix is equal to its transpose.
$Q^{T}=Q^{-1}$

Moreover, from matrix algebra we remember that if A has n linearly independent columns, then the first n columns of Q form an orthonormal basis for the space of columns of A. In general, the first k columns of Q form an orthonormal basis for the subspace of the first k columns of A, for each $1\leq k\leq n$ .

An identity matrix is a square matrix in which the elements of the main diagonal are 1 and all other elements are zeros.

Now you can restate the initial problem $Ax=B$ this way:

$QRx=B$

$Rx=Q^{-1}B$

$or$

$Rx=Q^{T}B$

We use the same variables as the LU decomposition example and using the function qr of scipy.linalg we calculate Q and R. Finally we assign to the variable y the value of $BQ^{T}$ .


""" QR decomposition with scipy """
import scipy.linalg as linalg
import numpy as np
# same matrix A and B as in LU decomposition
A = np.array([
[2., 1., 1.],
[1., 3., 2.],
[1., 0., 0]])
B = np.array([4., 5., 6.])
Q, R = linalg.qr(A) # QR decomposition with qr function
y = np.dot(Q.T, B) # Let y=Q'.B using matrix multiplication
x = linalg.solve(R, y) # Solve Rx=y
print x

Note that $Q.T$ is simply the transpose of Q, which is in turn equal to the inverse of Q.