使用QR法求解矩阵特征值

求解特征值问题的很多方法都是基于一系列正交相似变换。因为如果P是任意非奇异矩阵，那么矩阵A和PAP^-1具有相同的特征值。更进一步，如果v是A的一个特征向量，那么Pv便是PAP^-1的一个特征向量。这样，如果P是正交矩阵，特征值问题的条件数就不会受到影响。这便是相似变换法的基础。

QR法^{[20，52，53]}是矩阵特征值计算中得到最广泛应用的分解法之一。它使用一系列正交相似变换^[13，28]，按如下公式

构造矩阵迭代系列A=A₀，A₁，A₂，…（译者注：一般情况下，此矩阵序列“基本”收敛于一个上三角阵或分块上三角阵，即主对角线或主对角线子块以下元素均收敛于0，从而可以直接读取矩阵的特征值）。

与LU分解法类似，矩阵A也可以被分解为两个矩阵：

A=QR （7.8）

式中，Q是一个酉矩阵，R是一个上三角矩阵。酋矩阵Q为满足下式的矩阵：

QQ^∗=Q^∗Q=I （7.9）

式中，（∗）表示共轭转置。

酋矩阵的例子有

同时，酋矩阵的逆就是其共轭转置矩阵

Q^-1=Q^∗

这个分解导致A的列向量[a₁，a₂，…，a_n]与Q的列向量[q₁，q₂，…，q_n]存在如下关系：

列向量a₁，a₂，…，a_n必须从左到右依次正交归一化为一组标准正交基q₁，q₂，…，q_n。

QR法在具体实施时，常见的做法是先将矩阵A变换为一个具有相同特征值的Hessenberg矩阵H，然后对矩阵H使用QR法。这样，矩阵H最终被变换为一个“基本”上三角矩阵，特征值就可以从其对角线上读出。Hessenberg矩阵大致上是一个上三角矩阵，只是其主对角线下方的次对角线元素非零。将矩阵A先约简为Hessenberg矩阵，再对Hessenberg矩阵实施QR算法的原因是，这样做比直接对矩阵A实施QR算法计算量要小得多。

将矩阵A约简为Hessenberg矩阵的方法之一是Householder法。可以证明，对每个n×n矩阵A，存在n-2个Householder矩阵H₁，H₂，…，H_n-2，使得

Q=H_n-2…H₂H₁

且

P=Q^∗AQ

是一个Hessenberg矩阵^[27]。矩阵H是一个Householder矩阵，如果满足

注意，Householder矩阵同样也是酋矩阵。选择向量v使其满足

v_i=a_i±e_i‖a_i‖₂ （7.11）

式中，符号的正负选择应该使得‖v‖₂不会太小，e_i是I的第i列，a_i是A的第i列。

例7.2 对矩阵A进行QR分解

解7.2 施加第1次变换的目的是使矩阵A的第1列在主对角线以下的元素变为0，因此根据上面的公式有

从而有

第2次变换将针对上述约简后的矩阵中去除了第1行和第1列的部分进行，

因此

从而得到

继续这个过程得到(www.xing528.com)

从而得到

可以证明A=QR，并且Q^∗=Q^-1。

QR分解中的消去运算可以认为是Gauss消去法的一种替代算法。然而，其所需的乘法和除法次数是Gauss消去法的两倍。因此，QR分解很少被用来求解线性方程组，但它在特征值计算中发挥了重要的作用。

尽管特征值问题归结为一组简单代数方程的求解问题，

det（A-λI）=0

但实际上求解这组代数方程是困难的。计算特征方程的根或一个矩阵的零空间并不很适合于计算机。事实上，不存在一般性的直接方法可以在有限步内求出特征值。因此必须通过迭代计算的方法来产生一系列逐渐逼近矩阵特征值的近似值。

QR法通常被用来计算满矩阵的特征值和特征向量。如Francis所提出的那样[13]，QR法构造一系列相似变换

A_k=Q_k^∗A_k-1Q_kQ_k^∗Q_k=I （7.12）

式中，A_k与矩阵A相似。对A重复进行QR分解使得副对角线上的元素逐渐变为0。在最终的收敛结果中，特征值按照模值降序出现在A_k的对角线上。

例7.3 求例7.2中的矩阵的特征值和特征向量。

解7.3 第1个目标是要通过QR法求出矩阵A的特征值。根据例7.2已有的结果，得到第1次QR分解产生的矩阵Q₀为

将给定矩阵A标为A₀，第1次更新的矩阵A₁由下式计算得到

对A₁进行QR分解得到Q1

矩阵A₂为

注意，在对角线下的元素逐渐减小到0。这个过程一直进行下去，最终得到的矩阵A^∗为

特征值在A^∗的对角线上，并且是按模值降序排列的。因此特征值是

下一步是求出对应于每个特征值的特征向量。考虑到对每个特征值及其对应的特征向量有

Av_i=λ_iv_ii=1，…，n （7.14）

式（7.14）也可以被写成

Av_i-λ_iv_i=0

也就是说，矩阵A-λ_iI是奇异矩阵。因此，它只有3行（或列）是独立的。根据这个事实，一旦求出特征值，就可以确定特征向量。既然A-λ_iI不是满秩的，则特征向量v_i的元素之一可以任意选择。因此第1步，将A-λ_iI进行分块：

式中，a₁₁是一个标量，a_1，2n是一个1×（n-1）的向量，a_2n，1是一个（n-1）×1的向量，a_2n，2n是一个（n-1）×（n-1）的秩为（n-1）的矩阵。然后令v_i（1）=1，并求解特征向量的剩余部分。

现在更新v_i（1）为

这样对应于λ_i的特征向量为

最后一步是对特征向量进行归一化，因此

因此，对应于特征值构成的向量

Λ=[18.0425-6.4172-0.5269-0.0983]

相对应的特征向量为