广义最小残量方法

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

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

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

GMRES方法[编辑]

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

 Ax = b \,

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

这个问题的nKrylov子空间

 K_n = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{n-1}b \} \,

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

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

 q_1, q_2, \ldots, q_n \,

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

Arnoldi过程也产生一个 (n+1)\timesn阶上Hessenberg矩阵\tilde{H}_n满足

 AQ_n = Q_{n+1} \tilde{H}_n \,

因为Q_n是正交的,我们有

 \| Ax_n - b \| = \| \tilde{H}_ny_n - \beta e_1 \| \,

其中

 e_1 = (1,0,0,\ldots,0) \,

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

 \beta = \|b-Ax_0\| \,

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

 r_n = \tilde{H}_n y_n - \beta e_1

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

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

  1. 做一步Arnoldi迭代方法;
  2. 寻找使得||rn||最小的 y_n
  3. 计算 x_n = Q_n y_n
  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正定的,则

 \|r_n\| \leq \left( 1-\frac{\lambda_{\mathrm{min}}(A^T + A)}{2 \lambda_{\mathrm{max}}(A^T + A)} \right)^{n/2} \|r_0\|

其中\lambda_{\mathrm{min}}(M)\lambda_{\mathrm{max}}(M)分别为矩阵M的最小和最大特征值

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

 \|r_n\| \leq \left( \frac{\kappa_2^2(A)-1}{\kappa_2^2(A)} \right)^{n/2} \|r_0\|

其中\kappa_2(A)记为A在Euclidean范数下的条件数

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

 \|r_n\| \le \inf_{p \in P_n} \|p_n(A)\| \le \kappa_2(V) \inf_{p \in P_n} \max_{\lambda \in \sigma(A)} |p(\lambda)| \,

其中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方法的其中一部分是求解向量y_n使得

 \| \tilde{H}_n y_n - e_1 \| \,

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

 \Omega_n \tilde{H}_n = \tilde{R}_n

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

 \tilde{R}_n = \begin{bmatrix} R_n \\ 0 \end{bmatrix}

其中R_n是一个n\timesn阶三角(方)矩阵。

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

\tilde{H}_{n+1} = \begin{bmatrix} \tilde{H}_n & h_n \\ 0 & h_{n+1,n} \end{bmatrix}

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

 \begin{bmatrix} \Omega_n & 0 \\ 0 & 1 \end{bmatrix} \tilde{H}_{n+1} = \begin{bmatrix} R_n & r_k \\ 0 & \rho \\ 0 & \sigma \end{bmatrix}

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

 G_n = \begin{bmatrix} I_{n-1} & 0 & 0 \\ 0 & c_n & s_n \\ 0 & -s_n & c_n \end{bmatrix}

其中

 c_n = \frac{\rho}{\sqrt{\rho^2+\sigma^2}} \quad\mbox{and}\quad s_n = \frac{\sigma}{\sqrt{\rho^2+\sigma^2}}

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

 \Omega_{n+1} = G_n \begin{bmatrix} \Omega_n & 0 \\ 0 & 1 \end{bmatrix}

事实上,

 \Omega_{n+1} \tilde{H}_{n+1} = \begin{bmatrix} R_n & r_n \\ 0 & r_{nn} \\ 0 & 0 \end{bmatrix} \quad\text{其中}\quad r_{nn} = \sqrt{\rho^2+\sigma^2}

是一个三角矩阵。

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

 \| \tilde{H}_n y_n - e_1 \| = \| \Omega_n (\tilde{H}_n y_n - e_1) \| = \| \tilde{R}_n y_n - \Omega_n e_1 \|

\Omega_ne_1

 \tilde{g}_n = \begin{bmatrix} g_n \\ \gamma_n \end{bmatrix}

其中gnRn和γnR,则

 \| \tilde{H}_n y_n - e_1 \| = \| \tilde{R}_n y_n - \Omega_n e_1 \| = \left\| \begin{bmatrix} R_n \\ 0 \end{bmatrix} y - \begin{bmatrix} g_n \\ \gamma_n \end{bmatrix} \right\|

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

 y_n = R_n^{-1} g_n

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

注记[编辑]

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

参考[编辑]