Commit bbfba5a0 authored by Conor McCoid's avatar Conor McCoid
Browse files

Extrap: initial commit on MPE algo and example

parent 51146a6d
function x_out=ALGO_extrap_MPE_v1_20211022(F,X)
% MPE performs minimal polynomial extrapolation using Arnoldi iteration
[d,k] = size(F);
q0 = F(:,1); Hinv = 1/norm(q0); q0 = Hinv*q0;
Q = zeros(d,k); Q(:,1) = q0; x_out=zeros(d,k); x_out(:,1)= X(:,1);
for j=2:k
Q(:,j) = F(:,j);
h = zeros(j-1,1);
for i=1:j-1
h(i) = Q(:,j)'*Q(:,i);
Q(:,j) = Q(:,j) - h(i)*Q(:,i);
end
hj = norm(Q(:,j)); Q(:,j) = Q(:,j)/hj;
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;
end
% if nargin==1
% x_out=x_out(:,end);
% 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;
A = rand(d); b = rand(d,1); x_true = -(A - eye(d)) \ b;
X = b;
for i=1:n
X(:,i+1) = A*X(:,i) + b;
end
F = X(:,2:end) - X(:,1:n);
x_out = ALGO_extrap_MPE_v1_20211022(F,X(:,1:n));
x = x_out(:,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);
for i=1:n
Error_list(i) = norm(x_out(:,i) - x_true);
end
semilogy(1:n,Error_list,'b*--')
\ 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