function [xout, fv, info] = PriRFed_SGD(problem, x, options)
        
    fv  = problem.cost(x);
    gf  = problem.grad(x);
    ngf = problem.M.norm(x, gf);
    fvs  = fv;
    ngfs = ngf;

    if ~isfield(options, 'stepsize_type')
        options.stepsize_type = 'fix';
    elseif ~isfield(options, 'lambda')
        options.lambda = 0.9;
    end
    
    fprintf('t:%6d|fv:%.6e(%.6e)|ngf:%.6e|ngf/ngf0:%.6e|stepsize:%.6e\n', 0, fvs(end),fvs(end)/fvs(1),ngfs(end),ngfs(end)/ngfs(1),options.stepsize);
    
    t = 0;
    xtmp = cell(options.maxiter, 1);
    num_gfs = 0; cc = 1; stepsize0 = options.stepsize;
    while t < options.maxiter        
        if (strcmp(options.stepsize_type, 'decay')) && (mod(t, options.decaysteps) == 0) && (t>1)                        
            cc = cc + 1;
            % options.stepsize = options.lambda * options.stepsize;
            options.stepsize = stepsize0 / (options.lambda + cc);
        end

        if options.option == 1
            St = 1:problem.nagents;         
        elseif options.option == 2
            St = randperm(problem.nagents, options.st);
        end
                
        localx  = cell(length(St), 1);
        weights = zeros(length(St), 1);

        count = 1;
        % active agents locally update
        for i = reshape(St, 1, length(St)) % must be a row array
            % single agent
            agent_prob = problem.agent{i};
            agent_opts = options.agent{i};             

            weights(count) = agent_prob.lncostterms;
            [localx{count}, num_gf] = RSGD(agent_prob, x, agent_opts);
            count = count + 1;
        end

        % the server aggratates
        x = tangentMean(x, localx, weights, problem.M);
        num_gfs(end+1) = num_gfs(end) + num_gf(end);
        t = t+1;

        xtmp{t} = x;
               
        xout = x;
        fv   = problem.cost(xout);
        gf   = problem.grad(xout);        
        ngf  = problem.M.norm(xout, gf);
        fvs  = [fvs;fv];
        ngfs = [ngfs;ngf];
        fprintf('t:%6d|fv:%.6e(%.6e)|ngf:%.6e|ngf/ngf0:%.6e|stepsize:%.6e\n', t, fvs(end),fvs(end)/fvs(1),ngfs(end),ngfs(end)/ngfs(1),options.stepsize);        
    end 

    info.fvs  = fvs;    
    info.ngfs = ngfs;
    info.xs   = xtmp;
    info.num_gfs = reshape(num_gfs, length(num_gfs),1);
end

function output = tangentMean(x, localx, weights, M)
    st = length(weights);
    tmp = M.zerovec(x);
    weights = weights ./ sum(weights);
    for i = 1:st     
        tmp = tmp + weights(i) * M.invretr(x, localx{i});
        % tmp = M.lincomb(x, 1., tmp, weights(i), M.invretr(x, localx{i}));
    end
    output = M.retr(x, tmp);
end


function [x, ngfs] = RSGD(prob,x,opts)    
    lncostterms = prob.lncostterms; % number of local samples
    lmaxiter = opts.lmaxiter;       % maximum local iterations
    batch_size = opts.batch_size;   % batch size
    stepsize = opts.stepsize;  
    ngfs = 0;
    for iter = 1:lmaxiter
        % Shuffle the data at the start of each epoch
        indices = randperm(lncostterms);
        ngf = 0;
        for i = 1:batch_size:lncostterms
            % mini-batch indices
            batchEnd = min(i + batch_size - 1, lncostterms);            
            batch_indices = indices(i:batchEnd);
            
            % computing gradient 
            grad = prob.partialgrad(x, batch_indices); ngf = ngf + length(batch_indices);

            % updating iterate
            x = prob.M.retr(x, grad, -stepsize);
        end
        ngfs(end+1) = ngfs(end) + ngf;
    end
end