function [ Q, fnew, iter, time] = NPG_BB(B, Q_0, s, opts)
% This function solves the following optimization problem by nonmonotone proximal gradient with majorization.
% min 1/2 d^2(BQ,D),
% s.t. Q belong to C, C is a set of all r by r orthogonal matrix
% D = {B : max_{j=1,..,n} ||B_j||_0 <= s}, where B_j is the j-th row of B, B is an n by n matrix

% Input 
% B - n by n matrix 
% Q_0 - r by r orthogonal matrix
% s - sparsity level
% opts - contains the following structure variables, default values are in brackets:
%         
%           freq    - frequency for displaying output [5000]
%           maxiter - maximum number of iterations [inf]
%           disp    - display [1]: 1 - display per freq iterations
%                                  0 - no display
%           M       - memory in nonmonotone linesearch [5]
%

% Output
% 
% Q - approximate solution at termination
% iter - number of iterations
% fnew - function value at termination
% time - CPU time

tstart = tic;

% Parameters
if isfield(opts,'maxiter'), maxiter = opts.maxiter; else maxiter = 5000; end %5000, otherwise declare failure
if isfield(opts,'freq'), freq = opts.freq; else freq =5000; end
if isfield(opts,'disp'), disp = opts.disp; else disp = 1; end
if isfield(opts,'M'), M = opts.M; else M = 5; end
if isfield(opts,'tol'), tol = opts.tol; else tol = 1e-20; end

% Set matrix handles
Ahandle = @(x) B*x;
AThandle = @(x) B'*x;

% Initialization
eta = 2;
sigma = 1e-4; % line search parameter
step_init = 1; % line search parameter
Q = Q_0;
iter = 1;

% Compute function value
tmp = Ahandle(Q);
tmp_1 = tmp - sparsematrix(tmp,s) ;
fnew = 0.5*norm(tmp_1,'fro')^2;

% compute gradient
g = AThandle(tmp_1);

% Set up for nonmonotone line search
f1 = -inf(M,1);
f1(1) = fnew;
fmax = max(f1);

%% Main interation
while 1 == 1
    
    stepsize = step_init;
    z_2 = Q - stepsize*g;
    [U,~,V] = svd(z_2,'econ');
    Qnew = U*V';
    
    tmp = Ahandle(Qnew);
    tmp_1 = tmp - sparsematrix(tmp,s);
    fnew = 0.5*norm(tmp_1,'fro')^2;
    
    count = 0;
    while fnew > fmax - sigma/2*norm(Qnew - Q,'fro')^2 && fnew > min(tol,1e-20) && stepsize > 0
        % Line search
        stepsize = stepsize/eta;
        z_2 = Q - stepsize*g;
        [U,~,V] = svd(z_2,'econ');
        Qnew = U*V';
        tmp = Ahandle(Qnew);
        tmp_1 = tmp - sparsematrix(tmp,s);
        fnew = 0.5*norm(tmp_1,'fro')^2;
        count = count+1;
    end
    
    
    %% Display
    if mod(iter,freq) == 0 && disp == 1
        fprintf('iter = %5d, fval = %5.4e, ls_used = %3d, stepsize = %g\n', iter, fnew, count, stepsize);
    end
    
    
    % Brazilai-Borwein stepsize
    ds = Qnew - Q;
    gnew = AThandle(tmp_1);
    dy = gnew - g;
    dsdy = ds(:)'*dy(:);
    if  dsdy >= min(tol,1e-12)
        step_init = max(min(norm(ds,'fro')^2/dsdy,1e8),1e-8);

    else
        step_init = max(min(2.5*stepsize, 1e8), 1e-8); 
    end
    
    
    %% Update gradient
    Q = Qnew;
    g = gnew;
    
    %% Check for termination
    if iter > maxiter
        fprintf('NPG fails\n')
        break
    end
    if sqrt(fnew*2) < tol
        fprintf('Successful termination\n')
        break
    end
    if stepsize < 1e-10
        fprintf('Terminate due to lack of progress, NPG fails\n')
        break
    end
    
    
    %% Update iterate
    f1(mod(iter,M)+1,1) = fnew;
    fmax = max(f1);
    iter = iter + 1;
end
time = toc(tstart);
fprintf('iter = %d, fnew = %g, time = %g\n', iter, fnew, time)
