广义最小残量方法
在数学上,广义最小残量方法(一般简称GMRES)是一个求解线性方程组 数值解的迭代方法。这个方法利用在Krylov子空间中有着最小残量的向量来逼近解。Arnoldi迭代方法被用来求解这个向量。
GMRES方法由Yousef Saad和Martin H. Schultz在1986年提出。[1]
GMRES方法 编辑
需要求解的线性方程组记为
- 。
假设矩阵A是n n阶的可逆的。进一步,假设b是标准化的,即||b|| = 1 (在这篇文章中,||·||是Euclidean范数)。
这个问题的m阶Krylov子空间是
- 。
GMRES通过使得残量Axm − b最小的向量xm ∈ Km来逼近Ax = b的精确解。
但是,向量b, Ab, …, Am−1b几乎是线性相关的。因此,用Arnoldi迭代方法求得的这组Km的标准正交基
来取代上面的那组基。所以,向量xm ∈ Km写成xm = Qmym,其中ym ∈ Rm且Qm是由q1, …, qm组成的n m矩阵。
Arnoldi过程也产生一个 (m+1) m阶上Hessenberg矩阵 满足
- 。
因为 是正交的,我们有
- ,
其中
是Rm+1的标准基的第一个向量,并且
- ,
其中 是初始向量(通常是零向量)。因此,求使得残量
- 。
的范数最小的 。这是一个m阶线性最小二乘问题。
这就是GMRES方法。在迭代的每一步中:
- 做一步Arnoldi迭代方法;
- 寻找使得||rm||最小的 ;
- 计算 ;
- 如果残量不够小,重复以上过程。
在每一步迭代中,必须计算一次矩阵向量积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是正定的,则
- ,
其中 和 分别为矩阵 的最小和最大特征值。
如果A是对称的并且是正定的,则
- 。
其中 记为A在Euclidean范数下的条件数。
一般情况下,其中A是非正定的,则
- ,
其中Pm记为次数不超过m且p(0) = 1的多项式的集合,V是A的谱分解中的矩阵,而σ(A)是A的谱。粗略的说,当A的特征值聚集在远离原点的区域且A离正规不太远时,收敛速度较快。[2]
所有的不等式只界定残量,而不是实际误差(精确解和当前迭代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方法的其中一部分是求解向量 使得
最小。这个可以通过计算QR分解来实现:找到一个(m+1) (m+1)阶正交矩阵Ωm和一个(m+1) m上三角矩阵 满足
- 。
三角矩阵的行数比列数多1,所以它的最后一行由零组成。因此,它能被分解为
- ,
其中 是一个m m阶三角(方)矩阵。
QR分解能够简单的进行下去(update),从一步迭代到下一步迭代。因为每次的Hessenberg矩阵只在一行零元素和一列元素上有所不同:
- ,
其中hm = (h1m, … hmm)T。这意味着,Hessenberg矩阵左乘上Ωm的扩大矩阵(通过并上零元素和单位元素),所得到的是类似于三角矩阵的矩阵:
这个矩阵可以三角化,如果σ为零。为了修正这个矩阵,需要进行Givens旋转
其中
- 。
通过这个Givens旋转,我们构造
- 。
事实上,
是一个三角矩阵。
给出了QR分解,最小值问题就容易解决了。注意到
- 。
记 为
其中gm ∈ Rm和γm ∈ R,则
- 。
使得这个表达式最小的向量y为
- 。
再一次,向量 能够简单的进行下去(update)。[3]
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 <= 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()
注记 编辑
参考 编辑
- 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 (页面存档备份,存于互联网档案馆)