Krylov子空间 FOM GMRES CG算法

krylov子空间代码

投影方法

给定初值 $x^{(0)}\in {\mathbb{R}}^n$,采用仿射空间 $x^{(0)}+\mathcal{K}$，其中 $\mathcal{K}$ 是krylov子空间。找到一个 $\tilde{x}\in x^{(0)}+\mathcal{K}$，令$\tilde{x}=x^{(0)}+\hat{x}$ 使得 $$r_0-A\hat{x}\perp \mathcal{L}\text{\hspace{1em}(1)}$$此时可以证明 $\tilde{x}$ 是 $ Ax=b$ 的近似解。

krylov子空间

A ∈ R n × n , r ∈ R n , K m ( A , r ) ≜ s p a n { r , A r , . . . , A m − 1 r } ⊆ R n A\in \mathbb{R}^{n\times n},r\in\mathbb{R}^n,\mathcal{K}_m(A,r)\triangleq span\{r,Ar,...,A^{m-1}r\}\subseteq\mathbb{R}^n A∈Rn×n,r∈Rn,Km​(A,r)≜span{r,Ar,...,Am−1r}⊆Rn

求krylov子空间的方法：Arnoldi过程

Arnoldi过程是一种特殊的求正交基的方法，求出来的基是 V m = s p a n { r , A r , . . . , A ( m − 1 ) r } V_m=span\{r,Ar,...,A^{(m-1)}r\} Vm​=span{r,Ar,...,A(m−1)r}，同时会计算出$H_{m+1,m}$的矩阵，$H_m=H(1:m,m)$是Hessenberg矩阵。$V_m$和 $H$之间的关系如下：$$A{V}_m=V_{m+1}{H}_{m+1,m}=V_m{H}_{m,m}+h_{m+1,m}{v}_{m+1}{e}_m^T$$$$V_m^{T}A{V}_m=H_m,V_m^{T}r=V_m^T\beta v_1=\beta e_1,\beta =\Vert r{\Vert }_2\text{\hspace{1em}(2)}$$

算法1:基于Gram-Schmidt正交化的Arnoldi过程

function [V,H]=arnoldi_GS(A,r,m)

%输入初始的向量r,K维数m,原始矩阵A;
%输出一组标准正交基Vm，hessenberg matrix H;
%一种特殊的求正交基的方法，求出来的基是span{r,Ar,...,A^(m-1)r};

n=size(A,1);
V=zeros(n,m+1);
H=sparse(m+1,m);
V(:,1)=r/norm(r);

for j=1:m 
    w=A*V(:,j); 
    for i=1:j
        H(i,j)=dot(w,V(:,i));
    end
    for i=1:j
        w=w-H(i,j)*V(:,i);
    end
    H(j+1,j)=norm(w);
    if H(j+1,j)==0
        break;
    else
        V(:,j+1)=w/H(j+1,j);
    end
end       
end

算法2:基于MGS的Arnoldi过程

function [V,H]=arnoldi_GMS(A,r,m)

%输入初始的r,和K维数m,原始矩阵A;
%输出一组标准正交基Vm，hessenberg matrix H;
%一种特殊的求正交基的方法，求出来的基是span{r,Ar,...,A^(m-1)r};

n=size(A,1);
V=zeros(n,m+1);
H=sparse(m+1,m);
V(:,1)=r/norm(r);

for j=1:m 
    w=A*V(:,j); 
    for i=1:j
        H(i,j)=dot(w,V(:,i));
        w=w-H(i,j)*V(:,i);
    end
    H(j+1,j)=norm(w);
    if H(j+1,j)==0
        break;
    else
        V(:,j+1)=w/H(j+1,j);
    end
end      
end

计算近似解

设 K m = { v 1 , . . . , v m } , L m = { w 1 , . . . , w m } , V = [ v 1 , . . . , v m ] , W = [ w 1 , . . . , w m ] \mathcal{K}_m=\{v_1,...,v_m\},\mathcal{L}_m=\{w_1,...,w_m\},V=[v_1,...,v_m],W=[w_1,...,w_m] Km​={v1​,...,vm​},Lm​={w1​,...,wm​},V=[v1​,...,vm​],W=[w1​,...,wm​]
$$\hat{x}=Vy，\tilde{x}=x^{(0)}+\hat{x}\text{\hspace{1em}(3)}$$由（1）知$$r_0-AVy\perp w_i,i=1,2...,m$$写成矩阵形式 $$W^{T}AVy=W^T{r}_0\text{\hspace{1em}(4)}$$解出 $y$ 带入（2）就是近似解。

FOM方法

我们令${\mathcal{L}}_m={\mathcal{K}}_m,W=V$,（4）式变成$V^{T}AVy=V^T{r}_0$, 由（2）知$y=\beta H_m^{-1}{e}_1$，故$\tilde{x}=x^{(0)}+\beta V_m{H}_m^{-1}{e}_1$

