QR算法求特征值python实现--基于Householder变换

import numpy as np

#定义householder变换

def householder(x,i=0):

    x=x.reshape(len(x),1)

    ei=np.zeros((len(x),1))

    ei[i]=1

    y=np.linalg.norm(x,ord=2)*ei

    if x[i]>0:

        w=(x+y)/np.linalg.norm(x+y)

    else:

        w=(x-y)/np.linalg.norm(x-y)

    H=np.eye(len(x))-2*np.dot(w,w.T)

    return H

#进行QR分解

def QRFact(A):

    dim=len(A)

    Ri=A.copy()

    for i in range(dim-1):

        x=Ri[i:,i]

        Hi=householder(x)

        Ri[i:,i:]=np.dot(Hi,Ri[i:,i:])

        Qi=np.eye(dim)

        Qi[i:,i:]=Hi

        if i==0:

            Q=Qi

        else:

            Q=np.dot(Qi,Q)

    D=np.asmatrix(np.diag(np.where(np.diag(Ri)<0,-1,1)))

    R=D*np.asmatrix(Q.T)*D.T

    return Q,R

#定义迭代求特征值

def eig_QR(A,eps=1e-6):

    Ak=A.copy()

    flag=1

    n=0

    while flag:

        Ak0=Ak.copy()

        Qk,Rk=QRFact(Ak)

        Ak=Rk*Qk

        if(np.sum(np.diag(np.abs(Ak-Ak0)))

            flag=0

        n=n+1

    return np.diag(Ak),n

#主程序

if __name__=='__main__':

    A=np.array([[3,2],[4,5]],dtype='float')

    eig=eig_QR(A,1e-6)

    print('eig:',eig)

    eigval, vec = np.linalg.eig(A)

    print(eigval)

本文来自互联网用户投稿，文章观点仅代表作者本人，不代表本站立场，不承担相关法律责任。如若转载，请注明出处。 如若内容造成侵权/违法违规/事实不符，请点击【内容举报】进行投诉反馈！