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)