6.6 Hessenberg法求特征值

admin 管理员组

文章数量: 887006

6.6 Hessenberg法求特征值

文章目录

• 1 Gram-Schmidt正交化的缺点
• 2 Hessenberg矩阵
• 3 海森堡化简(Hessenberg reduction)
• 4 Givens rotation
• 5 多次Givens rotation（QR）
• 6 循环QR直至收敛
• 7 Python实现
• 8 代码测试

1 Gram-Schmidt正交化的缺点

  Gram-Schmidt正交化算法，有以下两个缺点：
  1、每一步正交分解的算法复杂度都是 O ( n 3 ) O(n^3) O(n3)，而且收敛到舒尔分解需要QR分解几次，乐观的情况下都需要n次。那么总体复杂度就是 O ( n 4 ) O(n^4) O(n4)。
  2、实践证明，当几个特征值特别接近时，收敛得特别慢，也就是说需要很多步骤才能收敛到舒尔分解。这种情况下算法复杂度就很恐怖了。

2 Hessenberg矩阵

  Hessenberg矩阵得名于Karl Hessenberg，柏林大学博士。

  Hessenberg矩阵译为海森堡矩阵，是下对角线以下为0的矩阵。举个例子：
( 1 1 1 1 1 1 1 1 1 1 0 1 1 1 1 0 0 1 1 1 0 0 0 1 1 ) \begin{pmatrix} 1 & 1& 1& 1&1\\ 1 & 1& 1& 1&1\\ 0 & 1& 1& 1&1\\ 0 & 0& 1& 1&1\\ 0 & 0& 0& 1&1\\ \end{pmatrix} ​11000​11100​11110​11111​11111​ ​
  Hessenberg法求特征值分两步，第一步是通过有限次Household reflector变换得到海森堡阵，但是我们知道海森堡阵不是拟上三角矩阵，不能用于求特征值。但是海森堡阵是原矩阵的相似矩阵。然后用一个双层循环求出相似的拟上三角矩阵，外层循环是QR分解，内层循环是Givens rotation。等于是说多次QR分解才能得到相似的拟上三角矩阵，但是每次QR分解都要执行很多次Givens旋转。得到拟上三角矩阵后，结果就出来了。

3 海森堡化简(Hessenberg reduction)

  海森堡化简的方法是Householder reflector。怎么个化简法呢？我们知道Householder变换最重要的是选择 u u u向量，u向量由以下公式给出，公式中的sign函数是取符号函数，也就是正数返回1，负数返回-1：
x = ( x 1 x 2 x 3 ⋮ x n ) ρ = − s i g n ( x 1 ) z = ( x 1 − ρ ∥ x ∥ x 2 x 3 ⋮ x n ) u = z ∥ z ∥ x=\begin{pmatrix}x_1\\x_2\\x_3\\\vdots\\x_n\\\end{pmatrix}\\ \rho=−sign(x_1)\\ z=\begin{pmatrix}x_1-\rho \| x \| \\x_2\\x_3\\\vdots\\x_n\\\end{pmatrix}\\ u=\frac z {\|z\|} x= ​x1​x2​x3​⋮xn​​ ​ρ=−sign(x1​)z= ​x1​−ρ∥x∥x2​x3​⋮xn​​ ​u=∥z∥z​  选好后，按以下步骤，先计算 P 1 P_1 P1​，按如下方法：
P 1 = ( 1 0 0 I − 2 u 1 u 1 T ) A 1 = P 1 A P 1 ⋮ A n − 2 = P n − 2 A n − 1 P n − 2 P_1=\begin{pmatrix}1&0\\0&I-2u_1u_1^T\end{pmatrix}\\ A_1=P_1AP_1\\ \vdots\\ A_{n-2}=P_{n-2}A_{n-1}P_{n-2} P1​=(10​0I−2u1​u1T​​)A1​=P1​AP1​⋮An−2​=Pn−2​An−1​Pn−2​
  P的通项公式如下：
P i = ( I 0 0 I − 2 u i u i T ) P_i = \begin{pmatrix} I&0\\ 0 & I-2u_iu_i^T \end{pmatrix} Pi​=(I0​0I−2ui​uiT​​)
  需要注意的是u向量的选择问题。每次选择，u向量的长度都越来越短。假如是矩阵的第一列，u向量来源于这个向量x：
x 1 = ( A 2 , 1 A 3 , 1 A 4 , 1 ⋮ A n , 1 ) x_1=\begin{pmatrix} A_{2,1}\\ A_{3,1}\\ A_{4,1}\\ \vdots\\ A_{n,1} \end{pmatrix} x1​= ​A2,1​A3,1​A4,1​⋮An,1​​ ​
  而 x x x向量的通项公式如下：
x i = ( A i + 1 , 1 A i + 2 , 1 A i + 3 , 1 ⋮ A n , 1 ) x_i=\begin{pmatrix} A_{i+1,1}\\ A_{i+2,1}\\ A_{i+3,1}\\ \vdots\\ A_{n,1} \end{pmatrix} xi​= ​Ai+1,1​Ai+2,1​Ai+3,1​⋮An,1​​ ​
   A n − 2 A_{n-2} An−2​就是和A相似的海森堡矩阵。海森堡矩阵用Gram-Schmidt法QR分解后，RQ的结果还是个海森堡矩阵，永远结束不了QR分解循环。
  我举个例子：
