Householder similarity transformation of matrix in Python
''' d,c = householder(a). Householder similarity transformation of matrix [a] to the tridiagonal form [c\d\c]. p = computeP(a). Computes the acccumulated transformation matrix [p] after calling householder(a). ''' from numpy import dot,diagonal,outer,identity from math import sqrt def householder(a): n = len(a) for k in range(n-2): u = a[k+1:n,k] uMag = sqrt(dot(u,u)) if u[0] < 0.0: uMag = -uMag u[0] = u[0] + uMag h = dot(u,u)/2.0 v = dot(a[k+1:n,k+1:n],u)/h g = dot(u,v)/(2.0*h) v = v - g*u a[k+1:n,k+1:n] = a[k+1:n,k+1:n] - outer(v,u) \ - outer(u,v) a[k,k+1] = -uMag return diagonal(a),diagonal(a,1) def computeP(a): n = len(a) p = identity(n)*1.0 for k in range(n-2): u = a[k+1:n,k] h = dot(u,u)/2.0 v = dot(p[1:n,k+1:n],u)/h p[1:n,k+1:n] = p[1:n,k+1:n] - outer(v,u) return p