基于Gram-Schmidt正交法的广义极小残量法(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

 

 

              

              

posted @ 2012-09-19 23:26  liang_l  阅读(3628)  评论(0编辑  收藏  举报