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

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

GMRES方法
需要求解的线性方程组记为
: Ax = b \, 。
假设矩阵A是n\timesn阶的可逆的。进一步，假设b是标准化的，即||b|| = 1 (在这篇文章中，||·||是Euclidean范数)。

这个问题的m阶Krylov子空间是
: K_m = \operatorname{span} \, \{ b, Ab, A^2b, \ldots, A^{m-1}b \} \, 。
GMRES通过使得残量Axm − b最小的向量xm ∈ Km来逼近Ax = b的精确解。

但是，向量b, Ab, …, Am−1b几乎是线性相关的。因此，用Arnoldi迭代方法求得的这组Km的标准正交基
: q_1, q_2, \ldots, q_m \,
来取代上面的那组基。所以，向量xm ∈ Km写成xm = Qmym，其中ym ∈ Rm且Qm是由q1, …, qm组成的n\timesm矩阵。

Arnoldi过程也产生一个 (m+1)\timesm阶上Hessenberg矩阵\tilde{H}_m满足
: AQ_m = Q_{m+1} \tilde{H}_m \, 。
因为Q_m是正交的，我们有
: \| Ax_m - b \| = \| \tilde{H}_my_m - \beta e_1 \| \, ，
其中
: e_1 = (1,0,0,\ldots,0) \,
是Rm+1的标准基的第一个向量，并且
: \beta = \|b-Ax_0\| \, ，
其中x_0是初始向量(通常是零向量)。因此，求使得残量
: r_m = \tilde{H}_m y_m - \beta e_1 。
的范数最小的x_m。这是一个m阶线性最小二乘问题。

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

做一步Arnoldi迭代方法；

寻找使得||rm||最小的 y_m ；

计算 x_m = Q_m y_m ；

如果残量不够小，重复以上过程。

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

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

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

在实验中，GMRES方法经常表现得很好。在特殊的情况下这能够被证明。如果A是正定的，则
: \|r_m\| \leq \left( 1-\frac{\lambda_{\mathrm{min}}(A^T + A)}{2 \lambda_{\mathrm{max}}(A^T + A)} \right)^{m/2} \|r_0\| ，
其中\lambda_{\mathrm{min}}(M)和\lambda_{\mathrm{max}}(M)分别为矩阵M的最小和最大特征值。

如果A是对称的并且是正定的，则
: \|r_m\| \leq \left( \frac{\kappa_2^2(A)-1}{\kappa_2^2(A)} \right)^{m/2} \|r_0\| 。
其中\kappa_2(A)记为A在Euclidean范数下的条件数。

一般情况下，其中A是非正定的，则
: \|r_m\| \le \inf_{p \in P_m} \|p_m(A)\| \le \kappa_2(V) \inf_{p \in P_m} \max_{\lambda \in \sigma(A)} |p(\lambda)| \, ，
其中Pm记为次数不超过m且p(0) = 1的多项式的集合，V是A的谱分解中的矩阵，而σ(A)是A的谱。粗略的说，当A的特征值聚集在远离原点的区域且A离正规不太远时，收敛速度较快。

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

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

迭代的开销以O(m2)增长，其中m是迭代次数。然而有时候，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，但他们没有达到最小的残量，因此对于这些方法残量不会单调递减。收敛性是不能保证的。

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

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

求解最小二乘问题
GMRES方法的其中一部分是求解向量y_m使得
: \| \tilde{H}_m y_m - e_1 \| \,
最小。这个可以通过计算QR分解来实现：找到一个(m+1)\times(m+1)阶正交矩阵Ωm和一个(m+1)\timesm上三角矩阵\tilde{R}_m满足
: \Omega_m \tilde{H}_m = \tilde{R}_m 。
三角矩阵的行数比列数多1，所以它的最后一行由零组成。因此，它能被分解为
: \tilde{R}_m = \begin{bmatrix} R_m \\ 0 \end{bmatrix} ，
其中R_m是一个m\timesm阶三角(方)矩阵。

QR分解能够简单的进行下去(update)，从一步迭代到下一步迭代。因为每次的Hessenberg矩阵只在一行零元素和一列元素上有所不同：
:\tilde{H}_{m+1} = \begin{bmatrix} \tilde{H}_m & h_m \\ 0 & h_{m+1,m} \end{bmatrix} ，
其中hm = (h1m, … hmm)T。这意味着，Hessenberg矩阵左乘上Ωm的扩大矩阵(通过并上零元素和单位元素)，所得到的是类似于三角矩阵的矩阵：
: \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旋转
: G_m = \begin{bmatrix} I_{m-1} & 0 & 0 \\ 0 & c_m & s_m \\ 0 & -s_m & c_m \end{bmatrix}
其中
: c_m = \frac{\rho}{\sqrt{\rho^2+\sigma^2}} \quad\mbox{and}\quad s_m = \frac{\sigma}{\sqrt{\rho^2+\sigma^2}} 。
通过这个Givens旋转，我们构造
: \Omega_{m+1} = G_m \begin{bmatrix} \Omega_m & 0 \\ 0 & 1 \end{bmatrix} 。
事实上，
: \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分解，最小值问题就容易解决了。注意到
: \| \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 \| 。
记\Omega_me_1为
: \tilde{g}_m = \begin{bmatrix} g_m \\ \gamma_m \end{bmatrix}
其中gm ∈ Rm和γm ∈ R，则
: \| \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为
: y_m = R_m^{-1} g_m 。
再一次，向量g_m能够简单的进行下去(update)。

Example code
Regular GMRES (python3)

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 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()

注记
参考

• 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. .
• 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.
• [http://www.netlib.org/linalg/html_templates/node29.html#SECTION00734000000000000000|J. 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