Commit d05638c7 authored by conmccoid's avatar conmccoid
Browse files

Extrap: added Givens rotation to MPE algo

parent 9bc6e052
......@@ -2,30 +2,40 @@ function [x_out,r_out]=ALGO_extrap_MPE_v3_20211119(F,X)
% MPE performs minimal polynomial extrapolation using Arnoldi iteration
% with Householder reflections
[d,k] = size(F);
R=F; H=zeros(d,k); %I=eye(d); Q=I;
[d,k] = size(F); k=min(d+1,k);
R=F; H=zeros(d,k); % Q=I;
x_out=zeros(d,k); r_out = zeros(k,1);
x_out(:,1)=X(:,1);r_out(1)=norm(F(:,1));
for i=1:k
if i<=d
w = SUB_extrap_Householder_v1_20211119(R(:,i),i);
% Q = Q - 2*Q*(w*w');
R = R - 2*(w*w')*R;
R = R - 2*(w*w')*R; % Q = Q - 2*Q*(w*w');
end
if i>1
ind=1:min(i-1,d);
H(ind,i-1)=R(ind,i)-R(ind,i-1);
if i<=d
H(i,i-1)=R(i,i);
end
u = zeros(min(i-1,d),1); u(1)=-R(1,1); u=H(ind,1:i-1)\u;
u = [1;u] - [u;0];
x_out(:,i) = X(:,1:i)*u;
r_out(i) = norm(F(:,1:i)*u);
if i==1
u=1;
elseif i==2
v = zeros(d+1,1); v(1)=-R(1,1); G=eye(d);
H(1,1)=R(1,2)-R(1,1); u=v(1)/H(1,1); u=[1;u]-[u;0];
elseif i>2
ind=1:i-1; indg=i-2:i-1;
H(ind,i-1)=G(ind,ind)*(R(ind,i)-R(ind,i-1)); H(i-1,i-2)=R(i-1,i-1);
g = SUB_extrap_Givens_v1_20211124(H,i-2);
H(indg,indg)=g*H(indg,indg); v(indg)=v(i-2)*g(:,1);
u=H(ind,ind)\v(ind); u = [1;u] - [u;0];
G(indg,:)=g*G(indg,:);
end
x_out(:,i) = X(:,1:i)*u;
r_out(i) = norm(F(:,1:i)*u);
end
end
% note: this is equivalent to previous formulation of HH reflections, with
% no significant change in results; also, unsure how to apply GMRES for k
% greater than d
\ No newline at end of file
% greater than d
function g = SUB_extrap_Givens_v1_20211124(A,i)
% GIVENS generates 2x2 Givens rotation matrix, zeroing out A(i+1,i)
g=A(i:i+1,i); g=g/norm(g);
g=[g(1), g(2);
-g(2), g(1)];
end
\ No newline at end of file
% Example - Extrapolation: fixed point iteration on linear function
% solve the equation x = Ax + b for random A and b
d = 10; n=2*d;
d = 100; n=2*d;
A = rand(d);
% A = A/norm(A);
A = A/norm(A);
b = rand(d,1); x_true = -(A - eye(d)) \ b;
X = b;
for i=1:n
......@@ -19,27 +19,21 @@ Error_abs= norm(x - x_true)
Error_rel= Error_abs/norm(x_true)
Residual = norm(A*x + b - x)
Error_v1 = zeros(1,n); res_MPE_v1 = zeros(1,n);
Error_v2 = zeros(1,n); res_MPE_v2 = zeros(1,n);
for i=1:n
Error_v1(i) = norm(x_v1(:,i) - x_true);
Error_v2(i) = norm(x_v2(:,i) - x_true);
res_MPE_v1(i) = norm(A*x_v1(:,i) + b - x_v1(:,i));
res_MPE_v2(i) = norm(A*x_v2(:,i) + b - x_v2(:,i));
end
Error_v1 = vecnorm(x_v1 - x_true); res_MPE_v1 = vecnorm(A*x_v1-x_v1+b);
Error_v2 = vecnorm(x_v2 - x_true);
[x_GMRES,~,~,~,res_GMRES]=gmres(A-eye(d),-b,d,0,n,eye(d),eye(d),X(:,1));
[x_GMRES,~,~,~,res_GMRES]=gmres(A-eye(d),-b,d,eps,d,eye(d),eye(d),X(:,1));
figure(1)
semilogy(1:n,Error_v1,'b*--',1:n,Error_v2,'k.--')
semilogy(1:n,Error_v1,'b*--',1:d+1,Error_v2,'k.--')
xlabel('Iteration')
ylabel('Error in solution')
title('Error vs iteration')
legend('MPE w/ Arnoldi','MPE w/ Householder')
legend('MPE w/ Arnoldi','MPE w/ HH,Givens')
figure(2)
semilogy(1:n,res_MPE_v1,'b*--',1:n,res_MPE_v2,'k.--',1:length(res_GMRES),res_GMRES,'ro')
semilogy(1:n,res_MPE_v1,'b*--',1:d+1,r_v2,'k.--',1:length(res_GMRES),res_GMRES,'ro')
xlabel('Iteration')
ylabel('Residual')
title('Residual vs iteration')
legend('MPE w/ Arnoldi','MPE w/ Householder','GMRES')
\ No newline at end of file
legend('MPE w/ Arnoldi','MPE w/ HH,Givens','GMRES')
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment