共轭梯度方法（CG）MATLAB编程及其和Gauss_Seidel方法的一个比较

共轭梯度方法（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分解，这个后面会写。