QR算法的变体

QR算法是数学软件中最成功和最强大的工具之一。

MATLAB^®核心库包括QR算法的几个变体。这些变体计算实对称矩阵、实非对称矩阵、实矩阵对、复矩阵对、复矩阵对以及各种类型矩阵的奇异值。这些核心库函数被用于各种MATLAB工具箱中，用于查找稀疏矩阵和线性算子的特征值和奇异值，查找多项式的零，求解特殊线性系统，评估数值稳定性，以及执行许多其他任务。

数十人参与了QR算法变体的开发。关于这一课题的第一篇论文来自1961年和1962年的J.G.F. Francis和1963年的Vera N. Kublanovskaya。但是，是J.H.威尔金森开发了第一个完整的QR算法的实现。威尔金森还发展了一个重要的收敛分析。威尔金森的书代数特征值问题他的两篇论文发表于1965年。这意味着我们将2015年作为实用QR算法的金婚纪念日来庆祝。

用于奇异值分解(SVD)的QR算法的变体由Gene Golub和Velvel Kahan在1965年发表，并在1969年由Golub和Christian Reinsch完善。

“QR”这个名字

QR是由表示正交矩阵的字母Q和表示直角三角形矩阵的字母R而来。有一个qr函数，但它计算的是QR分解，而不是QR算法。任何矩阵，无论是实矩阵还是复矩阵，方矩阵还是矩形矩阵，都可以分解成矩阵的乘积问用标准正交的列和矩阵R它只有在上面的三角形，或者说是右边的三角形，是非零的。你可能还记得格拉姆施密特过程，它做了几乎相同的事情，尽管在最初的形式中它在数值上是不稳定的。

一行程序

使用qr函数，QR算法的一个简单变体可以用一行MATLAB代码表示。让一个是一个广场,n——- - - - - -n矩阵,让我眼睛= (n, n)．然后给出QR迭代的一步

s = (n, n);[Q,R] = qr(A - s*I)

量\(s\)是移位;它加速收敛。当\(A\)趋于上三角矩阵时，\(s\)趋于特征值。如果在一行中输入这三个语句，可以使用向上箭头键进行迭代。

QR分解产生一个上三角\(R\)。

\[A-sI = QR \]

然后逆序乘法，\(RQ\)，恢复特征值，因为

\[RQ + sl = Q^{T}(A - sl)Q + sl = Q^{T}AQ \]

所以新的\(A\)和原来的\(A\)很相似。每次迭代有效地将一些“质量”从下三角形转移到上三角形，同时保留特征值。随着迭代的进行，矩阵开始接近上三角矩阵，特征值方便地显示在对角线上。

一个例子

为了说明这个过程，我们将使用MATLAB Gallery集合中的一个矩阵。

A = -149 -50 -154 537 180 546 -27 -9 -25

这一点也不明显，但是这个矩阵已经被构造成具有特征值1、2和3。我们的一行二维码的第一次迭代产生

A = 28.8263 -259.8671 773.9292 1.0353 -8.6686 33.1759 -0.5973 5.5786 -14.1578

矩阵现在更接近于上三角，但特征值仍然不明显。然而，在五个迭代之后，我们有了

A = 3.0321 -8.0851 804.6651 0.0017 0.9931 145.5046 -0.0001 0.0005 1.9749

我们开始看到特征值3 1 2出现在对角线上。再进行8次迭代

A = 3.0716 -7.6952 802.1201 0.0193 0.9284 158.9556 00 2.0000

将特征值2.0计算到显示精度，与之相邻的对角下元素为零。此时，有必要继续对左上2 × 2子矩阵进行迭代。

QR算法从来不会以这种简单的形式执行。在它之前，总是要将其简化成紧凑形式，其中次对角线以下的所有元素都为零。迭代保留了这种简化形式，并且分解可以更快地完成。移位策略更为复杂，对于不同形式的算法也不同。此外，简化形式对于迭代的收敛性是至关重要的。

对称矩阵

图1-3展示了QR算法的三个最重要的变体。这些图是从程序生成的输出中获取的快照eigsvdgui.m从MATLAB的数值计算．

最简单的变体涉及实对称矩阵。一个n——- - - - - -n实对称矩阵可以被简化为三对角形式n -Householder反射，它是一个保持特征值的相似变换序列。QR迭代适用于三对角形式。威尔金森提供了一个转移策略，使他能够证明全局收敛和局部三次收敛速度。即使存在舍入误差，该算法也能保证成功。

图1给出了一个初始对称矩阵，在约简到三对角线的情况下，在三对角线的情况下，在QR迭代的情况下，最后是特征值。实际上，由于矩阵是对称的，计算只在数组的一半上执行，但我们的图反映的结果显示了整个矩阵。

非对称矩阵

