5
\$\begingroup\$

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  
\$\endgroup\$
4
  • 2
    \$\begingroup\$ Are you aware of numpy.linalg.qr? \$\endgroup\$ Jul 29, 2017 at 12:29
  • 1
    \$\begingroup\$ @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. \$\endgroup\$
    – Imago
    Jul 29, 2017 at 14:47
  • \$\begingroup\$ 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. \$\endgroup\$
    – hpaulj
    Jul 29, 2017 at 23:58
  • 1
    \$\begingroup\$ 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] \$\endgroup\$
    – hpaulj
    Jul 30, 2017 at 0:08

1 Answer 1

3
\$\begingroup\$

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))
\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.