广义最小残量方法

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

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

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

${\displaystyle Ax=b\,}$。

假设矩阵A是n${\displaystyle \times }$n阶的可逆的。进一步，假设b是标准化的，即||b|| = 1 (在这篇文章中，||·||是Euclidean范数)。

这个问题的m阶Krylov子空间是

${\displaystyle K_{m}=\operatorname {span} \,\{b,Ab,A^{2}b,\ldots ,A^{m-1}b\}\,}$。

GMRES通过使得残量Ax_m − b最小的向量x_m ∈ K_m来逼近Ax = b的精确解。

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

${\displaystyle q_{1},q_{2},\ldots ,q_{m}\,}$

来取代上面的那组基。所以，向量x_m ∈ K_m写成x_m = Q_my_m，其中y_m ∈ R^m且Q_m是由q₁, …, q_m组成的n${\displaystyle \times }$m矩阵。

Arnoldi过程也产生一个 (m+1)${\displaystyle \times }$m阶上Hessenberg矩阵${\displaystyle {\tilde {H}}_{m}}$满足

${\displaystyle AQ_{m}=Q_{m+1}{\tilde {H}}_{m}\,}$。

因为${\displaystyle Q_{m}}$是正交的，我们有

${\displaystyle \|Ax_{m}-b\|=\|{\tilde {H}}_{m}y_{m}-\beta e_{1}\|\,}$，

其中

${\displaystyle e_{1}=(1,0,0,\ldots ,0)\,}$

是R^m+1的标准基的第一个向量，并且

${\displaystyle \beta =\|b-Ax_{0}\|\,}$ ，

其中${\displaystyle x_{0}}$是初始向量(通常是零向量)。因此，求使得残量

${\displaystyle r_{m}={\tilde {H}}_{m}y_{m}-\beta e_{1}}$。

的范数最小的${\displaystyle x_{m}}$。这是一个m阶线性最小二乘问题。

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

1. 做一步Arnoldi迭代方法；
2. 寻找使得||r_m||最小的${\displaystyle y_{m}}$ ；
3. 计算${\displaystyle x_{m}=Q_{m}y_{m}}$ ；
4. 如果残量不够小，重复以上过程。

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

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

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

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

${\displaystyle \|r_{m}\|\leq \left(1-{\frac {\lambda _{\mathrm {min} }(A^{T}+A)}{2\lambda _{\mathrm {max} }(A^{T}+A)}}\right)^{m/2}\|r_{0}\|}$，

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

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

${\displaystyle \|r_{m}\|\leq \left({\frac {\kappa _{2}^{2}(A)-1}{\kappa _{2}^{2}(A)}}\right)^{m/2}\|r_{0}\|}$。

其中${\displaystyle \kappa _{2}(A)}$记为A在Euclidean范数下的条件数。

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

${\displaystyle \|r_{m}\|\leq \inf _{p\in P_{m}}\|p_{m}(A)\|\leq \kappa _{2}(V)\inf _{p\in P_{m}}\max _{\lambda \in \sigma (A)}|p(\lambda )|\,}$，

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

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

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

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

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

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

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

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

GMRES方法的其中一部分是求解向量${\displaystyle y_{m}}$使得

${\displaystyle \|{\tilde {H}}_{m}y_{m}-e_{1}\|\,}$

最小。这个可以通过计算QR分解来实现：找到一个(m+1)${\displaystyle \times }$(m+1)阶正交矩阵Ω_m和一个(m+1)${\displaystyle \times }$m上三角矩阵${\displaystyle {\tilde {R}}_{m}}$满足

${\displaystyle \Omega _{m}{\tilde {H}}_{m}={\tilde {R}}_{m}}$。

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

${\displaystyle {\tilde {R}}_{m}={\begin{bmatrix}R_{m}\\0\end{bmatrix}}}$，

其中${\displaystyle R_{m}}$是一个m${\displaystyle \times }$m阶三角(方)矩阵。

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

${\displaystyle {\tilde {H}}_{m+1}={\begin{bmatrix}{\tilde {H}}_{m}&h_{m}\\0&h_{m+1,m}\end{bmatrix}}}$ ，

其中h_m = (h_1m, … h_mm)^T。这意味着，Hessenberg矩阵左乘上Ω_m的扩大矩阵(通过并上零元素和单位元素)，所得到的是类似于三角矩阵的矩阵：

${\displaystyle {\begin{bmatrix}\Omega _{m}&0\\0&1\end{bmatrix}}{\tilde {H}}_{m+1}={\begin{bmatrix}R_{m}&r_{k}\\0&\rho \\0&\sigma \end{bmatrix}}}$

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

${\displaystyle G_{m}={\begin{bmatrix}I_{m-1}&0&0\\0&c_{m}&s_{m}\\0&-s_{m}&c_{m}\end{bmatrix}}}$

