I was going through my old list of to-do items and noticed that I made a note to implement QR decomposition, from scratch, using Python. If you have a matrix A and apply QR decomposition, you get matrices Q and R where Q * R = A.

There are three main QR decomposition algorithms: Householder, modified Gram-Schmidt, Givens. Comparisons are difficult but in general, Householder has the highest precision, while GS and Givens have lower precision but are dramatically simpler to implement. I decided to use Householder.

One of the many tricky details about QR decomposition is that there are two versions. One version returns “complete” Q and R matrices. The second version returns “reduced” Q and R. There are many algorithms for QR decomposition; two common algorithms are Householder and Gram-Schmidt.

After a bit of searching the Internet for resources, I got a nice demo up and running of QR decomposition using the Householder algorithm:

C:\VSM\MatrixPseudoInverseQR: python qr_demo.py

Begin QR decomposition using NumPy

Source A:
[[4.0000 7.0000 1.0000]
 [6.0000 0.0000 3.0000]
 [8.0000 1.0000 9.0000]
 [2.0000 5.0000 6.0000]
 [1.0000 5.0000 4.0000]]

===================

from-scratch Q (complete):
[[-0.3636  0.5998  0.6426 -0.1606 -0.2634]
 [-0.5455 -0.2854  0.2952  0.5362  0.4964]
 [-0.7273 -0.2677 -0.3760 -0.4394 -0.2549]
 [-0.1818  0.4692 -0.5091  0.6279 -0.3055]
 [-0.0909  0.5167 -0.3153 -0.3152  0.7252]]

from-scratch R (complete):
[[-11.0000  -4.6364 -10.0000]
 [ -0.0000   8.8603   2.2162]
 [ -0.0000   0.0000  -6.1716]
 [ -0.0000  -0.0000   0.0000]
 [ -0.0000  -0.0000  -0.0000]]

-------------------

from-scratch Q (reduced):
[[-0.3636  0.5998  0.6426]
 [-0.5455 -0.2854  0.2952]
 [-0.7273 -0.2677 -0.3760]
 [-0.1818  0.4692 -0.5091]
 [-0.0909  0.5167 -0.3153]]

from-scratch R (reduced):
[[-11.0000  -4.6364 -10.0000]
 [ -0.0000   8.8603   2.2162]
 [ -0.0000   0.0000  -6.1716]]

===================

. . . using np.linalg.qr() gives same results

End QR demo

I validated my from-scratch QR decomposition function by comparing with the NumPy np.linalg.qr() function to make sure I got the same results.

My motivation for all of this is to eventually implement Moore-Penrose pseudo-inverse using QR (as opposed to my usual technique of using SVD). But that’s an exploration for another day. Good fun.

I enjoy implementing functions from scratch because it gives me full visibility into the inner workings of the system.

Click to enlarge.

Demo program.

# qr_demo.py
# see rosettacode.org/wiki/QR_decomposition

import numpy as np

np.set_printoptions(precision=4, suppress=True,
  floatmode='fixed')

print("\nBegin QR decomposition from scratch Python demo ");

# -----------------------------------------------------------

def my_qr(A, mode):
  # mode = 'complete' or 'reduced'
  # this version defined only for n_rows gte n_cols
  m, n = A.shape
  k = np.minimum(m,n)  # always k == n
  R = A.copy()
  Q = np.eye(m)

  if m == n: end = n-1
  else: end = n
  for i in range(end):
    H = np.eye(m)

    # Householder routine follows
    a = R[i:, i]
    norm_a = np.copysign(np.linalg.norm(a), a[0])
    v = a / (a[0] + norm_a)
    v[0] = 1.0
    h = np.eye(len(a))
    V = np.matmul(v[:, None], v[None, :])
    h -= (2.0 / np.dot(v, v)) * V

    H[i:, i:] = h  # copy h into lower right
    Q = np.matmul(Q, H)
    R = np.matmul(H, R)

  # 'reduced' : returns Q, R with dims (m, k), (k, n)
  if mode == 'complete':
    return Q, R
  elif mode == 'reduced':
    return Q[0:m,0:k], R[0:k,0:n]

  return Q, R

# -----------------------------------------------------------

# 5x3
A = np.array([
[4, 7, 1],
[6, 0, 3],
[8, 1, 9],
[2, 5, 6],
[1, 5, 4]], dtype=np.float64)

print("\nSource A: ")
print(A)

print("\n=================== ")

Q,R = my_qr(A, 'complete')
print("\nfrom-scratch Q (complete): ")
print(Q)

print("\nfrom-scratch R (complete): ")
print(R)

print("\n------------------- ")

Q,R = my_qr(A, 'reduced')
print("\nfrom-scratch Q (reduced): ")
print(Q)

print("\nfrom-scratch R (reduced): ")
print(R)

print("\n=================== ")

Q, R = np.linalg.qr(A, mode='complete')
print("\nnumpy Q (complete): ")
print(Q)

print("\nnumpy R (complete): ")
print(R)

print("\n------------------- ")

Q, R = np.linalg.qr(A, mode='reduced')
print("\nnumpy Q (reduced): ")
print(Q)

print("\nnumpy R (reduced): ")
print(R)

print("\n=================== ")

print("\nEnd QR demo ")