精华内容
下载资源
问答
  • QR分解求矩阵全部特征值

    万次阅读 2014-11-08 18:57:28
    QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式   将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值。  QR方法的计算步骤如下:    下面就依次进行介绍。  一. ...

    QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式

                                    

    将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值。

           QR方法的计算步骤如下:

              

              下面就依次进行介绍。


           一. 将一般矩阵化为上Hessenberg阵

            1.定义                                      

                     一个矩阵如果满足i>j+1时aij=0,则将这个矩阵成为上Hessenberg阵。上Hessenberg阵

             的形式如下:

                     

               2. Householder变换将一般矩阵转化为上Hessnberg阵

                           首先,选取Householder矩阵H1,使得经过H1相似变换后的矩阵H1AH1的元素a21下面的

                   元素全部为0,即a31, a41, ....., am1均为0,H1取如下形式

                                               

                        其中 为n-1阶HouseHolder矩阵。然后选取Householder矩阵H2,使得经过H2相似变换

                之后的矩阵H2(H1AH1)H2第二列中a32下面的a42, ....am2均为0。如此进行n-2次,可以构造

               n-2个householder矩阵H1,H2, Hn-2,使得 Hn-2....H2H2AH1H2....Hm-2 = H(H为上Hessenberg矩阵)。

                        对于一个n*m的矩阵A,第col次的H可以这样构造求得(col从0开始):

                          

                    其中,I为n*n的单位矩阵, v'表示矩阵v的转置, sign(x0)表示x0的符号的相反数( 当x0>0时sign=-1,当x<=0时为1),

                    ||x||表示向量x的长度, col等于所求的上hessenberg矩阵的序号,从0开始


          二. 用Givens变换对上hessnberg矩阵作QR分解

                      

                      

                           

                      

                              

                            此时有  H = R21' * R32' * ... * Rn(n-1)'R = QR。

                            多次计算H,直到H的变化小于一个较小的阈值时,停止迭代,此时H主对角线上的元素

                即为矩阵A的全部特征值。

                            下面举个例子来说明求解矩阵的全部特征值的过程。

                                     求矩阵的全部特征值

                                首先将A化成上hessenberg阵,取

                                  x = [0, 6, 4], 则 ||x|| = = 

                           则 w = [0, , 0] , v = w + 1 * x = [0, 6+, 4]

                           则 p = v*v'/v'*v =          

                           于是 H1 = I - 2*p =  

                                    所以 H = H1AH1 =

                                    H即为与A相似的上hessenberg矩阵。将H进行QR分解

                                    

                                                                                                             

                                         

                             

                  这个程序的完整代码可以到这里下载,http://download.csdn.net/detail/xxc1605629895/6473181

    展开全文
  • (转)QR分解求矩阵全部特征值

    千次阅读 2015-01-09 21:24:49
    QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式   将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值。  QR方法的计算步骤如下:    下面就依次进行介绍。  一. ...

    QR算法求矩阵全部特征值的基本思想是利用矩阵的QR分解通过迭代格式

                                    

    将A=A1化成相似的上三角阵,从而求出矩阵A的全部特征值。

           QR方法的计算步骤如下:

              

              下面就依次进行介绍。


           一. 将一般矩阵化为上Hessenberg阵

            1.定义                                      

                     一个矩阵如果满足i>j+1时aij=0,则将这个矩阵成为上Hessenberg阵。上Hessenberg阵

             的形式如下:

                     

               2. Householder变换将一般矩阵转化为上Hessnberg阵

                           首先,选取Householder矩阵H1,使得经过H1相似变换后的矩阵H1AH1的元素a21下面的

                   元素全部为0,即a31, a41, ....., am1均为0,H1取如下形式

                                               

                        其中 为n-1阶HouseHolder矩阵。然后选取Householder矩阵H2,使得经过H2相似变换

                之后的矩阵H2(H1AH1)H2第二列中a32下面的a42, ....am2均为0。如此进行n-2次,可以构造

               n-2个householder矩阵H1,H2, Hn-2,使得 Hn-2....H2H2AH1H2....Hm-2 = H(H为上Hessenberg矩阵)。

                        对于一个n*m的矩阵A,第col次的H可以这样构造求得(col从0开始):

                          

                    其中,I为n*n的单位矩阵, v'表示矩阵v的转置, sign(x0)表示x0的符号的相反数( 当x0>0时sign=-1,当x<=0时为1),

                    ||x||表示向量x的长度, col等于所求的上hessenberg矩阵的序号,从0开始


          二. 用Givens变换对上hessnberg矩阵作QR分解

                      

                      

                           

                      

                              

                            此时有  H = R21' * R32' * ... * Rn(n-1)'R = QR。

                            多次计算H,直到H的变化小于一个较小的阈值时,停止迭代,此时H主对角线上的元素

                即为矩阵A的全部特征值。

                            下面举个例子来说明求解矩阵的全部特征值的过程。

                                     求矩阵的全部特征值

                                首先将A化成上hessenberg阵,取

                                  x = [0, 6, 4], 则 ||x|| = = 

                           则 w = [0, , 0] , v = w + 1 * x = [0, 6+, 4]

                           则 p = v*v'/v'*v =          

                           于是 H1 = I - 2*p =  

                                    所以 H = H1AH1 =

                                    H即为与A相似的上hessenberg矩阵。将H进行QR分解

                                    

                                                                                                             

                                         

                             

                  这个程序的完整代码可以到这里下载,http://download.csdn.net/detail/xxc1605629895/6473181


    展开全文
  • -s V2[i,i+1] = s V2[i+1, i+1] = c ai2 = V2*ai2 ai2 = ai2*V1.T*V2.T print('迭代88次后得:') print(ai2) print('矩阵特征值为{:.7},{:.7},{:.7}'.format(ai2[0, 0],ai2[1, 1],ai2[2, 2]))
    import numpy as np
    import math
    
    
    def sgn1(x):
        if x > 0:
            return 1
        elif x == 0:
            return 0
        elif x < 0:
            return -1
    
    
    ai2 = np.mat([[-1, 2, 1],
                [2, -4, 1],
                [1, 1, -6]], dtype=float)
    n = ai2.shape[0]
    for i in range(n-2):
        c = -1*sgn1(ai2[i+1, i])*np.sqrt(np.sum(np.multiply(ai2[i+1:n, i], ai2[i+1:n, i])))
        lou = np.sqrt(2*c*(c-ai2[i+1, i]))
        l = ai2[i+1:n, i]
        l1 =l.copy()
        l1[0] = l1[0]-c
        b = l1/lou
        a = np.mat(np.zeros((i+1, 1)))
        u = np.vstack((a, b))
        I = np.mat(np.eye(n))
        H = I-2*u*u.T
        ai2 = H*ai2*H.T
    err = 1
    ai3 = ai2.copy()
    for t in range(88):
        i = 0
        sita = math.atan(ai2[i+1, i]/ai2[i, i])
        c = math.cos(sita)
        s = math.sin(sita)
        V1 = np.mat(np.eye(n))
        V1[i,i] = c
        V1[i+1,i] = -s
        V1[i,i+1] = s
        V1[i+1, i+1] = c
        ai2 = V1*ai2
        i = 1
        sita = math.atan(ai2[i+1, i]/ai2[i, i])
        c = math.cos(sita)
        s = math.sin(sita)
        V2 = np.mat(np.eye(n))
        V2[i,i] = c
        V2[i+1,i] = -s
        V2[i,i+1] = s
        V2[i+1, i+1] = c
        ai2 = V2*ai2
        ai2 = ai2*V1.T*V2.T
    print('迭代88次后得:')
    print(ai2)
    print('矩阵的特征值为{:.7},{:.7},{:.7}'.format(ai2[0, 0],ai2[1, 1],ai2[2, 2]))
    
    
    

    在这里插入图片描述

    展开全文
  • QR方法、幂法和反幂法求矩阵特征值 本文介绍求任意矩阵全部特征值QR方法,求部分特征值和特征向量的幂法,反幂法。
  • 矩阵特征值的数值求解方法中,一般有部分特征值和特征向量的幂法和反幂法,以及全部特征值QR分解方法。在此关注QR分解方法QR分解,是把矩阵分解为正交矩阵和三角矩阵的乘积,其中关键点在于正交矩阵的...

    前言

    在上一篇的PCA中使用了LinearAlgebra中封装的函数求取的特征值和特征向量。在此学习一下特征值的数值求解算法,主要参考[1]。矩阵特征值的数值求解方法中,一般有求部分特征值和特征向量的幂法和反幂法,以及求取全部特征值的QR分解方法。在此关注QR分解方法。注意内容较长,总体思路和方法来自参考文献[1]。

    • QR分解,是把矩阵分解为正交矩阵和三角矩阵的乘积,其中关键点在于正交矩阵的求法。

    • 正交矩阵的求解可以使用施密特正交化进行实现,需要理解施密特正交化的原理和思路。

    最后回溯:由施密特正交化->>QR分解->>矩阵特征值->>特征向量矩阵 。

    1. 施密特正交

    施密特方法是对一组向量进行正交化。给定一组输入的m维向量,目的找出正交坐标系统,获取由这些向量张成的空间。
    给定n个线性无关的向量,找到n个彼此垂直的单位向量。

    (1) 利用施密特正交求出正交矩阵Q

    假设存在nn个线性无关的mm维向量
    A1,A2,,AnA_{1},A_{2},\cdots,A_{n}

    • 11步:选定初始向量,归一化。

    y1=A1y_{1}=A_{1}

    p1=y1y12p_{1}=\dfrac{y_{1}}{||y_{1}||_{2}}

    • 22步:下一个向量例如A2A_{2},减去该向量在所有已知的正交向量集(p1pi,)(p_{1},p_{i},\cdots)上的投影分量,再进行归一化得到当前的正交向量。

    投影分量如何计算?利用向量的内积cos(θ)=<a,b>a2b2cos(\theta)=\dfrac{<a,b>}{||a||_{2}||b||_{2}},此处符号<a,b><a,b>为向量内积。

    由于正交向量pip_{i}的模长,也即是2范数为1,此时AjA_{j}在向量pip_{i}上的投影系数为<Aj,pi>Aj2\dfrac{<A_{j},p_{i}>}{||A_{j}||_{2}},也即是piTAjp_{i}^{T}A_{j},注意此为一个数值标量。进一步可得到,AjA_{j}在向量pip_{i}上的投影分量为pi(piTAj)p_{i}(p_{i}^{T}A_{j})

    y2=A2p1(p1TA2)y_{2}=A_{2}-p_{1}(p_{1}^{T}A_{2})

    p2=y2y22p_{2}=\dfrac{y_{2}}{||y_{2}||_{2}}

    • 以此类推,第kk步,

    yk=Akp1(p1TAk)pk1(pk1TAk)y_{k}=A_{k}-p_{1}(p_{1}^{T}A_{k})-\cdots -p_{k-1}(p_{k-1}^{T}A_{k})

    pk=ykykkp_{k}=\dfrac{y_{k}}{||y_{k}||_{k}}

    • 最终得出正交矩阵为Q=(p1pi,,pn)Q=(p_{1},p_{i},\cdots,p_{n}),注意此时的维度为mnm*n,即mmnn列。

    (2) 求出上三角矩阵R

    由以上求正交矩阵的过程,反推:
    A1=p1y12A_{1}=p_{1}*||y_{1}||_{2}

    A2=p2y22+p1(p1TA2)A_{2}=p_{2}*||y_{2}||_{2}+p_{1}*(p_{1}^{T}A_{2})

    Ak=pkyk2+p1(p1TAk)++pk1(pk1TAk)A_{k}=p_{k}*||y_{k}||_{2}+p_{1}*(p_{1}^{T}A_{k})+\cdots+p_{k-1}*(p_{k-1}^{T}A_{k})

    使用rijr_{ij}代表矩阵R中的元素,令主对角线rii=yi2r_{ii}=||y_{i}||_{2},令rik=piTAk,i<=k1r_{ik}=p_{i}^{T}A_{k},i<=k-1。则矩阵
    A=(A1,Ai,,An)=(p1,pi,,pn)[r11r12r1n0r22r2n00rin00rnn]A=(A_{1},A_{i},\cdots,A_{n})=(p_{1},p_{i},\cdots,p_{n})* \begin{bmatrix} r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \\ \cdots & \cdots &\cdots & \cdots \\ 0 & 0 &\cdots & r_{nn} \\ \end{bmatrix}

    也即是 A=QRA=QR 此时的方法被称为消减QR分解。

    #原始的消减QR分解,A_m*n=Q_m*n·R_n*n
    function schmidit_QR(A::Array)
        m,n=size(A)
        Q=zeros(m,n)
        R=zeros(n,n)
        for j in 1:n
            y=A[:,j]
            if j!=1
                for i in 1:j-1
                    R[i,j]=Q[:,i]'*A[:,j]
                    y=y-R[i,j].*Q[:,i]
                end
            end
            R[j,j]=sqrt(sum(y.^2))
            Q[:,j]=y./R[j,j]
        end
        return Q,R
    end
    

    (3) 改进的消减QR分解

    改进的消减QR分解的思路在于:改进矩阵R中的rijr_{ij}

    此时计算的不再是当前向量AjA_{j}在向量pip_{i}上的投影系数(piTAj)(p_{i}^{T}A_{j})。而是,当前向量AjA_{j}减去之前所有在其他向量上的投影分量得到向量yy,即是中间每一迭代的消减向量yy,然后再对向量yy在当前正交向量上求出投影系数作为rijr_{ij}

    代码如下,仅在第10进行更改。

    #改进的消减QR分解
    function schmidit_QRc(A::Array)
        m,n=size(A)
        Q=zeros(m,n)
        R=zeros(n,n)
        for j in 1:n
            y=A[:,j]
            if j!=1
                for i in 1:j-1
                    R[i,j]=Q[:,i]'*y #注意此处发生更改
                    y=y-R[i,j].*Q[:,i]
                end
            end
            R[j,j]=sqrt(sum(y.^2))
            Q[:,j]=y./R[j,j]
        end
        return Q,R
    end
    

    测试
    对比MATLAB中的消减qr函数qr(A,0)。

    A=[1 -4; 2 3; 2 2]
    Qc,Rc=schmidit_QRc(A)
    

    在这里插入图片描述

    2. 完全QR分解

    以上的消减QR分解得出的是,Amn=QmnRnnA_{m*n}=Q_{m*n}*R_{n*n},当n<mn<m时,矩阵QmnQ_{m*n}不是方阵。

    完全QR分解即是将以上矩阵分解为Amn=QmmRmnA_{m*n}=Q_{m*m}*R_{m*n}的形式,以更好的利用正交方阵的性质QmmT=Qmm1Q_{m*m}^{T}=Q_{m*m}^{-1},即单位正交矩阵的转置即是其逆矩阵。

    • 基于以下事实:
      (A1,Ai,,An)=(p1,pi,,pn,,pm)[r11r12r1n0r22r2n00rin00rnn000000](A_{1},A_{i},\cdots,A_{n})=(p_{1},p_{i},\cdots,p_{n},\cdots,p_{m})* \begin{bmatrix} r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \\ \cdots & \cdots &\cdots & \cdots \\ 0 & 0 &\cdots & r_{nn} \\ 0 & 0 &\cdots & 0 \\ 0 & 0 &\cdots & 0 \\ \end{bmatrix} 此时的RR矩阵最后具有mnm-n行全0向量。
      根据以上的表达式,关键点是,在矩阵QQ中如何求出额外的mnm-n个正交的单位列向量,而前n个列向量同消解QR分解得出的列向量一致。

    • 解决思路是,在输入的矩阵(A1,Ai,,An)(A_{1},A_{i},\cdots,A_{n})后在补mnm-n个与之前向量都线性无关的列向量,把其补全为mmm*m的方阵。

    • 如何确定向量是线性无关的。对于方针而言,求方阵的行列式,如果行列式不为零,则所有列向量线性无关。当然确定矩阵线性无关的方法,可以使用初等变换,或者利用最小二乘,在此使用行列式。

    #矩阵的行列式的计算,利用递归和划分矩阵方法。
    function mydet(A)
        m,n=size(A)
        if m!=n
            throw("Error: the input Array must be a square matrix")
        else
            if m==2
                return A[1,1]*A[2,2]-A[1,2]*A[2,1]
            else
                sum=0
                for i in 1:m
                    #对矩阵进行按列拆分
                    del_row=A[setdiff(1:end,i),:] #集合操作,删除第i行 #经典
                    cut_mat=del_row[:,setdiff(1:end,1)] #删除第1列
                    sum += (-1)^((1+i)%2)*A[i,1]*mydet(cut_mat)
                end
                return sum
            end
        end
    end
    
    
    #将矩阵补全为方阵,并利用行列式 判断列向量之间是否线性无关
    function repair_mat(A::Array)
        m,n=size(A)
        if m==n
            return A
        else
            res=zeros(m,m)
            while true
                rA=[A randn(m,abs(m-n))] #矩阵拼接
                if abs(mydet(rA))>eps(Float64)
                    res=rA
                    break
                end
            end
            return res
        end
    end
    
    #改进的完全QR分解,此时 A_m*n=Q_m*m·R_m*n
    function schmidit_QRfull(B::Array)
        A=repair_mat(B)
        m,n=size(A)
        Q=zeros(m,n)
        R=zeros(n,n)
        for j in 1:n
            y=A[:,j]
            if j!=1
                for i in 1:j-1
                    R[i,j]=Q[:,i]'*y
                    y=y-R[i,j].*Q[:,i]
                end
            end
            R[j,j]=sqrt(sum(y.^2))
            Q[:,j]=y./R[j,j]
        end
        mb,nb=size(B)
        return Q,R[:,1:nb] #此时的R矩阵只需要前nb列即可
    end
    

    验证
    MATLAB中使用qr(A)函数进行完全QR分解。

    A=[1 -4; 2 3; 2 2]
    Qf,Rf=schmidit_QRfull(A)
    

    在这里插入图片描述

    3. 矩阵QR分解的作用

    (1) 方便求解线性方程组Ax=bAx=b

    完全QR分解后:QRx=bQRx=b -->> Rx=QTbRx=Q^{T}b
    由于矩阵RR是上三角矩阵,此时可以回代得出所有的xix_{i}

    (2) 求解最小二乘 minAxb2min ||Ax-b||_{2}

    由于对于单位正交矩阵存在以下性质:
    Qx2=x2||Qx||_{2}=||x||_{2},向量乘以单位正交矩阵模长不变。
    Qx2=(Qx)T(Qx)=xTQTQx=xTx=x2 ||Qx||_{2}=(Qx)^{T}(Qx)=x^{T}Q^{T}Qx=x^{T}x=||x||_{2}

    此时,
    minAxb2=minQRxb2=minRxQTb2min ||Ax-b||_{2}=min||QRx-b||_{2}=min||Rx-Q^{T}b||_{2}

    注意到:
    RxQTb=[r11r12r1n0r22r2n00rin00rnn000000][x1xn]QT[b1bnbm]Rx-Q^{T}b= \begin{bmatrix}r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \\ \cdots & \cdots &\cdots & \cdots \\ 0 & 0 &\cdots & r_{nn} \\ 0 & 0 &\cdots & 0 \\ 0 & 0 &\cdots & 0 \\ \end{bmatrix} \begin{bmatrix}x_{1}\\ \cdots \\ x_{n}\end{bmatrix}-Q^{T} \begin{bmatrix}b_{1}\\ \cdots \\ b_{n}\\ \cdots\\ b_{m}\end{bmatrix}

    由于RxQTbRx-Q^{T}b的前nn项与自变量xix_{i}有关,而后mnm-n项为确定性的常数项,与xix_{i}无关。
    令:
    [d1dm]=QT[b1bm]\begin{bmatrix}d_{1}\\ \vdots\\ d_{m}\end{bmatrix}=Q^{T} \begin{bmatrix}b_{1}\\ \vdots \\ b_{m}\end{bmatrix}

    此时二范数最小化问题可以转化为前n项每一项都最小化的问题,也即是:
    [r11r12r1n0r22r2n00rin][x1xn][d1dn]=0\begin{bmatrix}r_{11}& r_{12}&\cdots & r_{1n}\\ 0 & r_{22}&\cdots & r_{2n} \\ 0 & 0 &\cdots & r_{in} \end{bmatrix} \begin{bmatrix}x_{1}\\ \vdots \\ x_{n}\end{bmatrix}-\begin{bmatrix}d_{1}\\ \vdots \\ d_{n}\end{bmatrix}=0

    回代以上线性方程组,可以使得did_{i}的前nn项为0,此时的最小二乘误差为:
    e22=dn+12++dm2||e||_{2}^{2}=d_{n+1}^{2}+\cdots+d_{m}^{2}

    (3) 求解特征值

    矩阵特征值的数值求解方法中,一般有求部分特征值和特征向量的幂法和反幂法,以及求取全部特征值的QR分解方法。在此关注QR分解方法。

    对于矩阵特征值的求解,QR分解中利用相似矩阵的性质。

    • 如果矩阵AAA0A_{0}相似(即存在可逆矩阵PP,使得A=PA0P1A=PA_{0}P^{-1},则矩阵AAA0A_{0}相似),则相似矩阵具有相同的特征值。

    非奇异矩阵特征值的性质。

    • 对于nn阶的实数非奇异矩阵,存在nn个特征值。但是特征值中可能存在重根特征值,或者成对的复数特征值。

    (3.1) 平移QR算法求解矩阵特征值

    定理:对于实数元素的方阵AA,存在正交矩阵QQ和实数舒尔形式的矩阵TT,满足A=QTTQA=Q^{T}TQ。其中,实数舒尔矩阵TT是主对角线存在2 x 2的矩阵块的矩阵。例如:
    [xxxxxxxxxxxxxxxxx]\begin{bmatrix} x & x & x & x & x \\ x & x & x & x & x\\ & & x & x & x\\ & & x & x & x\\ & & & & x \end{bmatrix}

    这种形式的矩阵特征值是其对角块矩阵的特征值,即是对应的对角线元素(其为1x1的块)或者2x2块矩阵的特征值。矩阵AA的舒尔分解也被称为“特征值发现分解”,如果能够实现这个分解,则能够找到相应的特征值[1]。

    完整的QR算法通过一系列相似变换迭代,将任意矩阵移动到它的舒尔分解形式。

    平移QR算法原理:

    • 对于非奇异方阵AA,定义A0=AA_{0}=A
      A0sI=Q1R1A_{0}-sI=Q_{1}R_{1} A1=R1Q1+sIA_{1}=R_{1}Q_{1}+sI
      其中II为单位矩阵,ss为标量系数。由于
      A1sI=R1Q1=Q1T(A0sI)Q1=Q1TA0Q1sIA_{1}-sI=R_{1}Q_{1}=Q_{1}^{T}(A_{0}-sI)Q_{1}=Q_{1}^{T}A_{0}Q_{1}-sI因此A0A1A_{0}和A_{1}是相似矩阵,因而,AkA1A_{k}和A_{1}也是相似矩阵。

    步骤

    • 在每一步中使用QR步骤,将任意矩阵移动到它的实数舒尔形式,然后检查最后一行。如果最后一行除了对角线元素anna_{nn}之外所有元素都很小,则这个值就是特征值。在余下的计算中去掉最后一行和最后一列进行收缩,然后在迭代寻找其他特征值 [1]。
    #eyes(n)产生n阶单位矩阵,eyes(n,m)产生对角线为1的矩阵
    function eyes(m::Int,index::Int=0)
        n=0
        if index==0
            n=m
        else
            n=index
        end
        I=zeros(m,n)
        [I[i,i]=1 for i in 1:min(m,n)]
        return I
    end
    
    # 平移QR算法求特征值,原理:逆向幂迭代和收缩技术推出平移QR
    # 仅计算方阵的实数特征值,并且A中没有大小相同的特征值
    # lam为返回的特征值
    function shiftQR(A::Array)
        tol=1e-14 #精度门限
        m=size(A,1)
        lam=zeros(m,1)
        n=m
        while n>1
            #abs()函数针对标量,参数为矢量时需要在函数后面加点号
            while sort(abs.(A[n,1:n-1]))[end]>tol #sort(abs(A[n,1:n-1]))[end] 与max(abs(A[n,1:n-1]))等价
                mu=A[n,n]
                q,r=schmidit_QRfull(A-mu*eyes(n))
                A=r*q+mu*eyes(n)
            end
            lam[n]=A[n,n]
            n=n-1
            A=A[1:n,1:n]
        end
        lam[1]=A[1,1] #存在错误
        return lam
    end
    
    #平移QR算法,能够计算方阵的实数和复数特征值
    function shiftQRc(A::Array) #A为方阵
        tol=1e-14
        maxiter=500
        m=size(A,1)
        lam=zeros(m,1)
        n=m
        while n>1
            iter=0
            while sort(abs.(A[n,1:n-1]))[end]>tol && iter<maxiter
                iter=iter+1 #记录qr的次数
                mu=A[n,n] #平移量
                q,r=schmidit_QRfull(A-mu*eyes(n))
                A=r*q+mu*eyes(n)
            end
            if iter<maxiter #具有1x1的块
                lam[n]=A[n,n]
                n=n-1
                A=A[1:n,1:n]
            else#迭代超过最大迭代次数,认为是2x2的块
                #lam为实数向量,需要转化为复数向量,否则后面报错,复数不能直接复制给实数变量
                lam=complex(lam)
                disc=(A[n-1,n-1]-A[n,n])^2+4*A[n,n-1]*A[n-1,n] #一元二次求根公式 a^2-4ac
                temp=sqrt.(disc+0*im)#需要添加0*im变为复数,才能开方
                lam[n]=(A[n-1,n-1]+A[n,n] +temp)/2
                lam[n-1]=(A[n-1,n-1]+A[n,n]-temp)/2
                n=n-2 #以2收缩
                A=A[1:n,1:n]
            end
        end
        if n>0
            lam[1]=A[1,1]
        end
        return lam
    end        
    

    测试
    在这里插入图片描述

    (3.2) 以上平移QR分解存在的问题

    以上的平移QR分解算法能够求解大多数方阵的特征值,但是对于具有重数的复数特征值,以上的算法无法将原始矩阵移动为舒尔形式。例如对于矩阵A=[0 0 0 1;0 0 -1 0;0 1 0 0;-1 0 0 0],在MATLAB环境下,求得的特征值为[0+1im;01im;0+1im;01im][0+1*im;0-1*im;0+1*im;0-1*im],而以上方法求出的特征值全为0。此时需要,首先对原始矩阵进行预处理,然后再使用以上算法求解。

    (3.3) 解决方案

    上海森伯格矩阵
    在矩阵中,如果元素aij=0i>j+1a_{ij}=0,i>j+1,则mnm*n矩阵是上海森伯格形式。如下示例:
    [xxxxxxxxxxxxx]\begin{bmatrix} x & x & x & x \\ x & x & x & x \\ & x & x & x \\ & & x & x \\ \end{bmatrix}

    需要注意的是,问题关键在于如何将原始矩阵转化为上海森伯格形式。
    在此,使用矩阵QR分解中的House Holder反射子,将原始矩阵的每一列转化为只有第一个元素非零,其余都为零的列向量。由于我们想要的是实数舒尔形式,例如矩阵第一列,此时需要保留第一列的前两个元素,这时初始选择的向量,应当截取第一列的第二个元素到最后作为最初映射的原始向量,例如A[2:end,1]A[2:end,1],此时第一个元素a11a_{11}不参与运算。

    投影矩阵和House Holder反射子
    基本思想:对于任意向量xx,我们期望把它转变成x轴上的向量,例如w=(x2,0,,0)Tw=(||x||_{2},0,\cdots,0)^{T},即仅在第一个坐标位置非零,并且模长与向量xx相等。

    假设存在矩阵HH,使得Hx=wHx=w,向量ww如上所示,此时的矩阵HH为转换矩阵。

    考虑等腰三角形,向量x,wx,w可以看做是三角形的腰长,根据等腰三角形性质,其底边向量与垂线向量相互垂直,此时的底边向量可以表示为wxw-x,垂线向量表示为w+xw+x,二者为正交向量。
    在这里插入图片描述
    投影矩阵
    对于任意非零向量,定义矩阵PP如下所示:
    P=vvTvTvP=\dfrac{vv^{T}}{v^{T}v}
    此时的矩阵PP是关于向量vv的投影矩阵。

    矩阵PP具有如下性质:
    P2=vvTvvT(vTv)(vTv)=vvTvTv=PP^{2}=\dfrac{vv^{T}vv^{T}}{(v^{T}v)(v^{T}v)}=\dfrac{vv^{T}}{v^{T}v}=P
    对于向量三角形,向量xx在底边v=wxv=w-x的投影为PxPx。根据几何关系,向量xx减去其2倍的投影得到的向量就是向量ww,即
    x2Px=(I2P)x=wx-2Px=(I-2P)x=w

    因此,可以得出变换矩阵H=I2PH=I-2P,此时的HH被称为House-Holder反射子,具有以下性质。

    由于矩阵HH是对称正交矩阵,
    HTH=HH=(I2P)(I2P)=I4P+4P2=IH^{T}H=HH=(I-2P)(I-2P)=I-4P+4P^{2}=I此时可以利用House-Holder反射子进行矩阵的QR分解。

    根据文献[1]478页,以5阶方阵A为例。矩阵AA左乘House-Holder算子H1H_{1}H1AH_{1}A,其第一列的后三个元素全为0,此时的H1AH_{1}AAA不再相似,需要在其右侧右乘矩阵H11=H1TH_{1}^{-1}=H_{1}^{T}保持相似性,即H1AH1TH_{1}AH_{1}^{T}。依次迭代(3=5阶-2)次,即是
    H3H2H1AH1TH2TH3T=(H3H2H1)A(H3H2H1)TH_{3}H_{2}H_{1}AH_{1}^{T} H_{2}^{T}H_{3}^{T}=(H_{3}H_{2}H_{1})A(H_{3} H_{2}H_{1})^{T} 以上形式的迭代矩阵同原始矩阵AA相似。
    注意,一般对于一个nn阶方阵AA,需要(n2)(n-2)个House-Holder算子将其转化为上海森伯格形式。

    实现代码

    #将矩阵转化为上海森伯格形式
    function hessen(A::Array)
        m,n=size(A)
        v=zeros(m,m)
        for k in 1:m-2
            x=A[k+1:m,k]
            v[1:m-k,k]=-1*sign(x[1]+eps())*sqrt(sum(x.^2))*eyes(m-k,1)-x
            v[1:m-k,k]=v[1:m-k,k]./sqrt(sum(v[1:m-k,k].^2)) #归一化
            A[k+1:m,k:m]=A[k+1:m,k:m]-2*v[1:m-k,k]*v[1:m-k,k]'*A[k+1:m,k:m] #投影求出正交向量
            A[1:m,k+1:m]=A[1:m,k+1:m]-2*A[:,k+1:m]*v[1:m-k,k]*v[1:m-k,k]' #在对剩余的矩阵右乘House-Holder矩阵,保持矩阵相似性
        end
        return A
    end
    

    感兴趣的同学,建议读一下参考文献,具有详细步骤示意图,更加直观。

    测试
    在这里插入图片描述
    需要注意的是,对于编写的hessen()函数在MATLAB环境下对测试矩阵A进行运算结果正确。然而,在Julia环境下,对于整数矩阵直接求解运算,报错,但是在加上精度数据eps()之后,再运算结果正确,目前还没有搞明白原因。注意矩阵和标量的运算,需要使用点号进行广播操作。

    总结

    总体而言,数值方法求解矩阵特征值,使用的都是迭代思想和矩阵的相似性质。House-Holder投影算子,启发于等腰三角形,十分巧妙。
    写到现在,使用QR分解,只求解出了矩阵的全部特征值,距离相应的特征向量仅有一步之遥。累了,后续休息休息,再进行更新和完善。
    参考文献
    [1] 裴玉茹,马庚宇译,数值分析(第二版),机械工业出版社。

    展开全文
  • 数值分析--矩阵QR分解的三种方法

    万次阅读 2015-08-27 10:33:26
    QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交相似变化成为Hessenberg矩阵,然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为...
  • 矩阵分析之QR分解

    万次阅读 2016-10-24 16:38:57
    1.定义QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交相似变化成为Hessenberg矩阵,然后再应用QR方法求特征值和特征向量,在分解过程中利用Gram−SchmidtGram-Schmidt正交化。...
  • QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交相似变化成为Hessenberg矩阵,然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规正交矩阵Q与上三角形矩阵R,所以称为...
  • QR方法

    千次阅读 2018-05-12 20:40:49
    它可以一般矩阵全部特征值和特征向量。为了使用这种方法我们首先要知道一些理论基础。 Householder变换 (课本上的话虽然严谨,但是枯燥。总是感觉大学课本上的东西是给会的人看的,这里我用自己语言聊一下...
  • 矩阵分解

    千次阅读 2015-03-04 17:05:51
    下面谈谈矩阵分解的应用 ... QR分解法是目前求一般矩阵全部特征值的最有效并广泛应用的方法,一般矩阵先经过正交相似变化成为Hessenberg矩阵,然后再应用QR方法求特征值和特征向量。它是将矩阵分解成一个正规
  • 文档 作业六 6-1 试验目的计算特征值实现算法 试验容随机产生一个 10 阶整数矩阵各数均在 -5 和 5 之间 1 用 MATLAB 函数 eig 求矩阵全部特征值 2 用幂法求 A 的主特征值及对应的特征向量 3 用基本 QR 算法求全部...
  • 求矩阵秩的全选主元高斯消去法、对称正定矩阵的乔里斯基分解与行列式的求、矩阵的三角分解 、一般实矩阵的QR分解、一般实矩阵的奇异分解 、求广义逆的奇异分解法、约化对称矩阵为对称三对角阵的豪斯荷尔德变换...
  • 利用带双步位移的QR分解法,求矩阵全部特征值和特征向量,利用到拟三角化,反幂法,在VC中调试通过
  • CMatrix

    热门讨论 2008-07-21 12:30:31
    约当法 对称正定矩阵的求逆 托伯利兹矩阵求逆的埃兰特方法 求行列式的全选主元高斯消去法 求矩阵秩的全选主元高斯消去法 对称正定矩阵的乔里斯基分解与行列式的求 矩阵的三角分解 一般实...
  • 5.4 赫申伯格矩阵全部特征值QR方法 5.5 求实对称矩阵特征值与特征向量的雅可比法 5.6 求实对称矩阵特征值与特征向量的雅可比过关法 第6章 线性代数方程组的求解 6.1 求解实系数方程组的全选主元高斯消去法 6.2 ...
  • 矩阵特征值与特征向量的计算75 2.1 对称三对角阵的全部特征值与特征向量75 2.2 求实对称矩阵全部特征值与特征向量的 豪斯荷尔德变换法80 2.3 赫申伯格矩阵全部特征值QR方法88 2.4 一般实矩阵的全部特征值95 ...
  • C++常用算法程序集

    热门讨论 2011-03-24 09:15:32
    2.3 赫申伯格矩阵全部特征值QR方法88 2.4 一般实矩阵的全部特征值95 2.5 求实对称矩阵特征值与特征向量的雅可比法102 2.6 求实对称矩阵特征值与特征向量的雅可比过关法109 第3章 线性代数方程组的求解115 3.1 ...
  • 2.3 赫申伯格矩阵全部特征值QR方法88 2.4 一般实矩阵的全部特征值95 2.5 求实对称矩阵特征值与特征向量的雅可比法102 2.6 求实对称矩阵特征值与特征向量的雅可比过关法109 第3章 线性代数方程组的求解115 3.1 ...
  • 3.17 赫申伯格矩阵全部特征值QR方法 3.18 求实对称矩阵特征值与特征向量的雅可比法 3.19 求实对称矩阵特征值与特征向量的雅可比过关法 第4章 线性代数方程组的求解 4.1 线性方程组类设计 4.2 全选主元高斯消去法...
  • 5.1 约化对称矩阵为对称三对角阵的豪斯荷尔德变换法 5.2 对称三对角阵的全部特征值与特征向量 5.3 约化一般实矩阵为赫申伯格矩阵的初等相似变换法 5.4 赫身伯格矩阵全部特征QR方法 5.5 求实对称矩阵特征值与...
  • 5.4 赫身伯格矩阵全部特征QR方法 5.5 求实对称矩阵特征值与特征向量的雅可比法 5.6 求实对称矩阵特征值与特征向量的雅可比过关法 第6章 线性代数方程组的求解 6.1 求解实系数方程组的全选主元高斯消去法 6.2 ...
  • 5.4 赫身伯格矩阵全部特征QR方法 5.5 求实对称矩阵特征值与特征向量的雅可比法 5.6 求实对称矩阵特征值与特征向量的雅可比过关法 第6章 线性代数方程组的求解 6.1 求解实系数方程组的全选主元高斯消去法 6.2 ...
  • 科学与工程数值计算算法

    热门讨论 2010-12-21 12:45:03
    2.17 赫申伯格矩阵全部特征值QR方法 2.18 求实对称矩阵特征值与特征向量的雅可比法 2.19 求实对称矩阵特征值与特征向量的雅可比过关法 第3章 线性代数方程组的求解 3.1 线性方程组类设计 3.2 全选主元...
  • 5.4 赫身伯格矩阵全部特征QR方法 5.5 求实对称矩阵特征值与特征向量的雅可比法 5.6 求实对称矩阵特征值与特征向量的雅可比过关法 第6章 线性代数方程组的求解 6.1 求解实系数方程组的全选主元高斯消去法 6.2 ...
  • C#数值计算算法编程 代码

    热门讨论 2009-12-11 11:31:58
    3.17 赫申伯格矩阵全部特征值QR方法 3.18 求实对称矩阵特征值与特征向量的雅可比法 3.19 求实对称矩阵特征值与特征向量的雅可比过关法 第4章 线性代数方程组的求解 4.1 线性方程组类设计 4.2 全选主元高斯消去法...
  • 5.4 赫身伯格矩阵全部特征QR方法 5.5 求实对称矩阵特征值与特征向量的雅可比法 5.6 求实对称矩阵特征值与特征向量的雅可比过关法 第6章 线性代数方程组的求解 6.1 求解实系数方程组的全选主元高斯消去法 6.2 ...

空空如也

空空如也

1 2
收藏数 37
精华内容 14
关键字:

qr方法求矩阵全部特征值