利用投影法求解向量函数的最小值问题

对于一个线性系统Ax=b而言，残差向量b-Ax的处理。

参考书籍：稀疏线性系统的迭代方法（科学出版社）P138~142

               Iterative Methods for Sparse linear systems（second edition）Yousef Saad

1. 最速下降法（Steepest Descent Method）

前提条件：A是SPD（对称正定矩阵）

其迭代过程如下：

1 Computer r= b-Ax and p=Ar
2 Until convergence, Do
3　　 a←(r,r)/(p,r)
4 　　x←x+ar
5 　　r←r-ap
6 　　Compute p:=Ar
7 EndDo

2. 极小残值迭代（Minimal Residual Iteration）

前提条件：A可以不是对称的，只需正定

其迭代过程如下：

1 Computer r= b-Ax and p=Ar
2 Until convergence, Do
3　　 a←(Ar,r)/(p,p)
4 　　x←x+ar
5 　　r←r-ap
6 　　Compute p:=Ar
7 EndDo

3. 残量范数最速下降法（Residual Norm Steepest Descent）

前提条件：A只需保证是一个非奇异矩阵。(In fact,the only requirement is that A be a nonsingular matrix)

其迭代过程如下：

1 Computer r:= b-Ax
2 Until convergence, Do
3 　　v:=A^Tr
4 　　Compute Av and /a:=||v||₂²/||Av||₂²
5 　　 x←x+av　　
6 　　r←r-aAv
7 EndDo

利用matlab实现如下：

1. 最速下降法（SD）

%实现利用SD迭代方法求解||b-A*x||的最小值
%参考书籍：稀疏线性系统的迭代方法（科学出版社）P140
%         Iterative Methods for Sparse linear systems（second edition）
%             Yousef Saad

function [x]= SD(A,b,x0)

r=b-A*x0;
TOL=0.000001;
x=x0;Num=1;
p=A*r;
while(norm(r)>TOL)
    alpha=r'*r/(p'*r);
    x=x+alpha*r;
    r=r-alpha*p;
    p=A*r;
    Num=Num+1;
end
disp('Num=');disp(Num)

测试算例：

>> clear
>> A=[2 1 1;1 2 1;1 1 2]
>> b=[14 18 20]';
>> x=[1 1 1]';
>> [x]= SD(A,b,x)

其结果显示如下：

Num=
    14


x =

    1.0000
    5.0000
    7.0000

 

2. 极小残值迭代（MR）

%实现利用MR迭代方法求解||b-A*x||的最小值
%参考书籍：稀疏线性系统的迭代方法（科学出版社）P140
%         Iterative Methods for Sparse linear systems（second edition）
%             Yousef Saad

function [x]= MR(A,b,x0)

r=b-A*x0;
TOL=0.000001;
x=x0;Num=1;
p=A*r;
while(norm(r)>TOL)
    v=A'*r;
    alpha=p'*r/norm(p)^2;
    x=x+alpha*r;
    r=r-alpha*p;
    p=A*r;
    Num=Num+1;
end
disp('Num=');disp(Num)

测试结果（算例与SD方法相同）：

Num=
    10


x =

    1.0000
    5.0000
    7.0000

 

3. 残量范数最速下降法（RNSD）

%实现利用RNSD方法求解||b-A*x||的最小值
%参考书籍：稀疏线性系统的迭代方法（科学出版社）P142
%         Iterative Methods for Sparse linear systems（second edition）
%             Yousef Saad

function [x]= RNSD(A,b,x0)

r=b-A*x0;
disp('r=')
disp(r)
TOL=0.000001;
x=x0;Num=1;
while(norm(r)>TOL)
    v=A'*r;
    Av=A*v;
    alpha=norm(v)^2/norm(Av)^2;
    x=x+alpha*v;
    r=r-alpha*Av;
    %disp('x=');disp(x)
    %disp(r)
    Num=Num+1;
    disp('Num=');disp(Num)
end

测试结果（算例与SD方法相同）：

Num=
    11


x =

    1.0000
    5.0000
    7.0000

对于一般的矩阵，如下述算例：

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

其结果为：

Num=
    56


x =

    2.7317
    2.4878
    2.0976

用方法(1)和方法(2)，将无法求解。

对于方法(3)，可将x0放在程序内部

%实现利用RNSD方法求解||b-A*x||的最小值
%参考书籍：稀疏线性系统的迭代方法（科学出版社）P142
%         Iterative Methods for Sparse linear systems（second edition）
%             Yousef Saad

function [x]= RNSD(A,b)
%初始化待求解的向量x
for i=1:length(b)
    xx0(i)=1;
end
%上述生成的xx0为一个行向量，需进行转置
x0=xx0'
r=b-A*x0;
disp('r=')
disp(r)
TOL=0.000001;
x=x0;Num=1;
while(norm(r)>TOL)
    v=A'*r;
    Av=A*v;
    alpha=norm(v)^2/norm(Av)^2;
    x=x+alpha*v;
    r=r-alpha*Av;
    %disp('x=');disp(x)
    %disp(r)
    Num=Num+1;
    disp('Num=');disp(Num)
end

在主窗口只需输入：

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