A 0 = ( 2.0 1.0 − 1.0 11.0 16.0 1.0 2.0 − 1.0 3.0 17.0 − 1.0 − 1.0 2.0 4.0 − 4.0 7.0 10.0 9.0 5.0 − 5.0 8.0 11.0 6.0 12.0 − 6.0 ) A 1 = ( 2.0 − 19.303 0.732 − 1.122 2.146 − 10.724 4.061 − 11.698 − 0.734 18.797 0.0 − 0.966 2.895 4.444 − 4.01 0.0 11.706 2.572 3.055 − 3.602 0.0 9.129 − 1.02 7.495 − 7.01 ) A 2 = ( 2.0 − 19.303 0.386 − 0.867 2.345 − 10.724 4.061 11.717 − 18.036 5.305 0.0 14.876 0.986 6.88 − 6.747 0.0 0.0 − 0.076 4.253 0.758 0.0 0.0 1.583 4.98 − 6.3 ) A 3 = ( 2.0 − 19.303 0.386 2.384 − 0.754 − 10.724 4.061 11.717 6.163 − 17.761 0.0 14.876 0.986 − 7.069 6.549 0.0 0.0 1.585 − 6.55 4.462 0.0 0.0 0.0 0.24 4.503 ) A_0=\begin{pmatrix}2.0 & 1.0 & -1.0 & 11.0 & 16.0\\ 1.0 & 2.0 & -1.0 & 3.0 & 17.0\\ -1.0 & -1.0 & 2.0 & 4.0 & -4.0\\ 7.0 & 10.0 & 9.0 & 5.0 & -5.0\\ 8.0 & 11.0 & 6.0 & 12.0 & -6.0\\ \end{pmatrix}\\ A_1= \begin{pmatrix}2.0 & -19.303 & 0.732 & -1.122 & 2.146\\ -10.724 & 4.061 & -11.698 & -0.734 & 18.797\\ 0.0 & -0.966 & 2.895 & 4.444 & -4.01\\ 0.0 & 11.706 & 2.572 & 3.055 & -3.602\\ 0.0 & 9.129 & -1.02 & 7.495 & -7.01\\ \end{pmatrix}\\ A_2= \begin{pmatrix}2.0 & -19.303 & 0.386 & -0.867 & 2.345\\ -10.724 & 4.061 & 11.717 & -18.036 & 5.305\\ 0.0 & 14.876 & 0.986 & 6.88 & -6.747\\ 0.0 & 0.0 & -0.076 & 4.253 & 0.758\\ 0.0 & 0.0 & 1.583 & 4.98 & -6.3\\ \end{pmatrix}\\ A_3= \begin{pmatrix}2.0 & -19.303 & 0.386 & 2.384 & -0.754\\ -10.724 & 4.061 & 11.717 & 6.163 & -17.761\\ 0.0 & 14.876 & 0.986 & -7.069 & 6.549\\ 0.0 & 0.0 & 1.585 & -6.55 & 4.462\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ A0​= ​2.01.0−1.07.08.0​1.02.0−1.010.011.0​−1.0−1.02.09.06.0​11.03.04.05.012.0​16.017.0−4.0−5.0−6.0​ ​A1​= ​2.0−10.7240.00.00.0​−19.3034.061−0.96611.7069.129​0.732−11.6982.8952.572−1.02​−1.122−0.7344.4443.0557.495​2.14618.797−4.01−3.602−7.01​ ​A2​= ​2.0−10.7240.00.00.0​−19.3034.06114.8760.00.0​0.38611.7170.986−0.0761.583​−0.867−18.0366.884.2534.98​2.3455.305−6.7470.758−6.3​ ​A3​= ​2.0−10.7240.00.00.0​−19.3034.06114.8760.00.0​0.38611.7170.9861.5850.0​2.3846.163−7.069−6.550.24​−0.754−17.7616.5494.4624.503​ ​

4 Givens rotation

  首先假设 c 2 + s 2 = 1 c^2+s^2=1 c2+s2=1，这就是一个圆的方程，所以 c = cos ⁡ ( θ ) , s = sin ⁡ ( θ ) c=\cos(\theta),s=\sin(\theta) c=cos(θ),s=sin(θ)。假设 i , j i,j i,j是矩阵A的两行， e i , e j e_i,e_j ei​,ej​是矩阵A的这两行舒尔分解后的两个特征向量。Givens rotation可以翻译为Givens旋转，是一个旋转矩阵，定义如下：
