# Implementation of the Householder Transformation

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

• Are you aware of numpy.linalg.qr? Jul 29, 2017 at 12:29
• @MrGrj, ye, I think I have seen something of sorts. However it is not my goal to solve a specific problem or computation, but rather just improve my understanding how algorithm works. Being able to write them down in clean code should be also part of the process imo. Jul 29, 2017 at 14:47
• What's a typical m? From an efficiency view point, the real question is how many times are you repeating the loop, a few times, or hundreds. Jul 29, 2017 at 23:58
• I think you can simplify the indexing of I with I[:,[0]]; Or the whole line with vi=(ai[i:l] + sign * abs * I[:,0])[:,None] Jul 30, 2017 at 0:08

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))