共軛梯度（CG）方法

簡單介紹

共軛梯度方法也是一種迭代方法，不同于Jacobi，Gauss-Seidel和SOR方法，理論上只要n步就能找到真解，實際計算中，考慮到舍入誤差，一般迭代3n到5n步，每步的運算量相當與矩陣乘向量的運算量，對稀疏矩陣特別有效。
共軛梯度方法對于求解大型稀疏矩陣是很棒的方法，但是這個方法看起來總不是太靠譜。這個方法也不是越迭代精度越高，有時候可能迭代多了，反而出錯，對迭代終止條件的選擇，要求還是很高的。
共軛梯度法收斂的快慢依賴于系數矩陣的譜分布情況，當特征值比較集中，系數矩陣的條件數很小，共軛梯度方法收斂得就快。“超線性收斂性”告訴我們，實際當中，我們往往需要更少的步數就能得到所需的精度的解。

基本思想和原理核心

考慮線性對稱正定方程組： $$Ax=b$$ 可以轉化為求解二次泛函 $$?(x)=\frac{1}{2}{x}^{T}Ax?b^{T}x。$$
的最小值問題，直接可以驗證${?}^{\prime }(x)=Ax?b$。
當然，選擇其他的二次泛函也能滿足這個條件導數等于這個的條件了，那是另外的一些方法了，這里不提。

比如說這里取$A=diag(2,2),b=[2,2{]}^T$，那么就有 $$\frac{1}{2}{x}^{T}Ax?bTx=x_1^2+x_2^2?2x_1?2x_2=(x_1?1{)}^2+(x_2?1{)}^2?2$$
這里可以看出它是一個尖朝下，經過源點的橢圓拋物面（當然，這個例子中橢圓是圓），且在$(1,1)$點達到最小值-2。如下圖所示（手繪的，比較丑，不過我已經很努力了）：

在介紹共軛梯度法前，我們得先說一個這樣一個二次泛函的性質，用大白話來說，就是這樣：

• 用一根xy平面上的直線穿過二次曲面在xy平面上的值形成的橢圓，那么這個直線上的二次泛函的值在直線被橢圓截出的線段的中點處取到最小。且這個最小值是唯一的，如圖。
  

• 最小值的梯度方向垂直于上面提到的那根直線。

• 最小值點位于位于這根直線的共軛超平面上（一維），不知道什么是共軛超平面，沒有涉及到公式推導計算，這里可先不用管，只要知道它是由這根直線的法平面和$A$以及真解有關的一個平面即可。

這個其實很好理解，對于二次泛函形成的函數圖像，我們沿著每個方向用一個平行于z軸的平面去切它，得到的是一個拋物面，放到坐標平面上，就是一根拋物線，拋物線的最小值點當然在中點處取到。注意，這里的負梯度方向，就是沿最速下降方向，它不一定能經過橢圓的圓心。

以上是以最簡單的二維的情況來說明這個問題，更高維的情況也如是。比如說三維的情況，我們可以想象有一個西瓜，它的內部每一點處有某一種性質，比如說甜度（或密度等等），假設它在正中心處甜度最高，最外面是瓜皮最不甜，依次往里甜度是呈拋物增長的。那么，我們拿一把刀，給它切一刀（不一定對半），會形成一個切面，那么這個切面是個橢圓，且在橢圓的正中間是最甜的（極值點），拿一根牙簽，從這個最甜的點垂直切面往里面一戳，這個方向就是這一點負梯度方向。當然，牙簽夠長，你也不一定能正好西瓜正中心。因此才有了從最速下降法到共軛梯度法的改進。

有了以上的思想，我們很容易就能推導出共軛梯度法的算法過程，不想擺公式，還是以切西瓜為例：

• 首先我們拿一根足夠長的牙簽沿著西瓜任意位置將它穿透，那么前面的性質告訴我們沒在西瓜中的牙簽的中點就是這根牙簽上的最值點。（這其實就是最速下降方法的步驟，最速下降法往往是沿著一個方向取到最小，步長以此來決定，接著換一個方向，同樣一個過程……）

• 接著找到最值點（牙簽中點）的梯度方向，最速下降法就是以梯度方向接著找最小值點，但我們現在不這么做。梯度方向和原來的牙簽的方向形成了一個面，我們試圖在這個面里面找一個更好的方向。前面的性質告訴我們，這個切面的中點是這個面的最小值點，那么我就應該以牙簽中點和這個切面的連線作為方向是最理想的。

