Implementation of the Householder Transformation

Question

I implemented the Householder transformation in Python, so that I can later use it in a QR decomposition. Unfortunately I haven't found a good concise source for reading up on the algorithm.

I am not really satisfied with my code - I mostly dislike its readability due to the number of various computations. I also find my usage of the identity matrices \$I\$ and \$J\$ inefficient. So I would be happy about any suggestions on improving this implementation.

def houseHolder2(m):    
    l = len(m)
    for i in range(0, l - 1):
        I = np.eye(l - i)
        ai = m[:, i] 
        abs = np.linalg.norm(ai[i:l]) 
        sign = np.sign(m[i, i])
        vi = np.array([ai[i:l]]).T + sign * abs * np.array([I[:, 0]]).T
        Qi = I - 2 * (vi @ vi.T) / (vi.T @ vi)
        J = np.eye(l)
        J[i:l, i:l] = Qi
        m = J@m
    return m  

Answer 1

This is an old question, but since people might get stuck here:

The basic reason for the inefficiency is that you try to actually compute \$Q\$ explicitly, which is not required when you want to perform the QR algorithm or try to solve a linear equation. All you need is the ability to compute \$Qb\$ and \$Q^tb\$ for vectors b. Thus, there's no need to compute the identity matrices.

The common procedure is to store the vectors \$v_i\$ and then compute the product \$Qb\$ using these vectors. I'm including my version of the householder algorithm with the corresponding function for solving the equation. The code works for non-quadratic matrices. I rewrote it from memory for a class, but Golub-van Loan should definitely have a better version.

def householder(A):
    (n,m)=A.shape
    p=min(n,m)
    alpha=np.zeros(m)
    for j in range(0,p):
        alpha[j]=np.linalg.norm(A[j:,j])*np.sign(A[j,j])
        if (alpha[j]!=0):
            beta=1/math.sqrt(2*alpha[j]*(alpha[j]+A[j,j]))
            A[j,j]=beta*(A[j,j]+alpha[j])
            A[j+1:,j]=beta*A[j+1:,j]
            for k in range(j+1,m):
                gamma=2*A[j:,j].dot(A[j:,k])
                A[j:,k]=A[j:,k]-gamma*A[j:,j]
    return A,alpha

def loese_householder(H,alpha,b):
    (n,m)=H.shape
    b=b.copy()
    x=np.zeros(n)
    # b=Q^t b.
    for j in range(0,n):
        b[j:]=b[j:]-2*(H[j:,j].dot(b[j:]))*H[j:,j]
    # Auflösen von Rx=b.
    for i in range(0,n):
        j=n-1-i
        b[j]=b[j]-H[j,j+1:].dot(x[j+1:])
        x[j]=-b[j]/alpha[j]
    return x


n=500
A=np.random.random([n,n])
b=np.random.random(n)
H,alpha=householder(A.copy())
x=loese_householder(H,alpha,b)
print(np.linalg.norm(A.dot(x)-b))