广义最小残量方法

维基百科,自由的百科全书
跳转至: 导航搜索

在数学上,广义最小残量方法(一般简称GMRES)是一个求解线性方程组 数值解的迭代方法。这个方法利用在Krylov子空间中有着最小残量的向量来逼近解。Arnoldi迭代方法被用来求解这个向量。

GMRES方法由Yousef Saad和Martin H. Schultz在1986年提出。[1]

GMRES方法[编辑]

需要求解的线性方程组记为

假设矩阵Amm阶的可逆的。进一步,假设b是标准化的,即||b|| = 1 (在这篇文章中,||·||是Euclidean范数)。

这个问题的nKrylov子空间

GMRES通过使得残量Axnb最小的向量xnKn来逼近Ax = b的精确解。

但是,向量b, Ab, …, An−1b几乎是线性相关的。因此,用Arnoldi迭代方法求得的这组Kn的标准正交基

来取代上面的那组基。所以,向量xnKn写成xn = Qnyn,其中ynRnQn是由q1, …, qn组成的mn矩阵。

Arnoldi过程也产生一个 (n+1)n阶上Hessenberg矩阵满足

因为是正交的,我们有

其中

Rn+1的标准的第一个向量,并且

其中是初始向量(通常是零向量)。因此,求使得残量

的范数最小的。这是一个n阶线性最小二乘问题

这就是GMRES方法。在迭代的每一步中:

  1. 做一步Arnoldi迭代方法;
  2. 寻找使得||rn||最小的
  3. 计算
  4. 如果残量不够小,重复以上过程。

在每一步迭代中,必须计算一次矩阵向量积Aqn。对于一般的m阶稠密矩阵,这要计算复杂度大约2m2浮点运算。但是对于稀疏矩阵,这个计算复杂度能减少到O(m)。进一步,关于矩阵向量积,在n次迭代中能进行O(n m)次浮点运算。

收敛性[编辑]

n次迭代获得在Krylov子空间Kn下的最小残量。因为每个子空间包含于下一个子空间中,所以残量单调递减。在第m次迭代后,其中m是矩阵A的阶数,Krylov空间Km是完整的Rm。因此,GMRES方法达到精确解。然而,问题在于:在极少的几次迭代后(相对于m),向量xn几乎已经是精确解的一个很好的逼近。

但是,在一般情况下这是不会发生的。事实上,Greenbaum,Pták和Strakoš的理论说明了对于每一个单调减少的序列a1, …, am−1, am = 0 ,能够找到一个矩阵A对于所有n满足||rn|| = an ,其中rn是上面所定义的残量。特别的,有可能找到一个矩阵,使得前m − 1次迭代的残量一直保持为常数,而只在最后一次迭代时达到零。

在实验中,GMRES方法经常表现得很好。在特殊的情况下这能够被证明。如果A正定的,则

其中分别为矩阵的最小和最大特征值

如果A对称的并且是正定的,则

其中记为A在Euclidean范数下的条件数

一般情况下,其中A是非正定的,则

其中Pn记为次数不超过np(0) = 1的多项式的集合,VA谱分解中的矩阵,而σ(A)是A。粗略的说,当A的特征值聚集在远离原点的区域且A正规不太远时,收敛速度较快。[2]

所有的不等式只界定残量,而不是实际误差(精确解和当前迭代xn之间的距离)。

GMRES方法的拓展[编辑]

同其他迭代方法一样,为了加快收敛,GMRES经常结合预处理方法。

迭代的开销以O(n2)增长,其中n是迭代次数。然而有时候,GMRES方法在k次迭代后重新开始,即xk又变回初始值。这样的方法叫做GMRES(k)。

与其他解法的比较[编辑]

对于对称矩阵,Arnoldi迭代方法变成Lanczos迭代方法。对应的Krylov子空间方法叫做Paige和Saunders的最小残量方法(MinRes)。不像非对称的情况,MinRes方法由三项循环关系(three-term recurrence relation)给出,并且同GMRES一样,使残量的范数最小。而对于一般矩阵,Krylov子空间方法不能由短的循环关系(short recurrence relation)给出。

另一类方法由非对称Lanczos迭代方法给出,特别的是BiCG方法。这个利用了three-term recurrence relation,但他们没有达到最小的残量,因此对于这些方法残量不会单调递减。收敛性是不能保证的。

第三类方法由CGSBiCGSTAB给出。这些也由three-term recurrence relation给出(因此,非最优)。而且可能过早的终止迭代了而没有达到收敛的目的。这些方法的想法是合适的选择迭代序列所产生的多项式。

对于所有矩阵,这三类方法都不是最好的;总有例使得一类方法好于另一类。因而,各种解法应该进行实际的试验,来决定对于给定的问题哪一种是最优的。

求解最小二乘问题[编辑]

GMRES方法的其中一部分是求解向量使得

最小。这个可以通过计算QR分解来实现:找到一个(n+1)(n+1)阶正交矩阵Ωn和一个(n+1)n三角矩阵满足

三角矩阵的行数比列数多1,所以它的最后一行由零组成。因此,它能被分解为

其中是一个nn阶三角(方)矩阵。

QR分解能够简单的进行下去(update),从一步迭代到下一步迭代。因为每次的Hessenberg矩阵只在一行零元素和一列元素上有所不同:

其中hn = (h1n, … hnn)T。这意味着,Hessenberg矩阵左乘上Ωn的扩大矩阵(通过并上零元素和单位元素),所得到的是类似于三角矩阵的矩阵:

这个矩阵可以三角化,如果σ为零。为了修正这个矩阵,需要进行Givens旋转

其中

通过这个Givens旋转,我们构造

事实上,

是一个三角矩阵。

给出了QR分解,最小值问题就容易解决了。注意到

其中gnRn和γnR,则

使得这个表达式最小的向量y

再一次,向量能够简单的进行下去(update)。[3]

注记[编辑]

  1. ^ Saad和Schultz
  2. ^ Trefethen & Bau, Thm 35.2
  3. ^ Stoer and Bulirsch, §8.7.2

参考[编辑]