function [x, fvs, ngfs] = DPRSGD(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;  
    sigma = opts.sigma;             % standard deviation of nosie
    if ~isfield(opts, 'tau')
        tau = 1;
    else
        tau = opts.tau; 
    end
    fvs = prob.cost(x);       % log the cost of each epoch
    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);
            grad = clipping(x, grad, tau, prob.M);

            % adding tangent Guassian noise to the gradient
            vereps = tangentSpaceGaussianSample(x, prob.M, sigma);
            dpgrad = prob.M.lincomb(x, 1., grad, 1., vereps);

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