G ( i , j , c , s ) = I + ( cos ⁡ ( θ ) − 1 ) e i e i T − cos ⁡ ( θ ) e i e j T + s e j e j T + ( cos ⁡ ( θ ) − 1 ) e j e j T = ( I 0 0 0 0 0 cos ⁡ ( θ ) 0 − sin ⁡ ( θ ) 0 0 0 I 0 0 0 sin ⁡ ( θ ) 0 cos ⁡ ( θ ) 0 0 0 0 0 I ) = ( I 0 0 0 0 0 c 0 − s 0 0 0 I 0 0 0 s 0 c 0 0 0 0 0 I ) G(i,j,c,s)=I+(\cos(\theta)-1)e_ie_i^T-\cos(\theta)e_ie_j^T+se_je_j^T+(\cos(\theta)-1)e_je_j^T\\ =\begin{pmatrix} I & 0 & 0 & 0 & 0\\ 0 & \cos(\theta) & 0 & -\sin(\theta) & 0 \\ 0 & 0 & I & 0 & 0\\ 0 & \sin(\theta) & 0 & \cos(\theta) & 0\\ 0 & 0 & 0 & 0 & I \end{pmatrix}\\ =\begin{pmatrix} I & 0 & 0 & 0 & 0\\ 0 & c & 0 & -s & 0 \\ 0 & 0 & I & 0 & 0\\ 0 & s & 0 & c & 0\\ 0 & 0 & 0 & 0 & I \end{pmatrix} G(i,j,c,s)=I+(cos(θ)−1)ei​eiT​−cos(θ)ei​ejT​+sej​ejT​+(cos(θ)−1)ej​ejT​= ​I0000​0cos(θ)0sin(θ)0​00I00​0−sin(θ)0cos(θ)0​0000I​ ​= ​I0000​0c0s0​00I00​0−s0c0​0000I​ ​
  公式中的矩阵的 cos ⁡ ( θ ) \cos(\theta) cos(θ)和 − sin ⁡ ( θ ) -\sin(\theta) −sin(θ)位于第 i i i行，同理 sin ⁡ ( θ ) \sin(\theta) sin(θ)和 cos ⁡ ( θ ) \cos(\theta) cos(θ)位于第j行。Givens旋转的几何意义如下，图片来自KTH（瑞典皇家理工学院）：

  图中 e i e_i ei​和 e j e_j ej​是两个特征向量， G G G将x在 e i e_i ei​和 e j e_j ej​组成的向量空间（平面）上逆时针旋转了 θ \theta θ。

5 多次Givens rotation（QR）

  我们已经得到了一个Hessenberg矩阵，现在就是对Hessenberg矩阵进行QR分解了，然后得到和Hessenberg矩阵相似的矩阵。这个矩阵就是我们想要的求特征值矩阵。在前面说过，对Hessenberg进行Gram-Schmidt分解，得到的还是Hessenberg矩阵，那怎么办呢？
  对Hessenberg矩阵进行N次Givens旋转，假设i为Givens旋转的次数，第一次旋转 i = 1 i=1 i=1，原始Hessenberg矩阵为 H H H,将 H H H赋值给 H 1 H_1 H1​，对于Hessenberg矩阵的QR分解公式如下：
H 1 = H r i = ( H i ) i , i 2 + ( H i ) i + 1 , i 2 c = cos ⁡ ( θ ) = ( H i ) i , i r i s = sin ⁡ ( θ ) = ( H i + 1 ) i , i r i G i = G ( i , i + 1 , c , s ) H i + 1 = G i T H i H = （ ∏ i = 1 n − 1 G i ） H n = Q R Q = ∏ i = 1 n − 1 G i R = H n H ′ = R Q H_1=H\\ r_i=\sqrt{(H_i)_{i,i}^2+(H_i)_{i+1,i}^2}\\ c=\cos(\theta)=\frac{(H_i)_{i,i}}{r_i}\\ s=\sin(\theta)=\frac{(H_{i+1})_{i,i}}{r_i}\\ G_i=G(i,i+1,c,s)\\ H_{i+1}=G_i^TH_i\\ H=（\prod_{i=1}^{n-1}G_i）H_n=QR\\ Q=\prod_{i=1}^{n-1}G_i\\ R=H_n\\ H'=RQ H1​=Hri​=(Hi​)i,i2​+(Hi​)i+1,i2​ ​c=cos(θ)=ri​(Hi​)i,i​​s=sin(θ)=ri​(Hi+1​)i,i​​Gi​=G(i,i+1,c,s)Hi+1​=GiT​Hi​H=（i=1∏n−1​Gi​）Hn​=QRQ=i=1∏n−1​Gi​R=Hn​H′=RQ
  得到的 H ′ H' H′与原先的Hessenberg矩阵 H H H相似。但是这个RQ不一定是拟上三角矩阵。这个公式是如此地复杂，我同样建议收藏本文，不要记忆，不要记忆，不要记忆！重要的事情说三遍。
  举个例子：
