基于QR分解迭代求解方阵特征值和特征向量 _矩阵

基于QR分解迭代求解方阵特征值和特征向量一、特征值与特征向量求解的难点
线性代数的知识告诉我们如果要求一个方阵的特征值，只需要求解如下的特征方程的根即可：
f ( λ ) = ( λ ? λ 1 ) n 1 ( λ ? λ 2 ) n 2 ? ( λ ? λ s ) n s f(\)=(\-\)^{n_1}(\-\)^{n_2}\cdots(\-\)^{n_s} f(λ)=(λ?λ1?)n1?(λ?λ2?)n2??(λ?λs?)ns?
但是在具体程序中，如何去求解一个高次的多项式方程的根本身就是一个难点，它的实现甚至要比求得特征值还要复杂。因此，线性代数中这种用来手算的方法很明显不适合程序实现。那么，我们该如何去编写程序求解特征值呢？
二、QR迭代法求解特征值的原理
经过搜集相关资料，已知在计算机编程领域求解特征值一般有两种比较常见的算法，一种是雅可比()迭代法，另一种就是我们要讲的QR迭代法，这两种算法本质上都是去不断迭代，通过逼近的方法去求得特征值。雅可比方法有一个弊端，就是他只能求解实对称矩阵的特征值，适用面比较窄，但是QR分解迭代法可以求解任意矩阵的特征值，而且实现上更加接近我们线性代数课程的知识，因此下面介绍QR分解迭代法求解特征值的原理。
**理论依据：**任意一个非奇异矩阵（满秩的方阵） A A A都可以分解为一个正交矩阵 Q Q Q和一个上三角矩阵 R R R的乘积，且当 R R R对角元符号确定时，分解是唯一的。QR分解是一种迭代方法，迭代格式如下：
A k = Q k R k A k + 1 = R k + 1 Q k + 1 A_k=\\ A_{k+1}=R_{k+1}Q_{k+1} Ak?=Qk?Rk?Ak+1?=Rk+1?Qk+1?
当$ A_K $基本上收敛到上三角矩阵的时候，迭代完成，此时主对角线元素就是特征值。
**特别的：**当 A A A是实对称矩阵时， A K A_K AK?是对角阵 Λ \ Λ， Q = Q k ? 1 Q k ? 2 ? Q 1 Q=Q_{k-1}Q_{k-2}\cdots Q_1 Q=Qk?1?Qk?2??Q1?就是其标准正交的特征向量矩阵，且有 Q T A Q = A k = Λ Q^TAQ=A_k=\ QTAQ=Ak?=Λ 。
那么我们该如何去理解上述原理呢？下面进行简要的推导。
不妨令 A 1 = A A_1=A A1?=A，对 A 1 A_1 A1?进行QR分解得到:
A 1 = Q 1 R 1 R 1 = Q 1 T A 1 = Q 1 T A A_1=\\ R_1=Q_1^TA_1=Q_1^TA A1?=Q1?R1?R1?=Q1T?A1?=Q1T?A
令 A 2 = R 1 Q 1 = Q 1 T A Q 1 A_2==Q_1^TAQ_1 A2?=R1?Q1?=Q1T?AQ1?，很明显 A ∽ A 2 A\ A_2 A∽A2?，对 A 2 A_2 A2?再进行QR分解：
A 2 = Q 2 R 2 R 2 = Q 2 T A 2 A_2=\\ R_2=Q_2^TA_2 A2?=Q2?R2?R2?=Q2T?A2?
继续令 A 3 = R 2 Q 2 = Q 2 T A 2 Q 2 = Q 2 T Q 1 T A Q 1 Q 2 A_3==Q_2^=Q_2^TQ_1^ A3?=R2?Q2?=Q2T?A2?Q2?=Q2T?Q1T?AQ1?Q2?，可以看到 A ∽ A 3 A\ A_3 A∽A3?，同理，我们不断去进行QR分解，最终就会得到：
A k = Q k ? 1 T Q k ? 2 T ? Q 1 T A Q 1 ? Q k ? 2 Q k ? 1 A ∽ A k A_k=Q_{k-1}^TQ_{k-2}^T\cdots Q_1^TAQ_1\cdots Q_{k-2}Q_{k-1}\\ A\ A_k Ak?=Qk?1T?Qk?2T??Q1T?AQ1??Qk?2?Qk?1?A∽Ak?
因此， A A A与 A k A_k Ak?拥有相同的特征值。特别的，当 A = A T A=A^T A=AT时，有 A k = Λ A_k=\ Ak?=Λ，将 Q 1 ? Q k ? 2 Q k ? 1 Q_1\cdots Q_{k-2}Q_{k-1} Q1??Qk?2?Qk?1?记作 Q Q Q，那么 Q T Q^T QT和 Λ \ Λ就分别是正交的特征向量矩阵和对角线元素为特征值的对角阵。
有了上述理论介绍，我们已经知道了如何去进行QR迭代求特征值，下面马上又有一个问题摆在了我们面前——如何用程序去做QR分解呢？
三、QR分解的初等变换法
QR分解有三种常用的算法：Gram-正交化（线性代数讲授的方法）、初等变换法、变换法。本文介绍初代变换法。此处的初等变换指的是：
将矩阵的某一行（列）的 k k k倍加到另一行（列）
我们有以下引理：
引理1：设 A A A是一个实矩阵，若 A A A列满秩，则 A T A A^TA ATA对称正定，因此 A T A A^TA ATA有唯一的分解 A T A = L D L T A^TA=LDL^T ATA=LDLT，其中 L L L是单位下三角矩阵， D D D是对角元全部为正的对角矩阵。