HouseHolderQR.py

import numpy as np
def qr(a):
    m=np.shape(a)[0]
    n=np.shape(a)[1]
    r=a.copy()

    ident=np.identity(m)
    q=np.identity(m)
    for i in range(0,n):
        x = r[i:, i].copy()
        w=x.copy()
        w[0]+=np.linalg.norm(x,2)*sign(x[0])
        v=w/np.linalg.norm(w,2)
        h=np.identity(m-i)-2*np.outer(v,v)
        ident = np.identity(m)
        h2=ident[i:,i:]=h
        h2 = ident
        r=h2.dot(r)

        q=h2.dot(q)
    q=q.conj().T

    return q, r

def sign(x):
    if x>=0:
        return 1
    else:
        return -1
    end


a=np.random.rand(4, 3)*4
(q,r)=qr(a)

print('a', a)
print('r:', r)
print('q:', q)
print('norm(qr-a):', np.linalg.norm(q.dot(r) - a))