实非对称矩阵的情况要复杂得多。初始减量用途n -户主相似变换创建一个海森伯格矩阵，它是上三角形加上一个“额外的”次对角线。然后使用双移策略的QR迭代。这保留了海森伯格形式，同时试图创建一个实舒尔形式，这是上三角形，除了2 × 2块对应于对角线上的复共轭特征值对。

非对称海森伯格QR算法并不是绝对可靠的。这是一个迭代过程，并不总是保证收敛。甚至在30年前，基本迭代的反例就已经为人所知。威尔金森引入了一个额外的“临时”转换来处理它们，但没有人能够证明一个完整的收敛定理。因此，在极少数情况下，MATLAB用户可能会看到这样的消息:

错误使用==> eig，解不会收敛

几年前，收到这条信息的人可能会认为这是不可避免的。但今天，大多数人会感到惊讶或烦恼;他们已经开始期望万无一失。

我们现在知道了一个4乘4的例子，它可能会导致真实的非对称QR算法在某些计算机上失败，即使有威尔金森的临时移位。矩阵是

\[A = begin{bmatrix} 0 & 2 & 0 & -1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 0 \\ 1 \end{bmatrix}\]

这是多项式(p(x)=x^4 - 2x^2 + 1)的同伴矩阵，这个陈述

根([1 0 -2 0 1])

的计算eig (A)．值\(\lambda = 1 \)和\(lambda = -1 \)都是特征值或多项式根，重数为2。对于真实的\(x\)，多项式\(p(x)\)永远不为负。这些二重根极大地减慢了迭代速度，以至于在某些计算机上，在检测到收敛之前，舍入误差的变化会干扰迭代。迭代可能永远徘徊，试图收敛，但当它接近时就偏离方向。

该表单的示例显示了类似的行为

\[begin{bmatrix}0 & 1 & 0 & 0 \\ 1 & 0 & -\delta & 0 \\ 0 & \delta & 0 & 1 & 0 \end{bmatrix} \]

其中，\(\delta\)很小，但又不至于小到可以忽略——比如，\(\delta = 10^{-8}\)。精确的特征值接近于一对二重根。威尔金森双移迭代使用每对的一个特征值。这种迭代确实改变了矩阵，但不足以得到快速的收敛。所以我们必须使用不同的双移基于重复下面2 × 2块的一个特征值。

这些不同的临时转换策略已经被整合到LAPACK的最新版本中，因此也被整合到MATLAB中。我们现在的情况是，我们不知道任何矩阵导致eig或根显示“不收敛信息，但我们没有证据证明这些不同修饰的非对称代码是绝对可靠的。

图2显示了一个初始非对称矩阵、还原到Hessenberg形式中途的情况、Hessenberg形式、QR迭代中途的情况以及最终的实Schur形式。对于这个特殊的矩阵，有四个实特征值和三个复共轭对，总共十个特征值。

奇异值

一个可能的矩形矩阵的奇异值是对称矩阵的特征值的平方根。这一事实可以用来激励和分析算法，但不应作为有限精度算法实际计算的基础。Golub-Kahan-Reinsch算法的初始阶段包括Householder反射在左右两边操作，以将矩阵缩减为双对角形式。这一阶段之后是QR算法的SVD变体在双对角线上操作。该算法也应用了对称三对角QR的威尔金森分析，保证了算法的全局收敛性。

图3给出了初始矩形矩阵、约简到双对角线形式的情况、双对角线形式、QR迭代到一半的情况以及包含奇异值的最终对角线形式。

QR算法应用

虽然QR算法计算特征值和奇异值是密切相关的，但其结果的应用通常是非常不同的。特征值经常被用来分析常微分方程组，其中行为作为时间的函数是重要的。另一方面，奇异值对于分析联立线性方程组的静态系统是有用的，其中方程组的数目往往不等于未知数的数目。

控制理论和控制设计自动化大量使用特征值。控制理论中研究的经典状态空间微分方程组是

{对齐}\ \[\开始点{x} & y = Ax +布鲁里溃疡\ \ & =残雪+ Du \{对齐}结束\]

利用QR算法计算\(A\)的特征值对于研究稳定性和可控性等问题至关重要。

在统计学中，SVD是一种获得主成分的数值可靠方法。主成分分析(PCA)是一种分析联立线性方程组的超定系统的技术

\ [Ax = b \]

其中\(A\)的行多于列。利用QR算法计算\(A\)的奇异值和向量，得到主分量。

为进一步阅读和查看

Golub, Gene H.和Charles F. Van Loan，矩阵计算，第四版，约翰霍普金斯大学出版社，1996,697页。

硅藻土、克里夫MATLAB的数值计算，第十章“特征值和奇异值”。

硅藻土、克里夫1976矩阵奇异值分解薄膜．

威尔金森,J.H代数特征值问题，牛津大学出版社，1965年，662页