QR法求解特征值特征向量

一 QR原理

理论依据：任意一个非奇异矩阵（满秩的方阵）A都可以分解为一个正交矩阵Q和一个上三角矩阵R的乘积，且当R对角元符号确定时，分解是唯一的。QR分解是一种迭代方法，迭代格式如下：
当Ak基本收敛到为上三角矩阵时，迭代完成，此时主对角元素就是特征值。

特别地：当A是对称阵的时候，Ak是对角阵Λ，Q=Qk-1Qk-2…Q1就是其正交特征向量矩,有Q^TAQ=Ak=Λ，即A正交对角化与Ak。

如何理解？我们看下图公式:

所以，QR迭代过程从数学的角度来想其实就是不断正交化的过程。

二 QR算法步骤

1.Householder变换进行QR分解

反射矩阵：任取单位向量w,反射矩阵H=E-2WW^T ,显然HH^T =E，H是正交阵

定理：任取两个模长相等的的向量x,y，一定存在一个反射矩阵H，使得Hx=y，
此时w=(x-y)/(|x-y|)（向量的差除以向量差的模）

应用：现在我们取矩阵的一列为x,m=|x|,y=m*[1,0,0,…0]^T 根据上面的定理求出H，使得Hx=y,是不是通过正交变化就把那一列化成了[m,0,0,0]^T ,这样就达到了将下三角元素全化为0的效果。看下图，举个例子来说明QR分解过程：

看懂上述过程就知道，Householder变换是利用了反射定理，经过n-1轮正交变换，将下三角元素全部化为0，从而得到上三角矩阵R，将所有H矩阵左乘运算再转置得到正交矩阵Q，即A=QR

我们看看QR分解的代码：

#QR分解
def qrSplit(A):
    n=A.shape[0]#A的维度
    Q=[[]]
    R=A
    for i in range(0,n-1):
        B=R
        if i!=0:
            #删除一行一列,得n-1阶子阵
            B=B[i:,i:]
        #取第一列向量
        x=B[:,0]
        #向量摸长
        m=np.linalg.norm(x)
        #生成一个模长为m，其余项为0的向量y
        y=[0 for j in range(0,n-i)]
        y[0]=m
        #计算householder反射矩阵
        #w = (x-y)/||x-y||
        w=