function [w, infos] = HSGD(problem, options)
 
% set dimensions and samples
d = problem.dim();
n = problem.samples();

% extract options
if ~isfield(options, 'step_init')
    step_init = 0.1;
else
    step_init = options.step_init;
end
step = step_init;

if ~isfield(options, 'step_alg')
    step_alg = 'fix';
else
    if strcmp(options.step_alg, 'decay')
        step_alg = 'decay';
    elseif strcmp(options.step_alg, 'decay-2')
        step_alg = 'decay-2';
    elseif strcmp(options.step_alg, 'fix')
        step_alg = 'fix';
    else
        step_alg = 'decay';
    end
end

if ~isfield(options, 'lambda')
    lambda = 0.1;
else
    lambda = options.lambda;
end

if ~isfield(options, 'tol_optgap')
    tol_optgap = 1.0e-12;
else
    tol_optgap = options.tol_optgap;
end

if ~isfield(options, 'batch_size')
    batch_size = 10;
else
    batch_size = options.batch_size;
end


if ~isfield(options, 'max_epoch')
    max_epoch = 100;
else
    max_epoch = options.max_epoch;
end

if ~isfield(options, 'w_init')
    w = randn(d,1);
else
    w = options.w_init;
end

if ~isfield(options, 'f_opt')
    f_opt = -Inf;
else
    f_opt = options.f_opt;
end

if ~isfield(options, 'permute_on')
    permute_on = 1;
else
    permute_on = options.permute_on;
end

if ~isfield(options, 'verbose')
    verbose = false;
else
    verbose = options.verbose;
end

if ~isfield(options, 'store_w')
    store_w = false;
else
    store_w = options.store_w;
end


% initialize
epoch = 0;
grad_calc_count = 0;

% store first infos
clear infos;
infos.iter = epoch;
infos.place=0;
infos.time = [0];
infos.time2 = [0];
infos.grad_calc_count = grad_calc_count;
f_val = problem.cost(w);
optgap = f_val - f_opt;
infos.optgap = optgap;
infos.cost = f_val;
if store_w
    infos.w = w;
end

% display infos
if verbose > 0
    fprintf('HSGD: Epoch = %03d, optgap = %.4e\n', epoch,  optgap);
end

k=options.k;
% set start time
% main loop
while (optgap > tol_optgap) && (epoch < max_epoch)
    
    % permute samples
    if permute_on
        perm_idx = randperm(n);
    else
        perm_idx = 1:n;
    end
    num_of_bachces = ...
    ceil(log((options.zeta-1)*n/(options.tau)+1)/log(options.zeta))+1; %floor(n / batch_size);
    start_index=1;
    
    for j=1:num_of_bachces
        %% calculate batch size and index 
        batch_size_s=floor(options.tau*(options.zeta^(j-1)));
            
        batch_size_s=min(batch_size_s,n-start_index+1);% to avoid the number beyond the total number

        indice_j = perm_idx(start_index:start_index+batch_size_s-1);
        start_index=start_index+batch_size_s;
        
        
        %% upate w
        grad =  problem.grad(w, indice_j);
        w = w - step * grad;
        w = hardthreshold(w, k, d);
         
         
        %% store info
        infos.time=[infos.time start_index-1+n*epoch];
        f_val = problem.cost(w);
        optgap = f_val - f_opt;
        infos.optgap = [infos.optgap optgap];
        infos.cost = [infos.cost f_val];
        
        if store_w
            infos.w = [infos.w w];
        end
        if start_index>n
            break;
        end
    end
    % update epoch
    epoch = epoch + 1;
    % display infos
    if verbose > 0
        fprintf('HSGHT: Epoch = %03d,  optgap = %.4e\n', epoch, optgap);
    end
end
end
