8.QR分解的python實現

import numpy as np
import math

#直到主對角線上的值變化很小時，結束循環
def is_same(a,b):
    print(a)
    print(b)
    n = len(a)
    for i in range(n):
        if(math.fabs(a[i]-b[i]) > 1e-9):
            return False
    return True

if __name__ == '__main__':
    a = np.array([0.65,0.28,0.02,0.15,0.67,0.18,0.12,0.36,0.52])
    n = int(math.sqrt(len(a)))
    a = a.reshape((n,n))
    #計算特征值和特征向量
    value,v = np.linalg.eig(a)
    # print("特征值:",value)
    # print("特征向量:",v)

    times = 0
    #diag矩陣的對角線元素
    while(True):
        # v1保存的是前一個矩陣的對角線元素
        v1 = np.diag(a)
        print(a)
        print(v1)
        #QR分解計算新的a
        q,r = np.linalg.qr(a)
        a = np.dot(r,q)
        times += 1
        if( is_same(np.diag(a),v1) ):
            break
        # print("正交陣:",q)
        # print("三角陣:",r)
        # print("近似陣",a)
        # print("==="*30)

    print("次數:",times," 近似值:",np.diag(a))
    print("精准特征值:",value)