function [t] = solve_mmpca(Sigma, m, X_c, varargin)
% Solves the multimatrix PCA problem using the convex-concave procedure.
% That is, find X(n,m) to maximize min_i(X_c'*Sigma(:,:,i)X_c) subject to
% X_c'*X_c = I.
% Sigma is a (n x n) x p matrix of covariance matrices
% m is the number of principle components
% X_c is the initial guess
% varargin is an optional length 4 array that has values
% tau_0 initial value of tau
% tau_inc amount to increase tau at each iteration
% tau_max, maximum value of tau
% tol stopping toleratnce.
numvarargs = length(varargin);
optargs = {0.5 1.05 1000 1e-3};
optargs(1:numvarargs) = varargin;
[tau, tau_inc, tau_max, tol] = optargs{:};

% preprcoess the data
[n,~,p] = size(Sigma);
Sigma_adj_sqrt = zeros(size(Sigma));
lambda = zeros(p,1);
for i = 1:p
    lambda(i) = max(eig(Sigma(:,:,i)))+1;
    Sigma_adj_sqrt(:,:,i) = sqrtm(lambda(i)*eye(n)-Sigma(:,:,i));            
end

% perform convex concave procedure
cvx_prev = inf;
iterations = 100;
for iter = 1:iterations
    cvx_begin quiet
    variable X(n,m)
    variables Z_1(m,m) Z_2(m,m) t
    expressions tr(p)
    minimize -t+tau*trace(Z_1+Z_2)
    subject to
    for i = 1:p
        tr(i) = 0;
        for j = 1:m
            tr(i) = tr(i) + square_pos(norm(Sigma_adj_sqrt(:,:,i)*X(:,j)));
        end
        t <= -tr(i)+lambda(i)*m;
    end
    %x^Tx <= I
    [[eye(m)+Z_1 transpose(X)];
        [X eye(n)]] == semidefinite(m+n)
    %x^Tx >= I
    -X_c'*X_c+X_c'*X+transpose(X)*X_c - eye(m) + Z_2 == semidefinite(m)
    Z_1 == semidefinite(m);
    Z_2 == semidefinite(m);
    cvx_end
    
    X_c = X;
    tau = min(tau*tau_inc,tau_max);
    if((norm(cvx_optval-cvx_prev) < 1e-3 && trace(Z_1+Z_2) < tol) || strcmp(cvx_status,'Failed'))
        break
    end
    cvx_prev = -t+tau*trace(Z_1+Z_2);
end
end