A = ( 2.0 1.0 − 1.0 11.0 16.0 1.0 2.0 − 1.0 3.0 17.0 − 1.0 − 1.0 2.0 4.0 − 4.0 7.0 10.0 9.0 5.0 − 5.0 8.0 11.0 6.0 12.0 − 6.0 ) H = ( 2.0 − 19.303 0.386 2.384 − 0.754 − 10.724 4.061 11.717 6.163 − 17.761 0.0 14.876 0.986 − 7.069 6.549 0.0 0.0 1.585 − 6.55 4.462 0.0 0.0 0.0 0.24 4.503 ) H 1 = ( 2.0 − 19.303 0.386 2.384 − 0.754 − 10.724 4.061 11.717 6.163 − 17.761 0.0 14.876 0.986 − 7.069 6.549 0.0 0.0 1.585 − 6.55 4.462 0.0 0.0 0.0 0.24 4.503 ) G 1 = ( 0.183 0.983 0.0 0.0 0.0 − 0.983 0.183 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 ) H 2 = ( 10.909 − 7.531 − 11.448 − 5.621 17.322 0.0 − 18.231 2.528 3.473 − 3.997 0.0 14.876 0.986 − 7.069 6.549 0.0 0.0 1.585 − 6.55 4.462 0.0 0.0 0.0 0.24 4.503 ) G 2 = ( 1.0 0.0 0.0 0.0 0.0 0.0 − 0.775 − 0.632 0.0 0.0 0.0 0.632 − 0.775 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 ) H 3 = ( 10.909 − 7.531 − 11.448 − 5.621 17.322 0.0 23.53 − 1.335 − 7.16 7.237 0.0 0.0 − 2.362 3.281 − 2.547 0.0 0.0 1.585 − 6.55 4.462 0.0 0.0 0.0 0.24 4.503 ) G 3 = ( 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 − 0.83 − 0.557 0.0 0.0 0.0 0.557 − 0.83 0.0 0.0 0.0 0.0 0.0 1.0 ) H 4 = ( 10.909 − 7.531 − 11.448 − 5.621 17.322 0.0 23.53 − 1.335 − 7.16 7.237 0.0 0.0 2.844 − 6.374 4.601 0.0 0.0 0.0 3.611 − 2.286 0.0 0.0 0.0 0.24 4.503 ) G 4 = ( 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 0.0 0.998 − 0.066 0.0 0.0 0.0 0.066 0.998 ) Q = ( 0.183 − 0.762 0.516 0.346 − 0.023 − 0.983 − 0.142 0.096 0.064 − 0.004 0.0 0.632 0.643 0.431 − 0.029 0.0 0.0 0.557 − 0.829 0.055 0.0 0.0 0.0 0.066 0.998 ) R = ( 10.909 − 7.531 − 11.448 − 5.621 17.322 0.0 23.53 − 1.335 − 7.16 7.237 0.0 0.0 2.844 − 6.374 4.601 0.0 0.0 0.0 3.619 − 1.982 0.0 0.0 0.0 0.0 4.645 ) A= \begin{pmatrix}2.0 & 1.0 & -1.0 & 11.0 & 16.0\\ 1.0 & 2.0 & -1.0 & 3.0 & 17.0\\ -1.0 & -1.0 & 2.0 & 4.0 & -4.0\\ 7.0 & 10.0 & 9.0 & 5.0 & -5.0\\ 8.0 & 11.0 & 6.0 & 12.0 & -6.0\\ \end{pmatrix}\\ H= \begin{pmatrix}2.0 & -19.303 & 0.386 & 2.384 & -0.754\\ -10.724 & 4.061 & 11.717 & 6.163 & -17.761\\ 0.0 & 14.876 & 0.986 & -7.069 & 6.549\\ 0.0 & 0.0 & 1.585 & -6.55 & 4.462\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ H_1= \begin{pmatrix}2.0 & -19.303 & 0.386 & 2.384 & -0.754\\ -10.724 & 4.061 & 11.717 & 6.163 & -17.761\\ 0.0 & 14.876 & 0.986 & -7.069 & 6.549\\ 0.0 & 0.0 & 1.585 & -6.55 & 4.462\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ G_1= \begin{pmatrix}0.183 & 0.983 & 0.0 & 0.0 & 0.0\\ -0.983 & 0.183 & 0.0 & 0.0 & 0.0\\ 0.0 & 0.0 & 1.0 & 0.0 & 0.0\\ 0.0 & 0.0 & 0.0 & 1.0 & 0.0\\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0\\ \end{pmatrix}\\ H_2= \begin{pmatrix}10.909 & -7.531 & -11.448 & -5.621 & 17.322\\ 0.0 & -18.231 & 2.528 & 3.473 & -3.997\\ 0.0 & 14.876 & 0.986 & -7.069 & 6.549\\ 0.0 & 0.0 & 1.585 & -6.55 & 4.462\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ G_2= \begin{pmatrix}1.0 & 0.0 & 0.0 & 0.0 & 0.0\\ 0.0 & -0.775 & -0.632 & 0.0 & 0.0\\ 0.0 & 0.632 & -0.775 & 0.0 & 0.0\\ 0.0 & 0.0 & 0.0 & 1.0 & 0.0\\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0\\ \end{pmatrix}\\ H_3= \begin{pmatrix}10.909 & -7.531 & -11.448 & -5.621 & 17.322\\ 0.0 & 23.53 & -1.335 & -7.16 & 7.237\\ 0.0 & 0.0 & -2.362 & 3.281 & -2.547\\ 0.0 & 0.0 & 1.585 & -6.55 & 4.462\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ G_3= \begin{pmatrix}1.0 & 0.0 & 0.0 & 0.0 & 0.0\\ 0.0 & 1.0 & 0.0 & 0.0 & 0.0\\ 0.0 & 0.0 & -0.83 & -0.557 & 0.0\\ 0.0 & 0.0 & 0.557 & -0.83 & 0.0\\ 0.0 & 0.0 & 0.0 & 0.0 & 1.0\\ \end{pmatrix}\\ H_4= \begin{pmatrix}10.909 & -7.531 & -11.448 & -5.621 & 17.322\\ 0.0 & 23.53 & -1.335 & -7.16 & 7.237\\ 0.0 & 0.0 & 2.844 & -6.374 & 4.601\\ 0.0 & 0.0 & 0.0 & 3.611 & -2.286\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ G_4= \begin{pmatrix}1.0 & 0.0 & 0.0 & 0.0 & 0.0\\ 0.0 & 1.0 & 0.0 & 0.0 & 0.0\\ 0.0 & 0.0 & 1.0 & 0.0 & 0.0\\ 0.0 & 0.0 & 0.0 & 0.998 & -0.066\\ 0.0 & 0.0 & 0.0 & 0.066 & 0.998\\ \end{pmatrix}\\ Q= \begin{pmatrix}0.183 & -0.762 & 0.516 & 0.346 & -0.023\\ -0.983 & -0.142 & 0.096 & 0.064 & -0.004\\ 0.0 & 0.632 & 0.643 & 0.431 & -0.029\\ 0.0 & 0.0 & 0.557 & -0.829 & 0.055\\ 0.0 & 0.0 & 0.0 & 0.066 & 0.998\\ \end{pmatrix}\\ R= \begin{pmatrix}10.909 & -7.531 & -11.448 & -5.621 & 17.322\\ 0.0 & 23.53 & -1.335 & -7.16 & 7.237\\ 0.0 & 0.0 & 2.844 & -6.374 & 4.601\\ 0.0 & 0.0 & 0.0 & 3.619 & -1.982\\ 0.0 & 0.0 & 0.0 & 0.0 & 4.645\\ \end{pmatrix}\\ A= ​2.01.0−1.07.08.0​1.02.0−1.010.011.0​−1.0−1.02.09.06.0​11.03.04.05.012.0​16.017.0−4.0−5.0−6.0​ ​H= ​2.0−10.7240.00.00.0​−19.3034.06114.8760.00.0​0.38611.7170.9861.5850.0​2.3846.163−7.069−6.550.24​−0.754−17.7616.5494.4624.503​ ​H1​= ​2.0−10.7240.00.00.0​−19.3034.06114.8760.00.0​0.38611.7170.9861.5850.0​2.3846.163−7.069−6.550.24​−0.754−17.7616.5494.4624.503​ ​G1​= ​0.183−0.9830.00.00.0​0.9830.1830.00.00.0​0.00.01.00.00.0​0.00.00.01.00.0​0.00.00.00.01.0​ ​H2​= ​10.9090.00.00.00.0​−7.531−18.23114.8760.00.0​−11.4482.5280.9861.5850.0​−5.6213.473−7.069−6.550.24​17.322−3.9976.5494.4624.503​ ​G2​= ​1.00.00.00.00.0​0.0−0.7750.6320.00.0​0.0−0.632−0.7750.00.0​0.00.00.01.00.0​0.00.00.00.01.0​ ​H3​= ​10.9090.00.00.00.0​−7.53123.530.00.00.0​−11.448−1.335−2.3621.5850.0​−5.621−7.163.281−6.550.24​17.3227.237−2.5474.4624.503​ ​G3​= ​1.00.00.00.00.0​0.01.00.00.00.0​0.00.0−0.830.5570.0​0.00.0−0.557−0.830.0​0.00.00.00.01.0​ ​H4​= ​10.9090.00.00.00.0​−7.53123.530.00.00.0​−11.448−1.3352.8440.00.0​−5.621−7.16−6.3743.6110.24​17.3227.2374.601−2.2864.503​ ​G4​= ​1.00.00.00.00.0​0.01.00.00.00.0​0.00.01.00.00.0​0.00.00.00.9980.066​0.00.00.0−0.0660.998​ ​Q= ​0.183−0.9830.00.00.0​−0.762−0.1420.6320.00.0​0.5160.0960.6430.5570.0​0.3460.0640.431−0.8290.066​−0.023−0.004−0.0290.0550.998​ ​R= ​10.9090.00.00.00.0​−7.53123.530.00.00.0​−11.448−1.3352.8440.00.0​−5.621−7.16−6.3743.6190.0​17.3227.2374.601−1.9824.645​ ​

