Commit eb72d8aa authored by conmccoid's avatar conmccoid
Browse files

Extrap: additional code for linear FPI example

parent 8805a59c
......@@ -16,7 +16,4 @@ for j=2:k
Hinv = [Hinv, -Hinv*h/hj; zeros(1,j-1), 1/hj];
u=Hinv(:,end); u=u/sum(u); x_out(:,j)=X(:,1:j)*u;
r_out(j)=norm(F(:,1:j)*u);
end
% if nargin==1
% x_out=x_out(:,end);
% end
\ No newline at end of file
end
\ No newline at end of file
function [x_out,r_out]=ALGO_extrap_MPE_v1_20211029(F,X)
% MPE performs minimal polynomial extrapolation using Arnoldi iteration
[d,k] = size(F);
q0 = F(:,1);
H(1,1)=norm(q0);
q0 = q0/H(1,1);
Q = zeros(d,k); Q(:,1) = q0;
x_out=zeros(d,k); x_out(:,1)= X(:,1);
r_out=zeros(1,k); r_out(1)=norm(F(:,1));
for j=2:k
Q(:,j) = F(:,j);
for i=1:j-1
H(i,j) = Q(:,j)'*Q(:,i);
Q(:,j) = Q(:,j) - H(i,j)*Q(:,i);
end
H(j,j)=norm(Q(:,j)); Q(:,j)=Q(:,j)/H(j,j);
u=zeros(j,1); u(end)=1; u=H\u; u=u/sum(u);
x_out(:,j)=X(:,1:j)*u;
r_out(j)=norm(F(:,1:j)*u);
end
\ No newline at end of file
function [x_out,r_out]=ALGO_extrap_MPE_v2_20211029(F,X)
% MPE performs minimal polynomial extrapolation using Arnoldi iteration
[d,k] = size(F);
x_out=zeros(d,k); x_out(:,1)=X(:,1);
r_out=zeros(1,k); r_out(1)=norm(F(:,1));
for j=2:k
u=[ones(1,j);F(:,1:(j-1))'*F(:,1:j)] \ [1; zeros(j-1,1)];
x_out(:,j)=X(:,1:j)*u; r_out(j)=norm(F(:,1:j)*u);
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=20;
d = 10; n=2*d;
A = rand(d); b = rand(d,1); x_true = -(A - eye(d)) \ b;
X = b;
for i=1:n
......@@ -9,28 +9,34 @@ for i=1:n
end
F = X(:,2:end) - X(:,1:n);
[x_out,r_out] = ALGO_extrap_MPE_v1_20211022(F,X(:,1:n));
x = x_out(:,end);
[x_v1,r_v1] = ALGO_extrap_MPE_v1_20211022(F,X(:,1:n));
[x_v2,r_v2] = ALGO_extrap_MPE_v2_20211029(F,X(:,1:n));
x = x_v1(:,end);
Error_abs= norm(x - x_true)
Error_rel= Error_abs/norm(x_true)
Residual = norm(A*x + b - x)
Error_list = zeros(1,n); res_MPE = zeros(1,n);
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_list(i) = norm(x_out(:,i) - x_true);
res_MPE(i) = norm(A*x_out(:,i) + b - x_out(:,i));
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
[x_GMRES,~,~,~,res_GMRES]=gmres(A-eye(d),-b);
[x_GMRES,~,~,~,res_GMRES]=gmres(A-eye(d),-b,d,0,n);
figure(1)
semilogy(1:n,Error_list,'b*--')
semilogy(1:n,Error_v1,'b*--',1:n,Error_v2,'k.--',1:length(res_GMRES),res_GMRES,'ro')
xlabel('Iteration')
ylabel('Error in solution')
title('Error vs iteration')
legend('MPE w/ Arnoldi','MPE','GMRES')
figure(2)
semilogy(1:n,res_MPE,'b*--',1:length(res_GMRES),res_GMRES,'ro')
semilogy(1:n,res_MPE_v1,'b*--',1:n,res_MPE_v2,'k.--',1:length(res_GMRES),res_GMRES,'ro')
xlabel('Iteration')
ylabel('Residual')
title('Residual vs iteration')
\ No newline at end of file
title('Residual vs iteration')
legend('MPE w/ Arnoldi','MPE','GMRES')
\ No newline at end of file
Supports Markdown
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