• 問題是這個切面的中點不好找，好在前面的性質告訴我們，這個切面的中點的梯度方向是垂直于這個切面的，把提到的這幾個條件聯立起來，其實很快就能找到牙簽中點和切面中點的連線方向，以及我們所需要的步長。

把以上的過程用公式擺開，就得到了共軛梯度法算法過程，如下算法框圖所示：

當然，共軛梯度算法并不是很快，后來人們提出了很多方法去加快這個速度，比如說各種預優方法等等，都是后話，這里不提。

程序代碼和結果

C代碼

#include<stdio.h> #include<stdlib.h> #include<math.h> #include<string.h> #define N 5 #define epsilon 0.005void main() {void matrixTimesVec(double A[N][N], double b[N], double Ab[N]);double scalarProduct(double vec1[], double vec2[]);void vecPlus(double vec1[], double vec2[], double vec[]);void numPlusVec(double num, double vec0[], double vec[]);int i, j;static double A[N][N] = { 0 };static double b[N] = { 1, 1, 1, 1, 1 };static double x0[N] = { 1, 1, 1, 1, 1 };double x[N], r[N], p[N], w[N], alpha, rho00, rho0, rho1, beta;//生成一個大規模稀疏矩陣A，這里以三對角為例。for (i = 1; i < N - 1; i++){A[i][i - 1] = 2;A[i][i] = 3;A[i][i + 1] = 1;}A[0][0] = 3; A[0][1] = 1;A[N - 1][N - 2] = 2; A[N - 1][N - 1] = 3;printf("\n要求解的示例方程組為：\n A ||| b ||| x0\n");for (i = 0; i < N; i++){for (j = 0; j < N; j++){printf("%f ", A[i][j]);}printf("||| %f||| %f\n", b[i], x0[i]);}//initmemcpy(x, x0, N*sizeof(double));double Ax0[N];matrixTimesVec(A, x0, Ax0);numPlusVec(-1.0, Ax0, Ax0);vecPlus(b, Ax0, r);rho0 = scalarProduct(r, r);rho00 = rho0;memcpy(p, r, N*sizeof(double));do{matrixTimesVec(A, p, w);alpha = rho0 / (scalarProduct(p, w));double temp[N];numPlusVec(alpha, p, temp);double x_temp[N];vecPlus(x, temp, x_temp);memcpy(x, x_temp, N*sizeof(double));numPlusVec(-alpha, w, temp);double r_temp[N];vecPlus(r, temp, r_temp);memcpy(r, r_temp, N*sizeof(double));rho1 = scalarProduct(r, r);beta = rho1 / rho0;numPlusVec(beta, p, temp);vecPlus(r, temp, p);rho0 = rho1;} while (rho1 > epsilon);printf("\n\n");printf("方程組的解為：\n");for (i = 0; i < N; i++)printf("%f\n", x[i]);getchar(); }void matrixTimesVec(double A[N][N], double b[N], double Ab[N]) {int i, j;for (i = 0; i < N; i++){Ab[i] = 0.0;for (j = 0; j < N; j++){Ab[i] = Ab[i] + A[i][j] * b[j];}} }double scalarProduct(double vec1[], double vec2[]) {double s = 0;int i;for (i = 0; i < N; i++){s = s + vec1[i] * vec2[i];}return s; } void vecPlus(double vec1[], double vec2[], double vec[]) {int i;for (i = 0; i < N; i++){vec[i] = vec1[i] + vec2[i];} } void numPlusVec(double num, double vec0[], double vec[]) {int i;for (i = 0; i < N; i++)vec[i] = num*vec0[i];}

MATLAB代碼

clc clear %% 共軛梯度算法并不是迭代的次數越多越精確，所以要合理地設置迭代終止的的條件。 epsilon = 0.005; N = 5; A = zeros(N); for i=2:N-1A(i,i-1)=2;A(i,i) = 3;A(i,i+1) =1; end A(1,1)=3; A(1,2) = 1; A(N,N-1) = 2; A(N,N) = 3; b = ones(N,1); x0 = b;%% 初始化 x = x0; r = b - A*x0; rho0 = conj(r')*r; rho00 = rho0; p=r;while(rho0>0.005)w= A*p;alpha = rho0/(conj(p')*w);x = x+alpha*p;r = r- alpha*w;rho1 = conj(r')*r;beta = rho1/rho0;p = r + beta*p;rho0 = rho1;xend A\b x

總結

以上是生活随笔為你收集整理的共轭梯度（CG）算法的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。