6 循环QR直至收敛

  这个与Gram-Schmidt法类似，直到得出拟上三角矩阵。我简单举个例子：
A = ( 2.0 1.0 − 1.0 11.0 16.0 1.0 2.0 − 1.0 3.0 17.0 − 1.0 − 1.0 2.0 4.0 − 4.0 7.0 10.0 9.0 5.0 − 5.0 8.0 11.0 6.0 12.0 − 6.0 ) H 0 = ( 2.0 − 19.303 0.386 2.384 − 0.754 − 10.724 4.061 11.717 6.163 − 17.761 0.0 14.876 0.986 − 7.069 6.549 0.0 0.0 1.585 − 6.55 4.462 0.0 0.0 0.0 0.24 4.503 ) H 1 = ( 9.403 − 14.477 − 5.593 4.159 17.083 − 23.131 − 4.187 − 2.584 7.354 6.764 0.0 1.798 − 1.722 6.812 4.159 0.0 0.0 2.016 − 3.13 − 1.778 0.0 0.0 0.0 0.308 4.635 ) H 2 = ( 10.861 − 22.344 − 5.631 − 2.116 0.441 − 13.984 − 4.941 − 8.713 − 2.546 − 17.565 0.0 0.377 − 7.401 − 4.76 − 5.441 0.0 0.0 1.537 2.258 2.956 0.0 0.0 0.0 0.544 4.224 ) H 3 = ( 18.604 − 7.906 − 3.359 6.916 12.417 − 16.333 − 12.507 − 9.063 8.736 7.467 0.0 0.141 − 6.55 8.487 2.976 0.0 0.0 0.272 2.547 2.505 0.0 0.0 0.0 1.26 2.907 ) H 4 = ( 17.079 − 18.036 − 3.649 1.374 4.205 − 9.643 − 10.896 − 9.513 − 15.274 − 8.377 0.0 0.064 − 6.984 − 8.721 0.517 0.0 0.0 0.129 4.318 2.017 0.0 0.0 0.0 0.645 1.482 ) H 5 = ( 22.167 − 0.627 − 1.341 9.822 6.35 − 9.025 − 15.949 − 9.876 13.429 3.209 0.0 0.025 − 6.856 8.794 − 1.832 0.0 0.0 0.077 4.497 1.558 0.0 0.0 0.0 0.177 1.141 ) H 6 = ( 20.118 − 14.053 − 2.459 4.237 4.512 − 5.659 − 13.884 − 9.845 − 16.212 − 4.749 0.0 0.011 − 6.971 − 8.614 2.173 0.0 0.0 0.052 4.657 1.402 0.0 0.0 0.0 0.042 1.08 ) H 8 = ( 22.763 3.744 − 0.233 8.522 5.552 − 4.65 − 16.523 − 10.042 14.557 3.22 0.0 0.005 − 6.913 8.69 − 2.239 0.0 0.0 0.034 4.606 1.386 0.0 0.0 0.0 0.01 1.067 ) A= \begin{pmatrix}2.0 & 1.0 & -1.0 & 11.0 & 16.0\\ 1.0 & 2.0 & -1.0 & 3.0 & 17.0\\ -1.0 & -1.0 & 2.0 & 4.0 & -4.0\\ 7.0 & 10.0 & 9.0 & 5.0 & -5.0\\ 8.0 & 11.0 & 6.0 & 12.0 & -6.0\\ \end{pmatrix}\\ H_0= \begin{pmatrix}2.0 & -19.303 & 0.386 & 2.384 & -0.754\\ -10.724 & 4.061 & 11.717 & 6.163 & -17.761\\ 0.0 & 14.876 & 0.986 & -7.069 & 6.549\\ 0.0 & 0.0 & 1.585 & -6.55 & 4.462\\ 0.0 & 0.0 & 0.0 & 0.24 & 4.503\\ \end{pmatrix}\\ H_1= \begin{pmatrix}9.403 & -14.477 & -5.593 & 4.159 & 17.083\\ -23.131 & -4.187 & -2.584 & 7.354 & 6.764\\ 0.0 & 1.798 & -1.722 & 6.812 & 4.159\\ 0.0 & 0.0 & 2.016 & -3.13 & -1.778\\ 0.0 & 0.0 & 0.0 & 0.308 & 4.635\\ \end{pmatrix}\\ H_2= \begin{pmatrix}10.861 & -22.344 & -5.631 & -2.116 & 0.441\\ -13.984 & -4.941 & -8.713 & -2.546 & -17.565\\ 0.0 & 0.377 & -7.401 & -4.76 & -5.441\\ 0.0 & 0.0 & 1.537 & 2.258 & 2.956\\ 0.0 & 0.0 & 0.0 & 0.544 & 4.224\\ \end{pmatrix}\\ H_3= \begin{pmatrix}18.604 & -7.906 & -3.359 & 6.916 & 12.417\\ -16.333 & -12.507 & -9.063 & 8.736 & 7.467\\ 0.0 & 0.141 & -6.55 & 8.487 & 2.976\\ 0.0 & 0.0 & 0.272 & 2.547 & 2.505\\ 0.0 & 0.0 & 0.0 & 1.26 & 2.907\\ \end{pmatrix}\\ H_4= \begin{pmatrix}17.079 & -18.036 & -3.649 & 1.374 & 4.205\\ -9.643 & -10.896 & -9.513 & -15.274 & -8.377\\ 0.0 & 0.064 & -6.984 & -8.721 & 0.517\\ 0.0 & 0.0 & 0.129 & 4.318 & 2.017\\ 0.0 & 0.0 & 0.0 & 0.645 & 1.482\\ \end{pmatrix}\\ H_5= \begin{pmatrix}22.167 & -0.627 & -1.341 & 9.822 & 6.35\\ -9.025 & -15.949 & -9.876 & 13.429 & 3.209\\ 0.0 & 0.025 & -6.856 & 8.794 & -1.832\\ 0.0 & 0.0 & 0.077 & 4.497 & 1.558\\ 0.0 & 0.0 & 0.0 & 0.177 & 1.141\\ \end{pmatrix}\\ H_6= \begin{pmatrix}20.118 & -14.053 & -2.459 & 4.237 & 4.512\\ -5.659 & -13.884 & -9.845 & -16.212 & -4.749\\ 0.0 & 0.011 & -6.971 & -8.614 & 2.173\\ 0.0 & 0.0 & 0.052 & 4.657 & 1.402\\ 0.0 & 0.0 & 0.0 & 0.042 & 1.08\\ \end{pmatrix}\\ H_8= \begin{pmatrix}22.763 & 3.744 & -0.233 & 8.522 & 5.552\\ -4.65 & -16.523 & -10.042 & 14.557 & 3.22\\ 0.0 & 0.005 & -6.913 & 8.69 & -2.239\\ 0.0 & 0.0 & 0.034 & 4.606 & 1.386\\ 0.0 & 0.0 & 0.0 & 0.01 & 1.067\\ \end{pmatrix}\\ A= ​2.01.0−1.07.08.0​1.02.0−1.010.011.0​−1.0−1.02.09.06.0​11.03.04.05.012.0​16.017.0−4.0−5.0−6.0​ ​H0​= ​2.0−10.7240.00.00.0​−19.3034.06114.8760.00.0​0.38611.7170.9861.5850.0​2.3846.163−7.069−6.550.24​−0.754−17.7616.5494.4624.503​ ​H1​= ​9.403−23.1310.00.00.0​−14.477−4.1871.7980.00.0​−5.593−2.584−1.7222.0160.0​4.1597.3546.812−3.130.308​17.0836.7644.159−1.7784.635​ ​H2​= ​10.861−13.9840.00.00.0​−22.344−4.9410.3770.00.0​−5.631−8.713−7.4011.5370.0​−2.116−2.546−4.762.2580.544​0.441−17.565−5.4412.9564.224​ ​H3​= ​18.604−16.3330.00.00.0​−7.906−12.5070.1410.00.0​−3.359−9.063−6.550.2720.0​6.9168.7368.4872.5471.26​12.4177.4672.9762.5052.907​ ​H4​= ​17.079−9.6430.00.00.0​−18.036−10.8960.0640.00.0​−3.649−9.513−6.9840.1290.0​1.374−15.274−8.7214.3180.645​4.205−8.3770.5172.0171.482​ ​H5​= ​22.167−9.0250.00.00.0​−0.627−15.9490.0250.00.0​−1.341−9.876−6.8560.0770.0​9.82213.4298.7944.4970.177​6.353.209−1.8321.5581.141​ ​H6​= ​20.118−5.6590.00.00.0​−14.053−13.8840.0110.00.0​−2.459−9.845−6.9710.0520.0​4.237−16.212−8.6144.6570.042​4.512−4.7492.1731.4021.08​ ​H8​= ​22.763−4.650.00.00.0​3.744−16.5230.0050.00.0​−0.233−10.042−6.9130.0340.0​8.52214.5578.694.6060.01​5.5523.22−2.2391.3861.067​ ​
  特征值为 22.315 , − 16.075 , 4.631 , − 6.938 , 1.067 22.315, -16.075, 4.631, -6.938, 1.067 22.315,−16.075,4.631,−6.938,1.067。

