矩阵特征值的数值解法

特征值的数值解

https://www.google.com/url?sa=t&rct=j&q=&esrc=s&source=web&cd=&cad=rja&uact=8&ved=2ahUKEwiKtIiKo5byAhU-ITQIHQDPDLAQFnoECC4QAw&url=https%3A%2F%2Fwww.math.tamu.edu%2F~dallen%2Flinear_algebra%2Fchpt6.pdf&usg=AOvVaw34g7kDAuOx6I6ibjzOVzkr

特征值问题的数值解的稳定性

矩阵元素的微小变化，并不会导致特征值的微小变化。如果矩阵有微小变化且矩阵对称，那么特征值也有微小变化。因此在应用中，对称矩阵的数值解非常精确。

stability theorem：矩阵A，误差矩阵E，A与E为实对称矩阵，则A的特征值与A+E的特征值以及E的关系为(λ1−λ1^)2+(λ2−λ2^)2+...+(λn−λn^)2≤∣∣E∣∣F2(\lambda_1-\hat{\lambda_1})^2+(\lambda_2-\hat{\lambda_2})^2+...+(\lambda_n-\hat{\lambda_n})^2 \leq ||E||^2_F(λ1−λ1^)2+(λ2−λ2^)2+...+(λn−λn^)2≤∣∣E∣∣F2。

由定理可得推论 stability corollary：(λk−λk^)2≤(λ1−λ1^)2+(λ2−λ2^)2+...+(λn−λn^)2≤∣∣E∣∣F2(\lambda_k-\hat{\lambda_k})^2 \leq (\lambda_1-\hat{\lambda_1})^2+(\lambda_2-\hat{\lambda_2})^2+...+(\lambda_n-\hat{\lambda_n})^2 \leq ||E||^2_F(λk−λk^)2≤(λ1−λ1^)2+(λ2−λ2^)2+...+(λn−λn^)2≤∣∣E∣∣F2，因此(λk−λk^)2≤∣∣E∣∣F2(\lambda_k-\hat{\lambda_k})^2 \leq ||E||^2_F(λk−λk^)2≤∣∣E∣∣F2。即如果矩阵有微小变化且矩阵对称，那么特征值绝对误差很小。
但尽管绝对误差∣λk−λk^∣|\lambda_k-\hat{\lambda_k}|∣λk−λk^∣很小，相对误差∣λk−λk^∣∣λk∣\frac{|\lambda_k-\hat{\lambda_k}|}{|\lambda_k|}∣λk∣∣λk−λk^∣也可能很大。

幂方法power method

在实际中想要找到A的近似矩阵A+E很难，因此在实际中会针对不同的矩阵如对称，三对角等一些特殊矩阵分别设计各自的算法。

满足两个条件的矩阵A，特征向量为V1,...,Vn{V_1,...,V_n}V1,...,Vn
X0=c1V1+...+cnVnX_0 = c_1V_1 + ... + c_nV_nX0=c1V1+...+cnVn
AX0=c1λ1V1+...+cnλnVnAX_0 = c_1\lambda_1 V_1 + ... + c_n\lambda_n V_nAX0=c1λ1V1+...+cnλnVn
AmX0=c1λ1mV1+...+cnλnmVnA^mX_0 = c_1{\lambda_1}^m V_1 + ... + c_n{\lambda_n}^m V_nAmX0=c1λ1mV1+...+cnλnmVn
当m很大时，AmX0λ1m≈c1V1\frac{A^mX_0}{\lambda_1^m} \approx c_1V_1λ1mAmX0≈c1V1
c1≠0c_1\neq 0c1=0，且X0X_0X0不与V1V_1V1正交

