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仅在其上部或右侧三角形中为非零。你可能还记得Gram-Schmidt过程,它做了几乎相同的事情,尽管它的原始形式在数值上是不稳定的。
一行程序
使用qr
函数,QR算法的一个简单变体可以用一行MATLAB代码表示。让一个
是一个广场,n
-借-n
矩阵,让I=眼睛(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}AQ\]
因此,新的\(A\)与原来的\(A\)相似。每次迭代有效地将一些“质量”从下三角形转移到上三角形,同时保留特征值。随着迭代的进行,矩阵开始接近上三角矩阵,特征值方便地显示在对角线上。
一个例子
为了说明这个过程,我们将使用MatlabGallery集合中的矩阵。
A=画廊(3)A=-149-50-154537180546-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出现在对角线上。再进行八次迭代就可以得到
A=3.0716-7.6952 802.1201 0.0193 0.9284 158.9556 0 0 2.0000
将特征值2.0计算到显示精度,与之相邻的对角下元素为零。此时,有必要继续对左上2 × 2子矩阵进行迭代。
QR算法从来没有以这种简单的形式进行过。它总是先简化为一种紧凑的形式,在这种形式中,次对角线下的所有元素都为零。迭代保留了这种简化的形式,并且可以更快地进行因式分解。移位策略更为复杂,并且对于各种形式的变换也有所不同此外,简化形式对于迭代的收敛性至关重要。
对称矩阵
图1-3展示了QR算法最重要的三种变体。这些数字是从程序生成的输出中获取的快照eigsvdgui.m
从…起用MATLAB进行数值计算.
最简单的变体涉及实对称矩阵。一n-借-n实对称矩阵可以通过下列方法简化为三对角形式:n -Householder反射,这是一系列保留特征值的相似变换。QR迭代适用于三对角形式。Wilkinson提供了一种移位策略,允许他证明全局收敛和局部三次收敛速度。即使存在舍入误差,该算法也能保证成功了。
图1显示了一个初始对称矩阵,在还原为三对角的过程中的情况,三对角,在QR迭代过程中的情况,最后是特征值。实际上,由于矩阵是对称的,因此计算仅在阵列的一半上执行,但我们的图反映了显示整个矩阵的结果。
非对称矩阵
实非对称矩阵的情况要复杂得多。初始减量用途n -Householder相似变换创建Hessenberg矩阵,该矩阵为上三角形加上“额外”次对角。然后使用带有双移位策略的QR迭代。这保留了Hessenberg形式,同时尝试创建一个真实的Schur形式,它是上三角形,除了与对角上的复共轭特征值对相对应的2×2块。
非对称Hessenberg QR算法并非绝对正确。它是一个迭代过程,并不总是保证收敛。即使在30年前,基本迭代的反例也是众所周知的。Wilkinson还引入了一个额外的“特殊”shift来处理它们,但没有人能够证明一个完整的收敛定理。因此,在极少数情况下,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])
要求计算环境影响评估(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}\)。精确特征值接近于一对双根。威尔金森双移位迭代使用每对中的一个特征值。此迭代确实会改变矩阵,但不足以获得快速收敛。因此,我们必须使用一个不同的双移位,基于重复较低的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算法计算特征值和奇异值是密切相关的,但其结果的应用通常是非常不同的。特征值经常被用来分析常微分方程组,其中行为作为时间的函数是重要的。另一方面,奇异值对于分析联立线性方程组的静态系统是有用的,其中方程组的数目往往不等于未知数的数目。
控制理论和控制设计自动化大量使用特征值。研究了控制理论中微分方程组的经典状态空间系统
\[\begin{align}\dot{x}&=Ax+Bu\\y&=Cx+Du\end{align}\]
利用QR算法计算\(A\)的特征值对于研究稳定性和可控性等问题至关重要。
在统计学中,奇异值分解是获得主成分的可靠方法。主成分分析(PCA)是一种分析超定线性方程组的技术
\[Ax=b\]
其中\(A\)的行多于列。利用QR算法计算\(A\)的奇异值和向量,得到主分量。
为进一步阅读和查看
Golub, Gene H.和Charles F. Van Loan,矩阵计算,第4版,约翰霍普金斯大学出版社,1996年,697页。
硅藻土、克里夫用MATLAB进行数值计算,第10章,“特征值和奇异值。”
硅藻土、克里夫1976矩阵奇异值分解膜.
威尔金森,J.H。,代数特征值问题,牛津大学出版社,1965年,662页。
注意:本文的一部分是基于“QR算法,”MATLAB新闻和注释1995年夏天,。