function [x, fvs, ngfs] = DPRSVRG(prob,x,opts)
    
    lncostterms = prob.lncostterms;      % number of local samples
    lmaxiter = opts.lmaxiter;            % maximum local iterations
    lmaxinneriter = opts.lmaxinneriter;  % local iterations
    batch_size = opts.batch_size;        % batch size    
    stepsize = opts.stepsize;  
    sigma = opts.sigma;             % standard deviation of nosie
    if ~isfield(opts, 'tau0')
        tau0 = 1;
    else
        tau0 = opts.tau0; 
    end
    if ~isfield(opts, 'tau1')
        tau1 = 1;
    else
        tau1 = opts.tau1; 
    end
    fvs = prob.cost(x);       % log the cost of each epoch   
    ngfs = 0;
    for iter = 1:lmaxiter        
        gradfull = prob.partialgrad(x, 1:lncostterms); % full gradient
        gradfull = clipping(x, gradfull, tau0, prob.M); ngf = lncostterms;
        xtil = x;
        % Shuffle the data at the start of each epoch
        indices = randperm(lncostterms);
        for iiter = 1:lmaxinneriter            
                % mini-batch indices
                batchEnd = min(iiter + batch_size - 1, lncostterms);            
                batch_indices = indices(iiter:batchEnd);
                
                % computing svrg
                grad1 = prob.partialgrad(x, batch_indices); ngf = ngf + length(batch_indices);
                grad1 = clipping(x, grad1, tau1, prob.M);

                grad2 = prob.partialgrad(xtil, batch_indices); ngf = ngf + length(batch_indices);
                grad2 = clipping(xtil, grad2, tau1, prob.M);

                gradtmp = prob.M.lincomb(x, 1., grad2, -1., gradfull);
                tpgradtmp = prob.M.transp(x, xtil, gradtmp);
                
                svrg = prob.M.lincomb(xtil, 1., grad1, -1., tpgradtmp);
    
                % adding tangent Guassian noise to the gradient
                vereps = tangentSpaceGaussianSample(xtil, prob.M, sigma);
                dpsvrg = prob.M.lincomb(xtil, 1., svrg, 1., vereps);
    
                % updating iterate
                xtil = prob.M.retr(xtil, dpsvrg, -stepsize);                  
        end
        x = xtil;
        ngfs(end+1) = ngfs(end) + ngf;
        fvs(end+1) = prob.cost(x);
    end
end