This article will discuss QR Decomposition in Python. In previous articles we have looked at LU Decomposition in Python and Cholesky Decomposition in Python as two alternative matrix decomposition methods. QR Decomposition is widely used in quantitative finance as the basis for the solution of the linear least squares problem, which itself is used for statistical regression analysis.
One of the key benefits of using QR Decomposition over other methods for solving linear least squares is that it is more numerically stable, albeit at the expense of being slower to execute. Hence if you are performing a large quantity of regressions as part of a trading backtest, for instance, you will need to consider very extensively whether QR Decomposition is the best fit (excuse the pun).
For a square matrix $A$ the QR Decomposition converts $A$ into the product of an orthogonal matrix $Q$ (i.e. $Q^TQ=I$) and an upper triangular matrix $R$. Hence:
\begin{eqnarray*} A = QR \end{eqnarray*}There are a few different algorithms for calculating the matrices $Q$ and $R$. We will outline the method of Householder Reflections, which is known to be more numerically stable the the alternative Gramm-Schmidt method. I've outlined the Householder Reflections method below.
Note, the following explanation is an expansion of the extremely detailed article on QR Decomposition using Householder Reflections over at Wikipedia.
A Householder Reflection is a linear transformation that enables a vector to be reflected through a plane or hyperplane. Essentially, we use this method because we want to create an upper triangular matrix, $R$. The householder reflection is able to carry out this vector reflection such that all but one of the coordinates disappears. The matrix $Q$ will be built up as a sequence of matrix multiplications that eliminate each coordinate in turn, up to the rank of the matrix $A$.
The first step is to create the vector $\mathbb{x}$, which is the $k$-th column of the matrix $A$, for step $k$. We define $\alpha = -sgn(\mathbb{x}_k)(||\mathbb{x}||)$. The norm $||\cdot||$ used here is the Euclidean norm. Given the first column vector of the identity matrix, $I$ of equal size to $A$, $\mathbb{e}_1 = (1,0,...,0)^T$, we create the vector $\mathbb{u}$:
\begin{eqnarray*} \mathbb{u} = \mathbb{x} + \alpha \mathbb{e}_1 \end{eqnarray*}Once we have the vector $\mathbb{u}$, we need to convert it to a unit vector, which we denote as $\mathbb{v}$:
\begin{eqnarray*} \mathbb{v} = \mathbb{u}/||\mathbb{u}|| \end{eqnarray*}Now we form the matrix $Q$ out of the identity matrix $I$ and the vector multiplication of $\mathbb{v}$:
\begin{eqnarray*} Q = I - 2 \mathbb{v} \mathbb{v}^T \end{eqnarray*}$Q$ is now an $m\times m$ Householder matrix, with $Q\mathbb{x} = \left( \alpha, 0, ..., 0\right)^T$. We will use $Q$ to transform $A$ to upper triangular form, giving us the matrix $R$. We denote $Q$ as $Q_k$ and, since $k=1$ in this first step, we have $Q_1$ as our first Householder matrix. We multiply this with $A$ to give us:
\begin{eqnarray*} Q_1A = \begin{bmatrix} \alpha_1&\star&\dots&\star\\ 0 & & & \\ \vdots & & A' & \\ 0 & & & \end{bmatrix} \end{eqnarray*}The whole process is now repeated for the minor matrix $A'$, which will give a second Householder matrix $Q'_2$. Now we have to "pad out" this minor matrix with elements from the identity matrix such that we can consistently multiply the Householder matrices together. Hence, we define $Q_k$ as the block matrix:
\begin{eqnarray*} Q_k = \begin{pmatrix} I_{k-1} & 0\\ 0 & Q_k'\end{pmatrix} \end{eqnarray*}Once we have carried out $t$ iterations of this process we have $R$ as an upper triangular matrix:
\begin{eqnarray*} R = Q_t ... Q_2 Q_1 A \end{eqnarray*}$Q$ is then fully defined as the multiplication of the transposes of each $Q_k$:
\begin{eqnarray*} Q = Q^T_1 Q^T_2 ... Q^T_t \end{eqnarray*}This gives $A=QR$, the QR Decomposition of $A$.
To calculate the QR Decomposition of a matrix $A$ with NumPy/SciPy, we can make use of the built-in linalg library via the linalg.qr function. This is significantly more efficient than using a pure Python implementation:
import pprint
import scipy
import scipy.linalg # SciPy Linear Algebra Library
A = scipy.array([[12, -51, 4], [6, 167, -68], [-4, 24, -41]]) # From the Wikipedia Article on QR Decomposition
Q, R = scipy.linalg.qr(A)
print "A:"
pprint.pprint(A)
print "Q:"
pprint.pprint(Q)
print "R:"
pprint.pprint(R)
The output of the QR decomposition includes $A$, $Q$ and $R$. As a basic sanity check we can see that $R$ is in fact an upper triangular matrix:
A:
array([[ 12, -51, 4],
[ 6, 167, -68],
[ -4, 24, -41]])
Q:
array([[ 0.85714286, -0.39428571, -0.33142857],
[ 0.42857143, 0.90285714, 0.03428571],
[-0.28571429, 0.17142857, -0.94285714]])
R:
array([[ 14., 21., -14.],
[ -0., 175., -70.],
[ 0., -0., 35.]])
You aren't likely to ever need a pure Python implementation of QR Decomposition (homework notwithstanding), but I feel that it is helpful to gain an understanding of the Householder Reflections algorithm, so I have written my own implementation:
from math import sqrt
from pprint import pprint
def mult_matrix(M, N):
"""Multiply square matrices of same dimension M and N"""
# Converts N into a list of tuples of columns
tuple_N = zip(*N)
# Nested list comprehension to calculate matrix multiplication
return [[sum(el_m * el_n for el_m, el_n in zip(row_m, col_n)) for col_n in tuple_N] for row_m in M]
def trans_matrix(M):
"""Take the transpose of a matrix."""
n = len(M)
return [[ M[i][j] for i in range(n)] for j in range(n)]
def norm(x):
"""Return the Euclidean norm of the vector x."""
return sqrt(sum([x_i**2 for x_i in x]))
def Q_i(Q_min, i, j, k):
"""Construct the Q_t matrix by left-top padding the matrix Q
with elements from the identity matrix."""
if i < k or j < k:
return float(i == j)
else:
return Q_min[i-k][j-k]
def householder(A):
"""Performs a Householder Reflections based QR Decomposition of the
matrix A. The function returns Q, an orthogonal matrix and R, an
upper triangular matrix such that A = QR."""
n = len(A)
# Set R equal to A, and create Q as a zero matrix of the same size
R = A
Q = [[0.0] * n for i in xrange(n)]
# The Householder procedure
for k in range(n-1): # We don't perform the procedure on a 1x1 matrix, so we reduce the index by 1
# Create identity matrix of same size as A
I = [[float(i == j) for i in xrange(n)] for j in xrange(n)]
# Create the vectors x, e and the scalar alpha
# Python does not have a sgn function, so we use cmp instead
x = [row[k] for row in R[k:]]
e = [row[k] for row in I[k:]]
alpha = -cmp(x[0],0) * norm(x)
# Using anonymous functions, we create u and v
u = map(lambda p,q: p + alpha * q, x, e)
norm_u = norm(u)
v = map(lambda p: p/norm_u, u)
# Create the Q minor matrix
Q_min = [ [float(i==j) - 2.0 * v[i] * v[j] for i in xrange(n-k)] for j in xrange(n-k) ]
# "Pad out" the Q minor matrix with elements from the identity
Q_t = [[ Q_i(Q_min,i,j,k) for i in xrange(n)] for j in xrange(n)]
# If this is the first run through, right multiply by A,
# else right multiply by Q
if k == 0:
Q = Q_t
R = mult_matrix(Q_t,A)
else:
Q = mult_matrix(Q_t,Q)
R = mult_matrix(Q_t,R)
# Since Q is defined as the product of transposes of Q_t,
# we need to take the transpose upon returning it
return trans_matrix(Q), R
A = [[12, -51, 4], [6, 167, -68], [-4, 24, -41]]
Q, R = householder(A)
print "A:"
pprint(A)
print "Q:"
pprint(Q)
print "R:"
pprint(R)
The output from the Householder method implemented in pure Python is given below:
A:
[[12, -51, 4], [6, 167, -68], [-4, 24, -41]]
Q:
[[0.8571428571428571, 0.39428571428571435, -0.33142857142857135],
[0.4285714285714286, -0.9028571428571429, 0.034285714285714114],
[-0.28571428571428575, -0.17142857142857126, -0.942857142857143]]
R:
[[13.999999999999998, 21.00000000000001, -14.000000000000004],
[-5.506706202140776e-16, -175.00000000000003, 70.0],
[3.0198066269804245e-16, -3.552713678800501e-14, 35.000000000000014]]
You can see that we get the same answers as above in the SciPy implementation, albeit with a few more significant figures! One has to be extremely careful in numerical algorithms when dealing with floating point arithmetic, but that is a discussion for another day.
