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
开心洋葱 , 版权所有丨如未注明 , 均为原创丨未经授权请勿修改 , 转载请注明Householder similarity transformation of matrix in Python！