7 Python实现

  这么复杂的代码，不愿拿Java写，为了省时间，就python了：

# _*_ coding:utf-8 _*_import math
import cmathMIN_VALUE = 1e-2def round_num(x):if isinstance(x, complex):return complex(round_num(x.real), round_num(x.imag))return round(x * 1000) / 1000.0class Matrix:# 矩阵，为了方便计算，二维数组的每个一维数组表示一个向量，也就是一列# 所以看起来会比较奇怪def __init__(self, lines):self.__lines = linesdef __mul__(self, other):if len(self.__lines) != len(other.__lines[0]):raise Exception("矩阵A列数%d != 矩阵B的行数%d" % (len(self.__lines[0]), len(other.__lines)))# 弄一个m行n列的新矩阵m = len(self.__lines[0])n = len(other.__lines)p = len(other.__lines[0])result = [[0] * n for _ in range(0, m)]# i 代表self的行for i in range(0, m):# j 代表 B 矩阵的列for j in range(0, n):# 第一个矩阵的行 与第二个矩阵列的乘积和# k 代表 A矩阵的列和B矩阵的行for k in range(0, p):mul = self.__lines[k][i] * other.__lines[j][k]result[j][i] += mulreturn Matrix(result)def __str__(self):s = ''rows = len(self.__lines)columns = len(self.__lines[0])for i in range(0, columns):s += "|"for j in range(0, rows):if j > 0:s += '\t's += str(round(self.__lines[j][i] * 1000) / 1000)s += "\t|\n"return sdef eigen_values(self):h = self.to_hessenberg()while not h.is_upper_triangle():q, r = h.givens_rotations()h = r * qreturn h.__eigen_values()def is_upper_triangle(self):for i in range(len(self.__lines)):vector = self.__lines[i]for j in range(i + 1, len(vector)):e = vector[j]if abs(e) > MIN_VALUE:if j > i + 1:return False# 这时候再判断是不是拟三角矩阵# 如果左上角是不是0if i > 0 and abs(self.__lines[i - 1][j - 1]) > MIN_VALUE:return False# 如果右下角是不是0if i < len(self.__lines) - 1 and abs(self.__lines[i + 1][j + 1]) > MIN_VALUE:return Falsereturn Truedef __eigen_values(self):eigen_values = [0 for _ in self.__lines]for i, vector in enumerate(self.__lines):# 如果左边不是0，那么求过了if i > 0 and abs(self.__lines[i - 1][i]) > MIN_VALUE:continue# 如果是下一个元素不是0if i < len(self.__lines) - 1 and abs(vector[i + 1]) > MIN_VALUE:# 按求根公式计算a = vector[i]b = self.__lines[i + 1][i]c = vector[i + 1]d = self.__lines[i + 1][i + 1]root = a * a + d * d - 2 * a * d + 4 * b * cif root > 0:sqrt = math.sqrt(root)else:sqrt = cmath.sqrt(root)a_d = a + deigen_values[i] = round_num((a_d + sqrt) / 2)eigen_values[i + 1] = round_num((a_d - sqrt) / 2)else:eigen_values[i] = round(vector[i] * 1000) / 1000.0return eigen_valuesdef givens_rotations(self):h = selfn = len(self.__lines)q = Nonefor i in range(0, n - 1):g, h = h.givens_rotation(i)q = g if q is None else q * greturn q, hdef givens_rotation(self, i):vector = self.__lines[i]n = len(self.__lines)j = i + 1r = math.sqrt(vector[i] ** 2 + vector[j] ** 2)c = vector[i] / rs = vector[j] / rg = Matrix.identity_matrix(n)g.__lines[i][i] = cg.__lines[i][j] = sg.__lines[j][i] = -sg.__lines[j][j] = cg_t = Matrix.identity_matrix(n)g_t.__lines[i][i] = cg_t.__lines[i][j] = -sg_t.__lines[j][i] = sg_t.__lines[j][j] = ch = g_t * selfreturn g, hdef to_hessenberg(self):n = len(self.__lines)a = selffor i in range(0, n - 2):p = a.household_reflector(i)a = p * a * preturn adef household_reflector(self, i):j = i + 1column = self.__lines[i]rho = 1 if column[j] < 0 else -1z = [column[k] for k in range(j, len(column))]# 长度必须两次运算，因为第一个元素已经改变z[0] = z[0] - rho * Matrix.length(z)z_length = Matrix.length(z)u = [e / z_length for e in z]# 对u进行转置ut = [[e] for e in u]n = len(self.__lines)p_array = Matrix.identity_matrix(n)# 生成每一个向量# 对于i左边的就是矩阵Isub_matrix = Matrix([u]) * Matrix(ut)sub_p = sub_matrix.__lines# i 和 i以后的需要处理# 这里需要仔细处理# 对于i = 0, j = 1for k in range(j, n):for row in range(j, n):# k 以1开始p_array.__lines[k][row] = p_array.__lines[k][row] - 2 * sub_p[k - j][row - j]return p_array@staticmethoddef length(vector):return math.sqrt(Matrix.square(vector))@staticmethoddef square(vector):return Matrix.dot(vector, vector)@staticmethoddef dot(vector1, vector2):l = 0for i in range(0, len(vector1)):l += vector1[i] * vector2[i]return l@staticmethoddef identity_matrix(n):p_array = [None] * n# 生成每一个向量for k in range(0, n):# 创建数组p_array[k] = [1 if e == k else 0 for e in range(0, n)]return Matrix(p_array)def to_latex(self):s = '\\begin{pmatrix}'rows = len(self.__lines)columns = len(self.__lines[0])for i in range(0, columns):for j in range(0, rows):if j > 0:s += ' & 's += str(round(self.__lines[j][i] * 1000) / 1000)s += "\\\\\n"return s + "\end{pmatrix}\\\\\n"

