QR方法求解特征值、特征向量(附Matlab、C代码)

QR方法特征分解

• • 1. 实对称矩阵的QR方法特征分解
  • 2. 代码实现与验证
  • 3. 参考资料

    1. 实对称矩阵的QR方法特征分解

    ​ 使用QR分解的方法求取实对称矩阵特征值、特征向量，首先通过householder变换将实对称矩阵转化为三对角矩阵，对三对角矩阵进行QR分解迭代

    [ × × × × × × × × × × × × × × × × ] → H = I − 2 w ∗ w T [ × × 0 0 × × × 0 0 × × × 0 0 × × ] \begin{bmatrix}\times&\times&\times&\times\\\times&\times&\times&\times\\\times&\times&\times&\times\\\times&\times&\times&\times\end{bmatrix}\xrightarrow{H=I-2w*w^T}\begin{bmatrix}\times&\times&0&0\\\times&\times&\times&0\\0&\times&\times&\times\\0&0&\times&\times\end{bmatrix} ​××××​××××​××××​××××​ ​H=I−2w∗wT ​ ​××00​×××0​0×××​00××​ ​

    

    针对Hessenberg矩阵多采用带位移的QR算法1(如: Rayleigh quotient shift，Wilkinson shift)

    { A k − σ ∗ I = Q k ∗ R k A k + 1 = R k ∗ Q k + σ ∗ I \left\{ \begin{array}{l} {A}_{k} - \sigma * I = {Q}_{k} * {R}_{k} \\ {A}_{k + 1} = {R}_{k} * {Q}_{k} + \sigma * I \end{array}\right. {Ak​−σ∗I=Qk​∗Rk​Ak+1​=Rk​∗Qk​+σ∗I​

    记 Q ^ k = Q 0 ∗ Q 1 ∗ … ∗ Q k {\widehat{Q}}_{k} = {Q}_{0} * {Q}_{1} * \ldots * {Q}_{k} Q ​k​=Q0​∗Q1​∗…∗Qk​ Rayleigh商加速取 σ = H n n \sigma = {H}_{nn} σ=Hnn​ ,(如下 H H H 是QR迭代中的某一步 A k {A}_{k} Ak​ )

    H = [ × × × × × × × × 0 × × × 0 0 α λ n ] \begin{array}{r} H = \left\lbrack \begin{matrix} \times & \times & \times & \times \\ \times & \times & \times & \times \\ 0 & \times & \times & \times \\ 0 & 0 & \alpha & {\lambda }_{n} \end{matrix}\right\rbrack \end{array} H= ​××00​×××0​×××α​×××λn​​ ​​

    由特征多项式知道,当 α = 0 \alpha = 0 α=0 时, λ n {\lambda }_{n} λn​ 是 H H H 的一个特征值,实际计算时,当 α 1 iter=iter+1; sigma=A(k,k); [Q,R]=QR_givens(A-sigma*eye(k));%Givens变换实现QR分解 A=R*Q+sigma*eye(k);%移位 D(1:k,k+1:end) = Q.'*D(1:k,k+1:end); hatQ(:,1:k) = hatQ(:,1:k)*Q;%记录变换矩阵 if(abs(A(k,k-1)) //对三对角矩阵进行Givens-QR分解 int i,j,k; for (i = 0; i