基于GramSchmidt正交法的广义极小残量法(GMRES)

总体的目标：求解线性方程组Ax=b

待解的问题：当矩阵的维数较大时，不能用常规的直接分解法或者是一般的迭代法

采取的方法：广义极小残量法(GMRES)

参考文献：Iterative Methods for Sparse linear systems（second edition） Yousef Saad P165

               矩阵论简明教程（第二版）科学出版社 徐仲等编 P150 6.3.2节 A+在解线性方程组中的应用

               （Moore-Penrose 逆的使用）

 相关链接：http://blog.sciencenet.cn/home.php?mod=space&uid=315774&do=blog&id=383334

                （里面有Gram-Schmidt正交法思想的由来）

采用matlab编程实现如下：

%实现基于FOM方法和Gram-Schmidt正交法做GMRES（广义极小残量法）
%参考书籍：稀疏线性系统的迭代方法（科学出版社）P165
%         Iterative Methods for Sparse linear systems（second edition）
%             Yousef Saad
function [x]=GMRESGS(A,b)
%初始化待求解的向量x
for i=1:length(b)
    xx0(i)=1;
end
%上述生成的xx0为一个行向量，需进行转置
x0=xx0'
r0=b-A*x0;
beta=norm2(r0);
v(:,1)=r0/beta;
[n,m]=size(A);%n为A的行，m为A的列 The row n and column m of the matrix
for j=1:m
    w(:,j)=A*v(:,j);
    for i=1:j
        h(i,j)=w(:,j)'*v(:,i);%h_ij=(w_j,v_i)
        w(:,j)=w(:,j)-h(i,j)*v(:,i);
    end
    h(j+1,j)=norm(w(:,j));
    if(h(j+1,j)==0)
        m=j;
        %Hmbar;
        for row=1:m
         for col=1:m
           Hmbar(row,col)=h(row,col);
         end
        end
        break;
    end
    v(:,j+1)=w(:,j)/h(j+1,j);
end
%Hmbar;
for row=1:m
    for col=1:m
        Hmbar(row,col)=h(row,col);
    end
end
Hmbar
for i=1:m
    e1(1,i)=0;
end
e1(1,1)=1;
e2=beta*e1'
y(:,m)=pinv(Hmbar)*e2;%pinv代表Moore-Penrose逆，来求解最小二乘解
x(:,m)=x0+v(:,1:m)*y(:,m);%此处v矩阵为m*(m+1),由m+1个正交向量构成，v(:,1:m)为m*m

测试算例：

>>clear
>>A=[1 2 3;4 2 1;2 5 1];
>>b=[14 18 20]';
>>[x]=GMRESGS(A,b)

最后一列即为结果

基于重启技术的GMRES方法（Matlab）

参考资料：稀疏线性系统的迭代方法（科学出版社）P172

%包含有重启技术，节省存储空间
%实现基于FOM方法和Gram-Schmidt正交法做GMRES（广义极小残量法）
%参考书籍：稀疏线性系统的迭代方法（科学出版社）P165
%         Iterative Methods for Sparse linear systems（second edition）
%             Yousef Saad
function [x]=GMRESGS(A,b)
%初始化待求解的向量x
for i=1:length(b)
    xx0(i)=1;
end
%上述生成的xx0为一个行向量，需进行转置
x0=xx0'
[n,m]=size(A);%n为A的行，m为A的列 The row n and column m of the matrix
TOL=0.00001;%误差控制
r0=b-A*x0;
beta=norm(r0);
v(:,1)=r0/beta;
Num=0;
while(beta>TOL) %残差的范数作为判断的条件
    for j=1:m
        w(:,j)=A*v(:,j);
        for i=1:j
            h(i,j)=w(:,j)'*v(:,i);%h_ij=(w_j,v_i)
            w(:,j)=w(:,j)-h(i,j)*v(:,i);Num=Num+1;
        end
        h(j+1,j)=norm(w(:,j));
        if(h(j+1,j)==0)
            m=j;
            %Hmbar;
            for row=1:m
                for col=1:m
                    Hmbar(row,col)=h(row,col);Num=Num+1;
                end
            end
            break;
        end
        v(:,j+1)=w(:,j)/h(j+1,j);
    end
    %Hmbar;
    for row=1:m
        for col=1:m
            Hmbar(row,col)=h(row,col);Num=Num+1;
        end
    end
    Hmbar
    for i=1:m
        e1(1,i)=0;Num=Num+1;
    end
    e1(1,1)=1;
    e2=beta*e1';
    y(:,m)=pinv(Hmbar)*e2;%pinv代表Moore-Penrose逆，来求解最小二乘解
    x=x0+v(:,1:m)*y(:,m);%此处v矩阵为m*(m+1),由m+1个正交向量构成，v(:,1:m)为m*m
    r0=b-A*x;%重启开始
    beta=norm(r0);
    v(:,1)=r0/beta;
    Num=Num+1;
end
disp(Num)

测试算例与上同

结果显示如下：

>> clear
>> A=[1 2 3;4 2 1;2 5 1];
>> b=[14 18 20]';
>> [x]=GMRESGS(A,b)

x0 =

     1
     1
     1


Hmbar =

    6.8389   -0.0096    2.5868
    0.8055   -1.5244   -1.7596
         0    1.9342   -1.3145

    19


x =

    2.7317
    2.4878
    2.0976