% created by Nils Hase 2016 (nilshase@math.uni-bremen.de)

% Tikhonov regularization for problems 
%
% [x_opt,alpha_opt] = arg min 1/2 ||Ax-y||_2^2 + alpha/2 ||x||_2^2
% 
% solved via a conjugate gradient method for the normal equation with
% accelaration by Frommer & Maaß (1999) coupled with Morozov's discrepancy
% principle to determine a suitable regularization parameter


function [x_opt,alpha_opt] = TP_CGLS_fast(A,y,alpha,delta,plotter)

tau = 1;
n = size(A,2);

discrepancy = zeros(size(alpha));
x_opt = zeros(n,1);
x = zeros(n,1);

if plotter == 1; figure(); end

for k = 1:length(alpha)
    % Calculate estimate to fixed regularization parameter using an
    % accelerated CGLS-method
    alpha_opt = alpha(k);
    
    z = y - A*x;
    r = A'*z - alpha_opt*(x);
    p = r;
    kk = 0;
    
    while (kk < sqrt(n))
        % CGLS iteration
        kk = kk + 1;
        x_old = x;        
        q = A*p;
        gamma = norm(r)^2/(norm(q)^2+alpha_opt*norm(p)^2);
        x = x + gamma*p;
        
        z = z - gamma*q;
        r_old = r;
        r = A'*z - alpha_opt*x;
        omega = norm(r)^2/norm(r_old)^2;
        p = r + omega*p;
            
        if (norm(z) - 1/(2*sqrt(alpha(k)))*norm(r) > tau*delta*sqrt(length(y)))
            % Stop because of acceleration (see Frommer & Maaß (1999))
            break;
        elseif (x_old == x)
            % Stop because convergence is reached
            break;
        elseif (kk > sqrt(n))
            % Stop because maximum number of iterations is reached
        end
    end
    x_opt = x;
    
    discrepancy(k) = norm(z)/sqrt(length(y));
    % Morozov's discrepancy principle
    if (discrepancy(k) < tau*delta)
        break;
    end        
end

if plotter == 1
    figure();
    plot(log10(alpha),discrepancy)
    hold on;
    plot([log10(alpha(1)),log10(alpha(end))],tau*delta*[1,1],'r--')
    title('Discrepancy curve CGLS-fast')
    hold off;
end

end
