在之前的文章【數(shù)值算法】共軛梯度法求解線性方程組中，我們指出，共軛梯度法是求解對稱正定系數(shù)的線性方程組的極為有效的方法，并且指出：對于n階線性方程組，通常最多n+1次迭代可以獲得收斂。在實際的有限元開發(fā)實踐中，需要求解的經常是非常大，極端大的線性方程組，此時，采用共軛梯度法，有時可能需要上萬次的迭代才能收斂，雖然每次收斂的計算量并不大，但是整體求解也會花費較多時間。

實際上，共軛梯度法能否快速收斂與系數(shù)矩陣A密切相關。文獻【1】指出，對于共軛梯度法存在以下定理：

如果A=I+B為nxn的對稱正定矩陣，且rank（B）=r，則共軛梯度法至多r+1步收斂。其中I指單位矩陣，rank是矩陣的秩。

從該定理可知，共軛梯度法能否快速收斂，主要取決于“A”是否“接近于”單位矩陣。特別地，當A就是單位矩陣時，rank(B)=0，此時一步即可收斂，當然這是很顯然的。

基于上述定理，多種基于共軛梯度法法衍生的預處理共軛梯度法（PCG）得以廣泛地應用。具體來說，實際過程如下：

假設要求解的方程是

預處理共軛梯度法的思想是求解

其中，

并且，C為對稱正定矩陣，且使得 是良態(tài)（即盡量“接近”于單位矩陣）。如果定義

結合原始的共軛梯度法，可以得到以下的預處理共軛梯度法（PCG）：

從上述算法中可以看出，其主要改動是需要首先求解

這一預處理方程，得到Zk的值后再代替原來的rk進行運算。采用不同的M，當然就會得到不同的Zk，最終影響共軛梯度法的收斂性。很顯然，對于整個原始方程的PCG求解效率來說，首先這個預處理方程需要比較容易求解，否則求解該預處理方程就已經花費了較多時間，則整個原始方程求解效率自然不會高。其次，M的選取需要確實能夠減少迭代步數(shù)。這兩方面在實際中往往是沖突的，例如，假設M取為簡單的單位矩陣，則Zk就是rk，幾乎不需要求解，此時PCG共軛梯度法就退化為普通的共軛梯度法，預處理相當于未生效；假設M取原始系數(shù)矩陣A，則很顯然，一步迭代就能得到最終解，但是求解該預處理方程，本質上就和求解原始方程無二。

在實踐中，尤其是有限元程序開發(fā)中，常見的M的選取至少有兩種：對角線預處理和不完備的Cholesky預處理。

（1）對角線預處理：對角線預處理指的就是M矩陣為原始系數(shù)矩陣A的對角線矩陣。在知名的開源有限元軟件FEAPpv中，就是采用的這種方法進行共軛梯度法預處理。一般情況下，這種方法求解預處理方程十分簡單，但是通常對于系數(shù)矩陣嚴格對角占優(yōu)（即對角線元素大于該行其他元素之和）的情況下才比較有效。例如之前文章中待求解的系數(shù)矩陣為：

該矩陣并不嚴格對角占優(yōu)，采用對角線預處理共軛梯度法求解結果如下：

采用對角線預處理，迭代次數(shù)并無太多改善。

但是如果系數(shù)矩陣A的對角線擴大10倍：

采用共軛梯度法求解：

而采用對角線預處理共軛梯度法求解：

      采用對角線預處理共軛梯度法僅需4次迭代即獲得收斂結果。

（2）不完備的Cholesky分解預處理

    不完備的Cholesky分解預處理指的是對系數(shù)矩陣A進行不完備的Cholesky分解。常規(guī)的Cholesky分解如下：

    其中,A是正定矩陣，L是下三角矩陣。而在不完備的Cholesky預處理共軛梯度法中，也是對A進行此分解，但是僅對于A中非0元素進行分解，A中的0元素在L中依然是0，或者甚至是求解過程中只在內存中存取A的非0元素。

   求得系數(shù)矩陣A的不完備的下三角矩陣L后，

   令

    再采用預處理共軛梯度法的具體算法獲得最終解。當然，由于L是下三角矩陣，因此預處理方程一般通過兩次“回代”（參考本公眾號文章[數(shù)值算法與編程]高斯消去法的回代部分）即可求解。

    以之前的文章共軛梯度法中的原始方程求解為例，采用不完備Cholesky分解預處理求解如下：

    

    僅需3次迭代，即獲得收斂解，而原來的常規(guī)共軛梯度法需要9次迭代。并且，matlab計算結果如下：

    通過對比不難看出，雖然僅僅是3次迭代，但是已經具備較高的求解精度。

    在實際開發(fā)中，共軛梯度法還有較多的發(fā)揮空間，比如，比如，知名有限元大師Thomas J.R. Hughes在1983年創(chuàng)立了一種基于共軛梯度法的element by element算法，這種方法不需要組裝整體剛度矩陣，而是通過逐個單元進行求解，求解效率很高，且由于不需要組裝整體剛度矩陣，計算過程中的內存需求顯著減少，并且免去了常規(guī)的采用稀疏矩陣存儲有限元剛度矩陣的組裝過程，實際上，相較于常規(guī)的矩陣，對于CSR，CSC等格式的有限元整體剛度稀疏矩陣組裝，并不是一件十分容易的事情。以上，就是共軛梯度法（二）之預處理共軛梯度法的全部內容，感謝您的閱讀！

【完】

==================================================================================

參考文獻：

     【1】《矩陣計算》，Gene H Golub，Charles F. Van Loan

==================================================================================

     

   歡迎關注公眾號 有限元術