文章目录
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}
1100011100111101111111111
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=
x1x2x3⋮xn
ρ=−sign(x1)z=
x1−ρ∥x∥x2x3⋮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=(100I−2u1u1T)A1=P1AP1⋮An−2=Pn−2An−1Pn−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=(I00I−2uiuiT)
需要注意的是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,1A3,1A4,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,1Ai+2,1Ai+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.01.02.0−1.010.011.0−1.0−1.02.09.06.011.03.04.05.012.016.017.0−4.0−5.0−6.0
A1=
2.0−10.7240.00.00.0−19.3034.061−0.96611.7069.1290.732−11.6982.8952.572−1.02−1.122−0.7344.4443.0557.4952.14618.797−4.01−3.602−7.01
A2=
2.0−10.7240.00.00.0−19.3034.06114.8760.00.00.38611.7170.986−0.0761.583−0.867−18.0366.884.2534.982.3455.305−6.7470.758−6.3
A3=
2.0−10.7240.00.00.0−19.3034.06114.8760.00.00.38611.7170.9861.5850.02.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)eieiT−cos(θ)eiejT+sejejT+(cos(θ)−1)ejejT=
I00000cos(θ)0sin(θ)000I000−sin(θ)0cos(θ)00000I
=
I00000c0s000I000−s0c00000I
公式中的矩阵的
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,i2c=cos(θ)=ri(Hi)i,is=sin(θ)=ri(Hi+1)i,iGi=G(i,i+1,c,s)Hi+1=GiTHiH=（i=1∏n−1Gi）Hn=QRQ=i=1∏n−1GiR=HnH′=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.01.02.0−1.010.011.0−1.0−1.02.09.06.011.03.04.05.012.016.017.0−4.0−5.0−6.0
H=
2.0−10.7240.00.00.0−19.3034.06114.8760.00.00.38611.7170.9861.5850.02.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.00.38611.7170.9861.5850.02.3846.163−7.069−6.550.24−0.754−17.7616.5494.4624.503
G1=
0.183−0.9830.00.00.00.9830.1830.00.00.00.00.01.00.00.00.00.00.01.00.00.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.2417.322−3.9976.5494.4624.503
G2=
1.00.00.00.00.00.0−0.7750.6320.00.00.0−0.632−0.7750.00.00.00.00.01.00.00.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.2417.3227.237−2.5474.4624.503
G3=
1.00.00.00.00.00.01.00.00.00.00.00.0−0.830.5570.00.00.0−0.557−0.830.00.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.2417.3227.2374.601−2.2864.503
G4=
1.00.00.00.00.00.01.00.00.00.00.00.01.00.00.00.00.00.00.9980.0660.00.00.0−0.0660.998
Q=
0.183−0.9830.00.00.0−0.762−0.1420.6320.00.00.5160.0960.6430.5570.00.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.017.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.01.02.0−1.010.011.0−1.0−1.02.09.06.011.03.04.05.012.016.017.0−4.0−5.0−6.0
H0=
2.0−10.7240.00.00.0−19.3034.06114.8760.00.00.38611.7170.9861.5850.02.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.04.1597.3546.812−3.130.30817.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.5440.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.06.9168.7368.4872.5471.2612.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.01.374−15.274−8.7214.3180.6454.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.09.82213.4298.7944.4970.1776.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.04.237−16.212−8.6144.6570.0424.512−4.7492.1731.4021.08
H8=
22.763−4.650.00.00.03.744−16.5230.0050.00.0−0.233−10.042−6.9130.0340.08.52214.5578.694.6060.015.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 cmath
MIN_VALUE = 1e-2
def round_num(x):
if isinstance(x, complex):
return complex(round_num(x.real), round_num(x.imag))
return round(x * 1000) / 1000.0
class Matrix:
# 矩阵，为了方便计算，二维数组的每个一维数组表示一个向量，也就是一列
# 所以看起来会比较奇怪
def __init__(self, lines):
self.__lines = lines
def __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] += mul
return 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 s
def eigen_values(self):
h = self.to_hessenberg()
while not h.is_upper_triangle():
q, r = h.givens_rotations()
h = r * q
return 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
# 这时候再判断是不是拟三角矩阵
# 如果左上角是不是0
if i > 0 and abs(self.__lines[i - 1][j - 1]) > MIN_VALUE:
return False
# 如果右下角是不是0
if i < len(self.__lines) - 1 and abs(self.__lines[i + 1][j + 1]) > MIN_VALUE:
return False
return True
def __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
# 如果是下一个元素不是0
if 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 * c
if root > 0:
sqrt = math.sqrt(root)
else:
sqrt = cmath.sqrt(root)
a_d = a + d
eigen_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.0
return eigen_values
def givens_rotations(self):
h = self
n = len(self.__lines)
q = None
for i in range(0, n - 1):
g, h = h.givens_rotation(i)
q = g if q is None else q * g
return q, h
def givens_rotation(self, i):
vector = self.__lines[i]
n = len(self.__lines)
j = i + 1
r = math.sqrt(vector[i] ** 2 + vector[j] ** 2)
c = vector[i] / r
s = vector[j] / r
g = Matrix.identity_matrix(n)
g.__lines[i][i] = c
g.__lines[i][j] = s
g.__lines[j][i] = -s
g.__lines[j][j] = c
g_t = Matrix.identity_matrix(n)
g_t.__lines[i][i] = c
g_t.__lines[i][j] = -s
g_t.__lines[j][i] = s
g_t.__lines[j][j] = c
h = g_t * self
return g, h
def to_hessenberg(self):
n = len(self.__lines)
a = self
for i in range(0, n - 2):
p = a.household_reflector(i)
a = p * a * p
return a
def household_reflector(self, i):
j = i + 1
column = self.__lines[i]
rho = 1 if column[j] < 0 else -1
z = [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左边的就是矩阵I
sub_matrix = Matrix([u]) * Matrix(ut)
sub_p = sub_matrix.__lines
# i 和 i以后的需要处理
# 这里需要仔细处理
# 对于i = 0, j = 1
for 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
@staticmethod
def length(vector):
return math.sqrt(Matrix.square(vector))
@staticmethod
def square(vector):
return Matrix.dot(vector, vector)
@staticmethod
def dot(vector1, vector2):
l = 0
for i in range(0, len(vector1)):
l += vector1[i] * vector2[i]
return l
@staticmethod
def 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 unittest
from com.youngthing.mathalgorithm.linearalgebra.hessenberg import Matrix
class 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]
扫一扫
1万+
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
