% The program demonstrates the use of GMRES algorithm for calculating
% numerical derivative from noisy data. 

% Generating the matrix describimg the forward problem (antiderivative)

n = 128;
col = ones(n,1);
row = zeros(1,n);
row(1) = 1;
A = (1/n)*toeplitz(col,row);

% Generating the noiseless data and the true derivative 

t = linspace(0,1,n)';
b0 = 3*t - 1;
b0 = max(zeros(n,1),b0);
b0 = min(b0,ones(n,1));
xtrue = 3*(t>1/3).*(t<2/3);

% Adding random noise. We assume that the norm of the noise vector can be
% assumed to be known

sigma = 0.001;
errvect = sigma*randn(n,1);
eta = norm(errvect);          %  Norm of the noise vector
b = b0 + errvect;

% Solving with GMRES. The iteration will stop if either the maximum number
% of iteration is reached, the discrepancy limit is attained or if the
% Arnoldi process stalls. 

save_iterates = 'yes';

maxit = 20;       %  Maximum number of iterations
tau = 1.5;        %  The fudge factor used in discrepancy principle
discr = tau*eta;  %  The discrepancy limit

x0 = zeros(n,1);   % Initial value for the iteration
r = b - A*x0;
beta = norm(r);
v = (1/beta)*r;
H = zeros(maxit+1,maxit);
V = zeros(n,maxit);
X = x0;
iterate = 'yes';  % Iterations stopped if iterate == 'no'
j = 0;
plotcolor = 'k';
lw = 1;
figure(1)
plot(t,xtrue,'g')
hold on

while (j<=maxit)&strcmp(iterate,'yes')
    j = j + 1;
    V(:,j) = v;
    w = A*v;
    for i = 1:j
        H(i,j) = w'*V(:,i);
        w = w - H(i,j)*V(:,i);
    end
    H(j+1,j) = norm(w);
    if abs(H(j+1,j)) < eps
        % The new direction spans no new directions, Arnoldi stalls
        iterate = 'no';
        disp(['Arnoldi stalls at iteration j =' num2str(j)]);
    else
        v = (1/H(j+1,j))*w;
    end
    e1 = zeros(j+1,1); e1(1,1) = 1;
    y = beta*(H(1:j+1,1:j)\e1);
    x = x0 + V(:,1:j)*y;
    % Checking if the discrepancy limit is attained
    dj = norm(b - A*x);
    if dj < discr
        % The discrepancy limit is reached; stop iteration
        iterate = 'no';
        disp(['Discrepancy limit attained at iteration j =' num2str(j)]);
        plotcolor = 'r';
        lw = 2;
    end
    if strcmp(save_iterates,'yes')
        % Save iterates
        X = [X,x];
    end
    % Plotting the result
    figure(1)
    plot(t,x,plotcolor,'LineWidth',lw)
    pause
end
title(['Number of iterations =' num2str(j)])
hold off

