function [alphamat,betamat,gammamat] = proposed_admm(cost_Xstr_Xttr,G1,G2,lambda_1,lambda_2)
	m = size(cost_Xstr_Xttr,1);
	n = size(cost_Xstr_Xttr,2);
	alphamat = ones(m,n)/(m*n);
	betamat = ones(n,m)/(m*n);
	gammamat = ones(m,n)/(m*n);
	D_1 = zeros(m,n);
	D_2 = zeros(m,n);
	% ADMM parameters
	iter_max_admm = 50;
	rho = 1e5;
	% rho_multi = 2;
	% rho_div = 2;
	% mu = 2;
	% rho_max = 1e6;
    % Problem description for subproblems
    problem.m = m;
    problem.n = n;
    problem.rho = rho;
    problem.lambda_1 = lambda_1;
    problem.lambda_2 = lambda_2;
    problem.G1 = G1;
    problem.G2 = G2;
    % Options for subproblems
    options.maxiter = 100;
    options.tolgradnorm = 1e-10;
    options.maxinner = 30;
    options.tolrelgradnorm = 1e-10;
    options.method = 'cg';
    options.verbosity = 0;
	for admm_iter = 1:iter_max_admm
		const1 = 0.5*(D_1 + D_2 + cost_Xstr_Xttr/problem.rho - (gammamat*G2)/n - (G1*betamat')/m);
	    problem.const1 = const1;
		alphamat = admm_multinomial_alpha(problem, alphamat(:), options);
		const2 = alphamat + D_1;
		problem.const2 = const2';
		betamat_old = betamat;
		betamat = admm_beta(problem, [], options);
		const3 = alphamat + D_2;
		problem.const3 = const3;
		gammamat_old = gammamat;
		gammamat = admm_gamma(problem, [], options);
		temp1 = (alphamat - (G1*betamat')/m);
		temp2 = (alphamat - (gammamat*G2)/n);
		temp3 = ((betamat-betamat_old)');
		temp4 = ((gammamat-gammamat_old));
		D_1 = D_1 + temp1;
		D_2 = D_2 + temp2;
		primal_residual = [temp1;temp2];
		norm_primal_residual = norm([temp1;temp2],'fro');
		norm_dual_residual = norm([temp3;temp4],'fro');
        b1 = (alphamat*ones(n,1) - (1/m)*ones(m,1));
        b2 = (alphamat'*ones(m,1) - (1/n)*ones(n,1));
		true_objective = trace(cost_Xstr_Xttr'*alphamat) + (lambda_1/2)*quad_form(b1,G1) + (lambda_2/2)*quad_form(b2,G2);
		%fprintf('[ADMM] %d, true obj: %f, primal_residual: %e, dual_residual: %e, rho: %e\n',admm_iter,true_objective,norm_primal_residual,norm_dual_residual,problem.rho);
		% % This scheme discussed in Boyd's ADMM book
		% if norm_primal_residual > mu*norm_dual_residual
		% 	rho = rho*rho_multi;
		% elseif norm_dual_residual > mu*norm_primal_residual
		% 	rho = rho/rho_div;
		% else
		% 	rho = rho;
		% end
		% This scheme discussed in Liegeois et al (2015)
		% if rho > rho_max
		% 	rho = rho_max;
		% end
		% problem.rho = rho;
	end
end