Householder similarity transformation of matrix in Python,,''' d,c = ho


''' 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,identityfrom math import sqrtdef 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

评论关闭