廣義最小殘量方法

在數學上，廣義最小殘量方法(一般簡稱GMRES)是一個求解線性方程組數值解的迭代方法。這個方法利用在Krylov子空間中有着最小殘量的向量來逼近解。Arnoldi迭代方法被用來求解這個向量。

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

GMRES方法

需要求解的線性方程組記為

Ax=b\,

。

假設矩陣A是n $\times$ n階的可逆的。進一步，假設b是標準化的，即||b|| = 1 (在這篇文章中，||·||是Euclidean範數)。

這個問題的m階Krylov子空間是

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的標準正交基

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 $\times$ m矩陣。

Arnoldi過程也產生一個 (m+1) $\times$ m階上Hessenberg矩陣 ${\tilde {H}}_{m}$ 滿足

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

。

因為 $Q_{m}$ 是正交的，我們有

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

，

其中

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

是R^m+1的標準基的第一個向量，並且

\beta =\|b-Ax_{0}\|\,

，

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

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

。

的範數最小的 $x_{m}$ 。這是一個m階線性最小二乘問題。

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

做一步Arnoldi迭代方法；
尋找使得||r_m||最小的 $y_{m}$ ；
計算 $x_{m}=Q_{m}y_{m}$ ；
如果殘量不夠小，重複以上過程。

在每一步迭代中，必須計算一次矩陣向量積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是正定的，則

\|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}\|\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方法的拓展( Restarted GMRES )

同其他迭代方法一樣，為了加快收斂，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方法的其中一部分是求解向量 $y_{m}$ 使得

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

最小。這個可以通過計算QR分解來實現：找到一個(m+1) $\times$ (m+1)階正交矩陣Ω_m和一個(m+1) $\times$ m上三角矩陣 ${\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 $\times$ m階三角(方)矩陣。

QR分解能夠簡單的進行下去(update)，從一步迭代到下一步迭代。因為每次的Hessenberg矩陣只在一行零元素和一列元素上有所不同：

{\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的擴大矩陣(通過並上零元素和單位元素)，所得到的是類似於三角矩陣的矩陣：

{\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 _{m}e_{1}$ 為

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

其中g_m ∈ R^m和γ_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)。^[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()

註記

^ Saad和Schultz
^ Trefethen & Bau, Thm 35.2
^ 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 （頁面存檔備份，存於網際網路檔案館）

[1] Saad和Schultz

[2] Trefethen & Bau, Thm 35.2

[3] Stoer and Bulirsch, §8.7.2

[1]

[2]

[3]