残差估计：$\Vert \tilde{r}{\Vert }_2=h_{m+1,m}|y(m)|$ 可以利用残差在循环中判断是否达到要求的解，这种方法是实用FOM

算法3:FOM

function [x,error]=FOM(A,x_0,m,b)
%FOM方法 初值和较小的m后就能算出来
%输出近似解x（n维）
%输入矩阵A，b,初值r，K空间的维数m；

r=b-A*x_0; 
beta=norm(r);
[V,H]=arnoldi_GS(A,r,m);
e1=eye(m,1);
y=H(1:m,:)\(beta*e1);
error=H(m+1,m)*abs(y(m,1));
x=x_0+V(:,1:size(y))*y;
end

算法4:实用FOM

function [x,error,m]=FOM_economic(A,x_0,e,b)

%FOM方法 给定误差 就能算出来
%输出近似解x（n维）
%输入矩阵A，b,初值r，误差e，K空间的初始维数m；

m=0;
r=b-A*x_0;
beta=norm(r);
V(:,1)=r/norm(r);

while(1)
    
    m=m+1;
    w=A*V(:,m); 
    for i=1:m
        H(i,m)=dot(w,V(:,i));
        w=w-H(i,m)*V(:,i);
    end  
    H(m+1,m)=norm(w);
    
    e1=eye(m,1);
    y=H(1:m,:)\(beta*e1);
    error=H(m+1,m)*abs(y(m));
    if error<e
        break;
    else
        V(:,m+1)=w/H(m+1,m);
    end
end
       
x=x_0+V(:,1:size(y))*y;

end

GMRES方法

我们令${\mathcal{L}}_m=A{\mathcal{K}}_m,W=AV$,（4）式变成当且仅当极小化问题的解$$mi{n}_{x\in x^{(0)}}\Vert b-Ax{\Vert }_2$$

算法5:实用GMRES

function [x,error,m]=GMRES_economic(A,x_0,e,b)

%给定误差 就能算出来
%输出近似解x（n维）
%输入矩阵A，b,初值r，误差e，K空间的初始维数m；

j=0;
r=b-A*x_0;
beta=norm(r);
V(:,1)=r/norm(r);
Eps=beta*eye(size(x_0,1),1);


 while(1)
    
    j=j+1;
    w(:,j)=A*V(:,j);
    
    for i=1:j
        H(i,j)=dot(w(:,j),V(:,i));
        w(:,j)=w(:,j)-H(i,j)*V(:,i);
    end  
    H(j+1,j)=norm(w(:,j));
    
    
    %以上是arnoldi过程
    
    %对H做givens变换 第几行就要做几次givens变换,第一次是给前两行, 最后givens变换是看最新生成的H(j,j)
    
    for i=1:j-1
        H(i:i+1,j)=[c(i),s(i);-s(i),c(i)]*H(i:i+1,j); 
    end
    
    if H(j+1,j)==0
        m=j;
        break;
    else
        V(:,j+1)=w(:,j)/H(j+1,j);
    end
    
    if abs(H(j,j))>abs(H(j+1,j))
        tau=H(j+1,j)/H(j,j);
        c(j)=1/sqrt(1+tau^2);
        s(j)=c(j)*tau;
    else
        tau=H(j,j)/H(j+1,j);
        s(j)=1/sqrt(1+tau^2);
        c(j)=s(j)*tau;
    end
    H(j,j)=c(j)*H(j,j)+s(j)*H(j+1,j);
    H(j+1,j)=0;
    
    
    Eps(j:j+1)=[c(j),s(j);-s(j),c(j)]*[Eps(j);0];
    
    error=abs(Eps(j+1))*beta;
    if abs(Eps(j+1))*beta<e
        m=j;
        break
    end
 end


y=H(1:m,1:m)\(Eps(1:m));
x=x_0+V(:,1:size(y))*y;

end

算法6:共轭梯度法

function [x,m, r2] = CG(A, b, error, max_n)

	
	n = size(A, 1);
	x = zeros(n, 1);
	r = b - A * x;
	p = r;
	
	r2     = r' * r;
	r2_eps = r2 * error * error;
	%r2_old = r2;
	m = 1;
	r_norm(m) = norm(r);
	
	while ((m < max_n) && (r2 > r2_eps))
		s = A * p;
		alpha = r2 / (p' * s);

		x = x + alpha * p;
		r = r - alpha * s;
		
		r2_old = r2;
		r2 = r' * r;
		
		beta = r2 / r2_old;
		p = r + beta * p;
		
		m = m + 1;
		r_norm(m) = norm(r, 2);
	end
	
end

算例结果

testA的规模是$2304\times 2304$ ，当误差小于$1e-9$结束迭代