廣義最小殘量方法

在數學上，廣義最小殘量方法(一般簡稱GMRES)是一個求解線性方程組 數值解的迭代方法。這個方法利用在Krylov子空間中有著最小殘量的向量來逼近解。Arnoldi迭代方法被用來求解這個向量。GMRES方法由Yousef Saad和Martin H. Schultz在1986年提出。GMRES方法.需要求解的線性方程組記為formula_1。假設矩陣"A"是"n"formula_2"n"階的可逆的。進一步，假設"b"是標準化的，即||"b"|| = 1 (在這篇文章中，||·||是Euclidean範數)。這個問題的"m"階Krylov子空間是formula_3。GMRES透過使得殘量"Ax""m" − "b"最小的向量"x""m" ∈ "K""m"來逼近"Ax" = "b"的精確解。但是，向量"b", "Ab", …, "A""m"−1"b"幾乎是線性相關的。因此，用Arnoldi迭代方法求得的這組"K""m"的標準正交基formula_4來取代上面的那組基。所以，向量"x""m" ∈ "K""m"寫成"x""m" = "Q""m""y""m"，其中"y""m" ∈ R"m"且"Q""m"是由"q"1, …, "q"m組成的"n"formula_2"m"矩陣。Arnoldi過程也產生一個 ("m"+1)formula_2"m"階上Hessenberg矩陣formula_7滿足formula_8。因為formula_9是正交的，我們有formula_10，其中formula_11是R"m"+1的標準基的第一個向量，並且formula_12 ，其中formula_13是初始向量(通常是零向量)。因此，求使得殘量formula_14。的範數最小的formula_15。這是一個"m"階線性最小二乘問題。這就是GMRES方法。在迭代的每一步中：在每一步迭代中，必須計算一次矩陣向量積"Aq""m"。對於一般的"n"階稠密矩陣，這要計算複雜度大約2"n"2次浮點運算。但是對於稀疏矩陣，這個計算複雜度能減少到O("n")。進一步，關於矩陣向量積，在"m"次迭代中能進行O("m" "n")次浮點運算。收斂性.第"m"次迭代獲得在Krylov子空間"K""m"下的最小殘量。因為每個子空間包含於下一個子空間中，所以殘量單調遞減。在第"n"次迭代後，其中"n"是矩陣"A"的階數，Krylov空間"K""n"是完整的R"n"。因此，GMRES方法達到精確解。然而，問題在於：在極少的幾次迭代後(相對於"n")，向量"x""m"幾乎已經是精確解的一個很好的逼近。但是，在一般情況下這是不會發生的。事實上，Greenbaum，Pták和Strakoš的理論說明了對於每一個單調減少的序列"a"1, …, "a""n"−1, "a""n" = 0 ，能夠找到一個矩陣"A"對於所有"m"滿足||"r""m"|| = "a""m" ，其中"r""m"是上面所定義的殘量。特別的，有可能找到一個矩陣，使得前"n" − 1次迭代的殘量一直保持為常數，而只在最後一次迭代時達到零。在實驗中，GMRES方法經常表現得很好。在特殊的情況下這能夠被證明。如果"A"是正定的，則formula_18，其中formula_19和formula_20分別為矩陣formula_21的最小和最大特徵值。如果"A"是對稱的並且是正定的，則formula_22。其中formula_23記為"A"在Euclidean範數下的條件數。一般情況下，其中"A"是非正定的，則formula_24，其中"P""m"記為次數不超過"m"且"p"(0) = 1的多項式的集合，"V"是"A"的譜分解中的矩陣，而σ("A")是"A"的譜。粗略的說，當"A"的特徵值聚集在遠離原點的區域且"A"離正規不太遠時，收斂速度較快。所有的不等式只界定殘量，而不是實際誤差(精確解和當前迭代"x""m"之間的距離)。GMRES方法的拓展( Restarted GMRES ).同其他迭代方法一樣，為了加快收斂，GMRES經常結合預處理方法。迭代的開銷以O("m"2)增長，其中"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方法的其中一部分是求解向量formula_25使得formula_26最小。這個可以透過計算QR分解來實現：找到一個("m"+1)formula_2("m"+1)階正交矩陣Ω"m"和一個("m"+1)formula_2"m"上三角矩陣formula_29滿足formula_30。三角矩陣的行數比列數多1，所以它的最後一行由零組成。因此，它能被分解為formula_31，其中formula_32是一個"m"formula_2"m"階三角(方)矩陣。QR分解能夠簡單的進行下去(update)，從一步迭代到下一步迭代。因為每次的Hessenberg矩陣只在一行零元素和一列元素上有所不同：formula_34 ，其中"h""m" = ("h"1"m", … "h""mm")T。這意味著，Hessenberg矩陣左乘上Ω"m"的擴大矩陣(透過並上零元素和單位元素)，所得到的是類似於三角矩陣的矩陣：formula_35這個矩陣可以三角化，如果σ為零。為了修正這個矩陣，需要進行Givens旋轉formula_36其中formula_37。透過這個Givens旋轉，我們構造formula_38。事實上，formula_39是一個三角矩陣。給出了QR分解，最小值問題就容易解決了。注意到formula_40。記formula_41為formula_42其中"g""m" ∈ R"m"和γ"m" ∈ R，則formula_43。使得這個表示式最小的向量"y"為formula_44。再一次，向量formula_45能夠簡單的進行下去(update)。Example code.Regular GMRES (python3).import numpy as npimport mathclass 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_matrixself.b = b_boundary_condition_vectorself.maximum_number_of_basis_used = maximum_number_of_basis_usedself.threshold = thresholddef initial_guess_input( self, x_input_vector_initial_guess: np.array([], dtype = float ) ):self.x = x_input_vector_initial_guesstry: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_usedr = 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_normself.e = [self.error]# initialize the 1D vectorssn = np.zeros( m )cs = np.zeros( m )e1 = np.zeros( m + 1 )e1[0] = 1.0beta = r_norm * e1# beta is the beta vector instead of the beta scalarH = np.zeros(( m+1, m+1 ))Q = np.zeros(( n, m+1 ))Q[:,0] = r / r_normfor 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 vectorbeta[ k+1 ] = -sn[k] * beta[k]beta[ k ] = cs[k] * beta[k]# calculate and save the errorsself.error = abs(beta[k+1]) / b_normself.e = np.append(self.e, self.error)if( self.error  np.ndarray:n = b.sizeif A[n-1, n-1] == 0.0:raise ValueErrorx = np.zeros_like(b)x[n-1] = b[n-1] / A[n-1, n-1]for i in range( n-2, -1, -1 ):bb = 0for j in range ( i+1, n ):bb += A[i, j] * x[j]x[i] = (b[i] - bb) / A[i, i]return xdef 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 = 100for 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()