function [betamatOpt] = admm_beta(problem, betainit, options)
    % Problem description
    m = problem.m;
    n = problem.n;
    const2 = problem.const2;
    G1 = problem.G1;
    
    % 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(n,m);
    
    % Initialization
    if ~exist('betainit', 'var')
        betainit = [];
    end
    
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);
    
    % Cost
    problem.cost = @cost;
    function [f, store] = cost(betamat, store)
        temp1 = const2 - (betamat*G1)/m;
        f = 0.5*norm(temp1,'fro')^2;
    end
    
    % Gradient
    problem.egrad = @egrad;
    function [grad, store] = egrad(betamat, store)
        temp1 = const2 - (betamat*G1)/m;
        grad = -(temp1*G1)/m;
    end
    
    %     % Hessian: may be explored for even faster convergence
    %     problem.ehess = @ehess;
    %     function [graddot, store] = ehess(betamat, betamatdot,store)
    %     end
    
    % Stats that we compute every iteration
    % options.statsfun = @mystatsfun;
    
    
    % Additional stopping criteria
    options.stopfun = @mystopfun;
    function stopnow = mystopfun(problem, betamat, 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
        [betamatOpt,~,infos] = trustregions(problem, betainit, options);
    elseif strcmpi('CG', options.method)
        % Riemannian conjugategradients
        options.beta_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;
        [betamatOpt,~,infos] = conjugategradient(problem, betainit, options);
    end
    temp1 = const2 - (betamatOpt*G1)/m;
    fOpt = 0.5*norm(temp1,'fro')^2;
end