function [gammamatOpt] = admm_gamma(problem, gammainit, options)
    % Problem description
    m = problem.m;
    n = problem.n;
    const3 = problem.const3;
    G2 = problem.G2;
    
    % Local defaults for options
    localdefaults.maxiter = 50; % Max iterations.
    localdefaults.tolgradnorm = 1e-6; % Absolute tolerance on Gradnorm.
    localdefaults.maxinner = 30; % Max inner iterations for the tCG step.
    localdefaults.tolrelgradnorm = 1e-10; % Gradnorm/initGradnorm tolerance.
    localdefaults.method = 'cg'; % Default solver is trustregions (TR).
    localdefaults.verbosity = 2; % Default is print all the details.
    localdefaults.maxinnerepoch = 50; % Default max iteration for inner solver.
    localdefaults.innertolg = 1e-10; % Default norm(grad,'inf') tolerance for inner solver.
    
    % Manifold geometry
    %problem.M = euclideanfactory(n,m);
    problem.M = positivefactory(m,n);
    
    % Initialization
    if ~exist('gammainit', 'var')
        gammainit = [];
    end
    
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);
    
    % Cost
    problem.cost = @cost;
    function [f, store] = cost(gammamat, store)
        temp1 = const3 - (gammamat*G2)/n;
        f = 0.5*norm(temp1,'fro')^2;
    end
    
    % Gradient
    problem.egrad = @egrad;
    function [grad, store] = egrad(gammamat, store)
        temp1 = const3 - (gammamat*G2)/n;
        grad = -(temp1*G2)/n;
    end
    
    %     % Hessian: may be explored later for even faster convergence
    %     problem.ehess = @ehess;
    %     function [graddot, store] = ehess(gammamat, gammamatdot,store)
    %     end
    
    % Stats that we compute every iteration
    % options.statsfun = @mystatsfun;
    
    
    % Additional stopping criteria
    options.stopfun = @mystopfun;
    function stopnow = mystopfun(problem, gammamat, info, last)
        stopnow = (last >= 3 && info(last).gradnorm/info(1).gradnorm < options.tolrelgradnorm);
    end
    
    
    % Gradient and Hessian checks
    % checkgradient(problem);
    % pause;
    %     checkhessian(problem);
    %     pause;
    
    % Solver
    if strcmpi('TR', options.method)
        % Riemannian trustregions
        [gammamatOpt,~,infos] = trustregions(problem, gammainit, options);
    elseif strcmpi('CG', options.method)
        % Riemannian conjugategradients
        options.gamma_type = 'H-S';
        options.linesearch = @linesearch;
        options.ls_contraction_factor = .2;
        options.ls_optimism = 1.1;
        options.ls_suff_decr = 1e-4;
        options.ls_max_steps = 25;
        [gammamatOpt,~,infos] = conjugategradient(problem, gammainit, options);
    end
    temp1 = const3 - (gammamatOpt*G2)/n;
    fOpt = 0.5*norm(temp1,'fro')^2;
end