Implementation

import numpy as np

def QRstep(A):
    # modify A in place and return H
    m = A.shape[0]  # number of rows of A
    w = np.zeros(m)
    z = A[:,0] # first column of A
    w[0] = np.linalg.norm(z)
    if  z[0] >= 0: w[0] *=-1  # choose - to avoid possible bad r.o. error
    v = w - z
    Pv = np.outer(v,v)/np.dot(v,v)
    H = np.eye(m) - 2*Pv
    #A = H@A  # incorrect because creates a new matrix in memory
              # and points the name A at it, so the memory location referred to by A
              # in the calling routine is not modified.
    A[:,:] = H@A  # overwrite HA on A
    return H

def QR(A):
    # modify A in place and return Q
    m,n = A.shape  # number of rows and columns of A
    H = np.eye(m)
    for i in range(n):
        # construct the "Hi"
        Htildei = QRstep(A[i:,i:])
        Hi = np.eye(m)
        Hi[i:,i:] = Htildei
        H = Hi@H
    return H.T   # return Q

Note: an efficient implementation would not actually form the H matrices, but would work only with the v vectors.

# Let's test it!

A = np.random.randn(5,3)
Asaved = A.copy()

Q = QR(A)  # will modify A in place

display(Asaved)
display(Q)
display(A) # this is the transformed A, called R
           #                                   0

# first check that Q is orthogonal

import matplotlib.pyplot as plt
plt.imshow((Q.T)@Q,interpolation='nearest',cmap='gray')
plt.colorbar()

array([[ 0.16291236,  1.24942785, -3.15583016],
       [-1.37319029,  0.26587118,  0.57615065],
       [ 0.41994405,  0.50289645, -2.06341379],
       [-0.83800519,  0.029154  ,  1.0927124 ],
       [-3.4505053 , -0.88767803, -0.85167823]])

array([[-0.04249549, -0.85403932, -0.25619827, -0.01392667,  0.45053242],
       [ 0.35819503, -0.38964036,  0.56830581, -0.48935686, -0.3967809 ],
       [-0.1095419 , -0.29142011, -0.42393359,  0.30367714, -0.7944403 ],
       [ 0.21859265, -0.14424203,  0.51033369,  0.81672953,  0.06264092],
       [ 0.90006015,  0.11430564, -0.41380042,  0.03269523,  0.06727684]])

array([[-3.83363855e+00, -8.05540320e-01,  3.88103802e-02],
       [ 9.90766715e-17, -1.42288064e+00,  2.81706510e+00],
       [ 3.81724700e-18, -4.16484782e-18,  2.92077116e+00],
       [-2.52188336e-18, -5.10047312e-17, -1.27786966e-17],
       [ 2.81742418e-16,  7.91706409e-17, -3.63185212e-16]])

Observe that the last matrix is indeed upper-triangular to machine precision, and Q^TQ is the identity matrix to machine precision:

Now check that QR is really a factorization of A:

# check that QR = A
Q@A - Asaved

array([[-2.22044605e-16,  2.22044605e-16, -4.44089210e-16],
       [-4.44089210e-16, -1.66533454e-16, -5.55111512e-16],
       [-2.77555756e-16, -1.11022302e-16,  4.44089210e-16],
       [ 1.11022302e-16,  3.81639165e-17,  0.00000000e+00],
       [ 0.00000000e+00, -2.22044605e-16,  3.33066907e-16]])

Yes, the difference is zero to machine precision. 🙂