其中

${\displaystyle c_{m}={\frac {\rho }{\sqrt {\rho ^{2}+\sigma ^{2}}}}\quad {\mbox{and}}\quad s_{m}={\frac {\sigma }{\sqrt {\rho ^{2}+\sigma ^{2}}}}}$。

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

${\displaystyle \Omega _{m+1}=G_{m}{\begin{bmatrix}\Omega _{m}&0\\0&1\end{bmatrix}}}$。

事实上，

${\displaystyle \Omega _{m+1}{\tilde {H}}_{m+1}={\begin{bmatrix}R_{m}&r_{m}\\0&r_{mm}\\0&0\end{bmatrix}}\quad {\text{其 中}}\quad r_{mm}={\sqrt {\rho ^{2}+\sigma ^{2}}}}$

是一个三角矩阵。

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

${\displaystyle \|{\tilde {H}}_{m}y_{m}-e_{1}\|=\|\Omega _{m}({\tilde {H}}_{m}y_{m}-e_{1})\|=\|{\tilde {R}}_{m}y_{m}-\Omega _{m}e_{1}\|}$。

记${\displaystyle \Omega _{m}e_{1}}$为

${\displaystyle {\tilde {g}}_{m}={\begin{bmatrix}g_{m}\\\gamma _{m}\end{bmatrix}}}$

其中g_m ∈ R^m和γ_m ∈ R，则

${\displaystyle \|{\tilde {H}}_{m}y_{m}-e_{1}\|=\|{\tilde {R}}_{m}y_{m}-\Omega _{m}e_{1}\|=\left\|{\begin{bmatrix}R_{m}\\0\end{bmatrix}}y-{\begin{bmatrix}g_{m}\\\gamma _{m}\end{bmatrix}}\right\|}$。

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

${\displaystyle y_{m}=R_{m}^{-1}g_{m}}$。

再一次，向量${\displaystyle g_{m}}$能够简单的进行下去(update)。^[3]

