function [alphamatOpt] = admm_multinomial_alpha(problem, alphainit, options)
    % Problem description
    m = problem.m;
    n = problem.n;
    rho = problem.rho;
    lambda_1 = problem.lambda_1;
    lambda_2 = problem.lambda_2;
    const1 = problem.const1;
    G1 = problem.G1;
    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(m*n,1);
    problem.M = multinomialfactory(m*n,1);
    
    % Initialization
    if ~exist('alphainit', 'var')
        alphainit = [];
    end
    
    if ~exist('options', 'var') || isempty(options)
        options = struct();
    end
    options = mergeOptions(localdefaults, options);
    
    % Cost
    problem.cost = @cost;
    function [f, store] = cost(alphavec, store)
        alphamat = reshape(alphavec,m,n);
        temp1 = alphamat*ones(n,1)-ones(m,1)/m;
        temp2 = alphamat'*ones(m,1)-ones(n,1)/n;
        temp3 = alphamat + const1;
        f = rho * norm(temp3,'fro')^2 + lambda_1/2 * (temp1'*(G1*temp1)) + lambda_2/2 * (temp2'*(G2*temp2));
    end
    
    % Gradient
    problem.egrad = @egrad;
    function [grad, store] = egrad(alphavec, store)
        alphamat = reshape(alphavec,m,n);
        temp1 = alphamat*ones(n,1)-ones(m,1)/m;
        temp2 = alphamat'*ones(m,1)-ones(n,1)/n;
        temp3 = alphamat + const1;
        gradmat = 2*rho*(temp3) + lambda_1*((G1*temp1)*ones(1,n)) + lambda_2*(ones(m,1)*(temp2'*G2));
        grad = gradmat(:);
    end
    
    %     % Hessian: may be explored for even faster convergence
    %     problem.ehess = @ehess;
    %     function [graddot, store] = ehess(alphavec, alphavecdot,store)
    %     end
    
    % Stats that we compute every iteration
    % options.statsfun = @mystatsfun;
    
    
    % Additional stopping criteria
    options.stopfun = @mystopfun;
    function stopnow = mystopfun(problem, alphavec, 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
        [alphavecOpt,~,infos] = trustregions(problem, alphainit, 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;
        [alphavecOpt,~,infos] = conjugategradient(problem, alphainit, options);
    end
    alphamatOpt = reshape(alphavecOpt,m,n);
    temp1 = alphamatOpt*ones(n,1)-ones(m,1)/m;
    temp2 = alphamatOpt'*ones(m,1)-ones(n,1)/n;
    temp3 = alphamatOpt + const1;
    fOpt = rho/2 * norm(temp3,'fro')^2 + lambda_1/2 * (temp1'*(G1*temp1)) + lambda_2/2 * (temp2'*(G2*temp2));
end