AmX0λ1mY≈c1V1Y\frac{A^mX_0}{\lambda_1^m}Y \approx c_1V_1Yλ1mAmX0Y≈c1V1Y
Am+1X0λ1m+1Y≈c1V1Y\frac{A^{m+1}X_0}{\lambda_1^{m+1}}Y \approx c_1V_1Yλ1m+1Am+1X0Y≈c1V1Y
YYY不与V1V_1V1正交，且AmX0λ1mY=Am+1X0λ1m+1Y≠0\frac{A^mX_0}{\lambda_1^m}Y =\frac{A^{m+1}X_0}{\lambda_1^{m+1}}Y \neq 0λ1mAmX0Y=λ1m+1Am+1X0Y=0
两式相除得Am+1X0YAmX0Y≈λ1\frac{A^{m+1}X_0Y}{A^mX_0Y} \approx \lambda_1AmX0YAm+1X0Y≈λ1

幂方法的由来就是因为A的幂。用于幂方法需要先设计出X0X_0X0与YYY。
m从1开始往下不断计算，直到满足∣λ1calc−λ1actual∣|\lambda_1^{calc} - \lambda_1^{actual}|∣λ1calc−λ1actual∣很小，但我们并不知道λ1actual\lambda_1^{actual}λ1actual。
对于实对称矩阵∣λ1calc−λ1actual∣≤AX∗AXX∗X−(λ1calc)2,X=AmX0,λ1calc=AX∗XX∗X|\lambda_1^{calc} - \lambda_1^{actual}| \leq \sqrt{\frac{AX*AX}{X*X}-(\lambda_1^{calc})^2}, X=A^mX_0,\lambda_1^{calc}=\frac{AX*X}{X*X}∣λ1calc−λ1actual∣≤X∗XAX∗AX−(λ1calc)2,X=AmX0,λ1calc=X∗XAX∗X，误差阈值需要事先确定，当AX∗AXX∗X−(λ1calc)2\sqrt{\frac{AX*AX}{X*X}-(\lambda_1^{calc})^2}X∗XAX∗AX−(λ1calc)2小于阈值时，停止计算。
对于实非对称矩阵En+1=λ1calcn−λ1calcn+1λ1calcn+1E_{n+1} = \frac{\lambda_1^{calc n} - \lambda_1^{calc n+1}}{\lambda_1^{calc n+1}}En+1=λ1calcn+1λ1calcn−λ1calcn+1，误差阈值需要事先确定，当En+1E_{n+1}En+1小于阈值时，停止计算。

计算出主特征值后，计算下一个特征值。
B=A−λ1V1∣V1∣V1∣V1∣TB = A- \lambda_1 \frac{V_1}{|V_1|} \frac{V_1}{|V_1|}^TB=A−λ1∣V1∣V1∣V1∣V1T，B的特征值为0,λ2,λ3,...,λn0, \lambda_2,\lambda_3,...,\lambda_n0,λ2,λ3,...,λn
B的特征值与A一样，对B使用幂方法就可以算出A的下一个特征值。
这种方法称为method of deflation。但是还需要求出主特征值对应的特征向量。

特征对eigenpair：特征值与对应的特征向量组成特征对。

QR method

A=QRA=QRA=QR，Q为正交矩阵，R为上三角矩阵，Q为实数方阵
QT=Pn−1Pn−2...P2P1Q^T = P_{n-1}P_{n-2}...P_2P_1QT=Pn−1Pn−2...P2P1, R=Pn−1Pn−2...P2P1AR=P_{n-1}P_{n-2}...P_2P_1AR=Pn−1Pn−2...P2P1A，P为正交矩阵也是householder矩阵

这个定义我没看懂，不影响对方法的理解。
$A_1,...,A_m,Q_1,...,Q_m,R_1,...,R_m$的定义如下：
step1: $A_1=A, Q_1=Q, R_1=R$
step2: $A_2=R_1Q_1, A_2=Q_2R_2$
step3: $A_3=R_2Q_2, A_3=Q_3R_3$
step4: $A_m=R_{m-1}Q_{m-1}, A_m=Q_mR_m$

m越大，AmA_mAm越接近一个三角矩阵，AmA_mAm的特征值越接近A

通过PkP_kPk将A逐步化为上三角矩阵

Pn−1..P1AP_{n-1}..P_1APn−1..P1A为一个上三角矩阵

PkP_kPk的构造

householder矩阵：H=I−2VVT,∣V∣=1H=I-2VV^T, |V|=1H=I−2VVT,∣V∣=1，householder矩阵也是正交矩阵