# from "https://github.com/J-N-ch/GMRES_py_restart/blob/master/GMRES_API/GMRES.py" import numpy as np import math  class GMRES_API(object):     def __init__( self,                   A_coefficient_matrix: np.array([], dtype = float ),                   b_boundary_condition_vector: np.array([], dtype = float ),                   maximum_number_of_basis_used: int,                   threshold = 1.0e-16 ):          self.A = A_coefficient_matrix         self.b = b_boundary_condition_vector         self.maximum_number_of_basis_used = maximum_number_of_basis_used         self.threshold = threshold      def initial_guess_input( self, x_input_vector_initial_guess: np.array([], dtype = float ) ):          self.x = x_input_vector_initial_guess          try:             assert len( self.x ) == len( self.b )          except Exception:              print(" The input guess vector's size must equal to the system's size !\n")             print(" The matrix system's size == ", len( self.b ))             print(" Your input vector's size == ", len( self.x ))             self.x = np.zeros( len( self.b ) )              print(" Use default input guess vector = ", self.x, " instead of the incorrect vector you given !\n")       def run( self ):          n = len( self.A )         m = self.maximum_number_of_basis_used          r = self.b - np.dot(self.A , self.x)         r_norm = np.linalg.norm( r )          b_norm = np.linalg.norm( self.b )          self.error = np.linalg.norm( r ) / b_norm         self.e = [self.error]                  # initialize the 1D vectors          sn = np.zeros( m )         cs = np.zeros( m )         e1 = np.zeros( m + 1 )         e1[0] = 1.0          beta = r_norm * e1          # beta is the beta vector instead of the beta scalar          H = np.zeros(( m+1, m+1 ))         Q = np.zeros((   n, m+1 ))         Q[:,0] = r / r_norm          for k in range(m):              ( H[0:k+2, k], Q[:, k+1] )    = __class__.arnoldi( self.A, Q, k)             ( H[0:k+2, k], cs[k], sn[k] ) = __class__.apply_givens_rotation( H[0:k+2, k], cs, sn, k)                          # update the residual vector             beta[ k+1 ] = -sn[k] * beta[k]             beta[ k   ] =  cs[k] * beta[k]              # calculate and save the errors             self.error = abs(beta[k+1]) / b_norm             self.e = np.append(self.e, self.error)              if( self.error <= self.threshold):                 break           # calculate the result         #y = np.matmul( np.linalg.inv( H[0:k+1, 0:k+1]), beta[0:k+1] )         #TODO Due to H[0:k+1, 0:k+1] being a upper tri-matrix, we can exploit this fact.          y = __class__.__back_substitution( H[0:k+1, 0:k+1], beta[0:k+1] )           self.x = self.x + np.matmul( Q[:,0:k+1], y )          self.final_residual_norm = np.linalg.norm( self.b - np.matmul( self.A, self.x ) )          return self.x       '''''''''''''''''''''''''''''''''''     '        Arnoldi Function         '     '''''''''''''''''''''''''''''''''''     @staticmethod     def arnoldi( A, Q, k ):         h = np.zeros( k+2 )         q = np.dot( A, Q[:,k] )         for i in range ( k+1 ):             h[i] = np.dot( q, Q[:,i])             q = q - h[i] * Q[:, i]         h[ k+1 ] = np.linalg.norm(q)         q = q / h[ k+1 ]         return h, q       '''''''''''''''''''''''''''''''''''''''''''''''''''''''''     '           Applying Givens Rotation to H col           '     '''''''''''''''''''''''''''''''''''''''''''''''''''''''''     @staticmethod     def apply_givens_rotation( h, cs, sn, k ):         for i in range( k-1 ):             temp   =  cs[i] * h[i] + sn[i] * h[i+1]             h[i+1] = -sn[i] * h[i] + cs[i] * h[i+1]             h[i]   = temp          # update the next sin cos values for rotation         cs_k, sn_k, h[k] = __class__.givens_rotation( h[k-1], h[k] )                  # eliminate H[ k+1, i ]         h[k + 1] = 0.0          return h, cs_k, sn_k      ##----Calculate the Given rotation matrix----##     # From "http://www.netlib.org/lapack/lawnspdf/lawn150.pdf"     # The algorithm used by "Edward Anderson"     @staticmethod     def givens_rotation( v1, v2 ):         if( v2 == 0.0 ):             cs = np.sign(v1)             sn = 0.0             r = abs(v1)         elif( v1 == 0.0 ):             cs = 0.0             sn = np.sign(v2)             r = abs(v2)         elif( abs(v1) > abs(v2) ):             t = v2 / v1              u = np.sign(v1) * math.hypot( 1.0, t )               cs = 1.0 / u             sn = t * cs             r = v1 * u         else:             t = v1 / v2              u = np.sign(v2) * math.hypot( 1.0, t )               sn = 1.0 / u             cs = t * sn             r = v2 * u         return cs, sn, r      # From https://stackoverflow.com/questions/47551069/back-substitution-in-python     @staticmethod     def __back_substitution( A: np.ndarray, b: np.ndarray) -> np.ndarray:         n = b.size         if A[n-1, n-1] == 0.0:             raise ValueError          x = np.zeros_like(b)         x[n-1] = b[n-1] / A[n-1, n-1]         for i in range( n-2, -1, -1 ):             bb = 0             for j in range ( i+1, n ):                 bb += A[i, j] * x[j]             x[i] = (b[i] - bb) / A[i, i]         return x       def final_residual_info_show( self ):         print( "x  =", self.x, "residual_norm =  ", self.final_residual_norm )       def main():      A_mat = np.array( [[1.00, 1.00, 1.00],                        [1.00, 2.00, 1.00],                        [0.00, 0.00, 3.00]] )      b_mat = np.array( [3.0, 2.0, 1.0] )      GMRES_test_itr2 = GMRES_API( A_mat, b_mat, 2, 0.01)      x_mat = np.array( [1.0, 1.0, 1.0] )     print("x  =", x_mat)      # GMRES with restart, 2 iterations in each restart ( GMRES(2) )     max_restart_counts = 100     for restart_counter in range(max_restart_counts):         GMRES_test_itr2.initial_guess_input( x_mat )          x_mat = GMRES_test_itr2.run()         print(restart_counter+1," : x  =", x_mat)      xx = np.matmul( np.linalg.inv(A_mat), b_mat )     print("ANS : xx =", xx)    if __name__ == '__main__':     main() 

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

• A. Meister, Numerik linearer Gleichungssysteme, 2nd edition, Vieweg 2005, ISBN 978-3-528-13135-7.
• Y. Saad, Iterative Methods for Sparse Linear Systems, 2nd edition, Society for Industrial and Applied Mathematics, 2003. ISBN 978-0-89871-534-7.
• Y. Saad and M.H. Schultz, "GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems", SIAM J. Sci. Stat. Comput., 7:856-869, 1986. doi:10.1137/0907058.
• J. Stoer and R. Bulirsch, Introduction to numerical analysis, 3rd edition, Springer, New York, 2002. ISBN 978-0-387-95452-3.
• Lloyd N. Trefethen and David Bau, III, Numerical Linear Algebra, Society for Industrial and Applied Mathematics, 1997. ISBN 978-0-89871-361-9.
• Dongarra et al. , Templates for the Solution of Linear Systems: Building Blocks for Iterative Methods （页面存档备份，存于互联网档案馆）, 2nd Edition, SIAM, Philadelphia, 1994
• https://github.com/J-N-ch/GMRES_py_restart （页面存档备份，存于互联网档案馆）