共軛梯度方法（CG）MATLAB程式設計及其和Gauss_Seidel方法的一個比較

共軛梯度法是一種迭代法，不同於Jacobi，Gauss-Seidel和SOR方法。理論上最多n步找到真解，實際計算中考慮到舍入誤差，一般迭代3n~5n步，每步的運算量相當於舉證乘向量的運算量。這種神奇的方法，對應稀疏矩陣特別有效。

% 我們來造一個大型稀疏正定矩陣，比較共軛梯度法和Gauss_Seidei方法迭代速度上的差異。
%為了方便，不妨構造一個三對角矩陣。
clc
clear
profile on
m =1000;%先來一個一千階的方程組，便於CG方法和Gauss_Seidel方法的比較。
e = ones (m,1);
A = spdiags ([ e,4*e,e ],[-1,0,1],m,m);
%full(A)
true_x = rand(m,1);
b = A * true_x;
x0 = rand(m,1);
eps = 1.0e-9;

%測試一千萬階的方程組，CG方法在本機需要若干秒。此時，Gauss_Seidel方法卡死。
[x_CG,~] = CG(A,b,x0,eps);
error_CG = norm(x_CG - true_x,2);

%來看看Gauss_Seidel方法的速度
[x_G,~] = Gauss_Seidel(A,b,x0,eps);
error_G = norm(x_G - true_x,2);

% MATLAB對於特殊的一千萬階大規模方程組的求解也幾乎是秒算，可見其矩陣運算功能的強大。
%不愧矩陣實驗室，畢竟是經過優化的。
x_matlab = A\b;
error_matlab = norm(x_matlab - true_x,2);
 profile viewer
 p = profile('info');
 profsave(p,'profile_results')

可以看得出，CG方法在解這類問題上甩Gauss_Seidel不僅僅是用幾條街來衡量的。這只是在一千階的情況下。

需要注意的是，Gauss_Seidel方法時間主要消耗在了用該方法矩陣A是否收斂的判斷上。實際當中，這一步是多餘的，因為我們可以通過若干步的迭代結果來判斷是否收斂，一般不收斂，很快就會接待到無窮，很少有解震盪不收斂的情況產生。

下面附上兩個子函式：

---

1、CG.m

function [x,steps] = CG(A,b,x0,eps)
r0 = b - A*x0;
p0 = r0;
if nargin == 3
    eps = 1.0e-6;
end
steps = 0;
while 1
    if abs(p0) < eps
        break;
    end
    steps = steps + 1;
    a0 = r0'*r0/(p0'*A*p0);%多次用到可以存一步。
    x1 = x0 + a0*p0;

    r1 = r0 -a0*A*p0;

    b0 = r1'*r1/(r0'*r0);
    %這裡的r'* r雖然後面可能還會用到，但是由於計算量不大，沒有必要再設個新變數將
    %其存下了，記憶體上的開銷划不來。
    p1 = r1 + b0*p0;

    %只是用到前後兩層的向量，所以節省記憶體開銷，計算完後面一層，可以往回覆蓋掉沒用
    %的變數。

    x0 = x1;
    r0 = r1;
    p0 = p1;

end
x = x0;
end

2、Gauss_Seidel.m

function [x,out] = Gauss_Seidel(A,b,x0,eps,iterNum)
%% 同Jacobi迭代，只是迭代函式f有所改變。
n = length(b);
if det(A) == 0
    error('No Unique Solution Or No Solution!');
end
for i = 1 : n
    j = 1;
    while A(i,i) == 0
        if A(j,i) ~= 0 && A(i,j) ~= 0
            tempV = A(j,:);
            A(j,:) = A(i,:);
            A(i,:) = tempV;
            tempS = b(j);
            b(j) = b(i);
            b(i) = tempS;
        else
            j = j + 1;
        end
        if j > n
            error('No Exchange Exist!');
        end
    end
end
D = diag(diag(A));
invD = diag(1./diag(A));
J = invD * (D - A);
invDb = invD * b;
H = eye(n) - tril(A)^-1*A;
if max(abs(eig(H))) >= 1
    error('Gauss_Seidel Algorithm Cannot Convergence！')
end

    function x = f(x)
        for l = 1:n
            temp_x = J(l,:)*x + invDb(l);
            x(l) = temp_x;
        end
    end
x = x0;
if nargin == 5
    for k = 1:iterNum-1
        x = f(x);
    end
    x_0 = x;
    x = f(x);
    out = norm(x-x_0,2);
else
    if nargin == 3
        eps = 1.0e-6;
    end
    out = 0;
    while 1
        x_0 = x;
        x = f(x);
        out = out + 1 ;
        if norm( x - x_0 ,2 ) < eps
            break;
        end
    end
end
end

對於實正定方程組而言，CG方法是收斂的。CG方法也有若干改進，比如用cholesky分解，這個後面會寫。