8 代码测试

import unittestfrom com.youngthing.mathalgorithm.linearalgebra.hessenberg import Matrixclass MyTestCase(unittest.TestCase):def test_household_reflector(self):matrix = Matrix([[2, 1, -1, 7, 8], [1, 2, -1, 10, 11], [-1, -1, 2, 9, 6], [11, 3, 4, 5, 12], [16, 17, -4, -5, -6]])print(matrix.to_latex())print(matrix.to_hessenberg())print(matrix.eigen_values())def test_givens(self):matrix = Matrix([[2, 1, -1, 7, 8], [1, 2, -1, 10, 11], [-1, -1, 2, 9, 6], [11, 3, 4, 5, 12], [16, 17, -4, -5, -6]])print("A=", matrix.to_latex())h = matrix.to_hessenberg()print("H=", h.to_latex())print(h.givens_rotations())def test_qr(self):matrix = Matrix([[2, 1, -1, 7, 8], [1, 2, -1, 10, 11], [-1, -1, 2, 9, 6], [11, 3, 4, 5, 12], [16, 17, -4, -5, -6]])print("A=", matrix.to_latex())print(matrix.eigen_values())if __name__ == '__main__':unittest.main()

  测试结果，特征值为：

[22.315, -16.075, 4.631, -6.938, 1.067]

本文标签： 66 Hessenberg法求特征值