QR算法的变体
Cleve Moler, MathWorks
QR算法是数学软件中最成功和最强大的工具之一。
MATLAB®核心库包括QR算法的几个变体。这些变量计算实对称矩阵、实非对称矩阵、实矩阵对、复矩阵、复矩阵对和各种类型矩阵的奇异值的特征值。这些核心库函数在各种MATLAB工具箱中用于查找稀疏矩阵和线性算子的特征值和奇异值,查找多项式的零点,求解特殊线性系统,评估数值稳定性,以及执行许多其他任务。
数十人对QR算法变体的开发做出了贡献。关于这一主题的第一篇论文是由J.G.F.弗朗西斯(1961年和1962年)和Vera N. Kublanovskaya(1963年)发表的。但是,是J.H.威尔金森开发了QR算法的第一个完整实现。威尔金森还提出了一个重要的收敛分析。威尔金森的书代数特征值问题他的两篇论文发表于1965年。这意味着我们将能够庆祝2015年作为实用QR算法的黄金周年。
用于奇异值分解(SVD)的QR算法的变体由Gene Golub和Velvel Kahan在1965年发表,并由Golub和Christian Reinsch在1969年完善。
名称“QR”
“QR”的名字来源于字母Q,用来表示正交矩阵,字母R,用来表示直角三角形矩阵。有一个qr
函数,但它计算QR分解,而不是QR算法。任何矩阵,无论是实矩阵还是复矩阵,方阵还是方阵,都可以因式分解成一个矩阵的乘积问用标准正交列和矩阵R它只有在它的上三角形或右三角形中是非零的。你可能还记得Gram Schmidt过程,它做了几乎相同的事情,尽管它最初的形式在数值上是不稳定的。
一行程序
使用qr
函数,QR算法的简单变体可以在一行MATLAB代码中表示。让一个
是方形的,n
——- - - - - -n
矩阵,让I = eye(n,n)
.然后给出QR迭代的一步
s = A(n,n);[Q,R] = qr(A - s*I);A = R*Q + s*I
量\(s\)是位移;它加速了融合。当\(A\)趋于上三角矩阵时,\(s\)趋于特征值。如果在一行中输入这三条语句,则可以使用向上箭头键进行迭代。
QR分解得到一个上三角形R。
\[A-sI = QR \]
然后逆序乘法,\(RQ\),恢复特征值,因为
\[RQ + sl = Q^{T}(A - sl)Q + sl = Q^{T} Q \]
所以新的\(A\)和原来的\(A\)很相似。每次迭代都有效地将一些“质量”从下三角形转移到上三角形,同时保留特征值。随着迭代的进行,矩阵开始接近上三角矩阵,特征值方便地显示在对角线上。
一个例子
为了说明这个过程,我们将使用MATLAB Gallery集合中的一个矩阵。
A =画廊(3)A = -149 -50 -154 537 180 546 -27 -9 -25
这一点都不明显,但是这个矩阵已经被构造成具有特征值1 2和3。我们的单行QR码的第一次迭代产生
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 -户主反射,这是一个相似变换序列,保留特征值。QR迭代适用于三对角线形式。威尔金森提供了一种转移策略,使他能够证明全局收敛和局部三次收敛率。即使存在舍入误差,该算法也能保证成功。
图1显示了一个初始对称矩阵,三对角化简一半的情况,三对角化简,QR迭代一半的情况,最后是特征值。实际上,由于矩阵是对称的,计算只在数组的一半上执行,但我们的图反映了结果,以显示整个矩阵。
非对称矩阵
实的非对称矩阵的情况要复杂得多。最初减少使用n -户主相似度转换来创建一个海森伯格矩阵,它是上三角加上一个“额外的”次对角线。然后使用具有双移位策略的QR迭代。在尝试创建一个真正的舒尔形式时,这保留了海森伯格形式,它是上三角形,除了对角线上对应于复共轭特征值对的2 × 2块。
非对称Hessenberg QR算法并非绝对正确。这是一个迭代过程,并不总是保证收敛。甚至在30年前,基本迭代的反例就已经为人所知。威尔金森引入了一个额外的“特别”位移来处理它们,但还没有人能够证明一个完整的收敛定理。所以,在很少的情况下,MATLAB用户可能会看到这样的消息:
错误使用==> eig,解决方案将不收敛
几年前,收到这条信息的人可能会认为这是不可避免的。但今天,大多数人会感到惊讶或恼火;他们已经开始期望绝对正确。
我们现在知道一个4乘4的例子,它可能会导致真实的、非对称的QR算法在某些计算机上失败,即使是威尔金森的临时移位。矩阵是
\[A = \begin{bmatrix} 0 & 2 & 0 & -1 \\ 1 & 0 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 \\ 0 & 0 \ 0 & 1 & 0 \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 & 0 & 1 & 0 \end{bmatrix} \]
其中\(\delta\)很小,但又不至于小到可以忽略——比如,\(\delta = 10^{-8}\)。准确的特征值接近于一对二重根。Wilkinson双移迭代使用每对中的一个特征值。这种迭代确实改变了矩阵,但不足以获得快速收敛。所以我们必须使用不同的双移基于重复一个低的2 × 2块的特征值。
这些不同的临时移位策略已被纳入LAPACK的最新版本,并因此纳入MATLAB。我们现在的情况是,我们不知道任何矩阵会引起eig
或根
要显示“不会收敛
消息,但是我们没有证据证明这些不同修饰的非对称代码是绝对正确的。
图2显示了一个初始的非对称矩阵,一半的情况下,约简到Hessenberg形式,Hessenberg形式,一半的情况,通过QR迭代,和最终的实Schur形式。对于这个特殊的矩阵,它恰好有四个实数特征值和三个复共轭对,总共有十个特征值。
奇异值
一个可能的矩形矩阵\(a \)的奇异值是对称矩阵\(a ^{T} a \)特征值的平方根。这个事实可以用来激励和分析一个算法,但它不应该是实际计算的基础与有限精度算术。Golub-Kahan-Reinsch算法的初始阶段涉及在左侧和右侧操作的Householder反射,以将矩阵简化为双对角线形式。这个阶段之后是QR算法在双对角线上运行的SVD变体。将Wilkinson对对称三对角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数值计算,第10章,“特征值和奇异值”。
硅藻土、克里夫1976年矩阵奇异值分解薄膜.
威尔金森,J.H代数特征值问题,牛津大学出版社,1965年,662页。
注:本文的一部分是基于“QR算法”MATLAB新闻和笔记, 1995年夏。
发布于2014 - 92224v00