%% Multi-matrix Principle Component Analysis Example
% This example uses the convex-concave procedure to solve a multi-matrix
% principle component analysis problem.  That is find the directions that
% maximize the minimum variance across a set of matrices.
%
% Details can be found in section 5.4 of 
% <http://www.stanford.edu/~boyd/papers/cvx_ccv Variations and Extension of
% the Convex-Concave Procedure>.
%
% This example uses CVX, which is available <http://www.cvxr.com/cvx/ here>.
% You can download the code <mmpca_example.m> along with the
% subroutine <solve_mmpca.m>.


%%
clear all
close all

%% Problem generation parameters
p = 3; % number of matrices
n = 10; % size of matrices
m = 4; % number of principle components
sigma_var = 1; % how much to vary the diagonal entries up to +/-sigma_var/2
delta_var = 1; % how much to vary the off diagonal entries up to+/- delta_var/2
rand('state',0);

%% Algorithm parameters
tau_0 = .5;
tau_inc = 1.05;
tau_max = 1e4;
trials = 5;

%% Problem generation

% construct a base case
temp1 = orth(rand(n));
Base = temp1*diag(rand(n,1))*temp1';
sigma_var = 1;
delta_var = 1;
lambda = zeros(p,1);

% perturb the matrix
for i = 1:p
    while(true)
        pert = delta_var*2*rand(n)+(1-delta_var);
        pert = (pert+pert')/2;
        pert(1:n:end) = sigma_var*rand(n,1)+(1-sigma_var/2);
        Sigma(:,:,i) = Base.*pert;
        % make sure the matrix stays positive definite
        if(min(eig(Sigma(:,:,i))) >= .001)
            break;
        end
    end
end

%% Calculate baseline

% calculate the value from taking the best result from peforming PCA on
% each matrix individually
Mins = zeros(p,1);
for i = 1:p
     [U S V] = svd(Sigma(:,:,i));
     X_test = U(:,1:m);
     Mins(i) = trace(X_test'*Sigma(:,:,i)*X_test);
end

% calculate the value by performing PCA on the average
Sigma_sum = zeros(n,n);
avg_obj = zeros(p,1);
for i = 1:p
    Sigma_sum = Sigma_sum+Sigma(:,:,i);
end
Sigma_sum = Sigma_sum./p;
[U S V] = svd(Sigma_sum);
X_guess = U(:,1:m);
for j = 1:p
    avg_obj(j) = trace(X_guess'*Sigma(:,:,j)*X_guess);
end

% give the best value found (which is a lower bound, and an upper bound on
% the best possible value.
fprintf('The best known value is %.4f.\nAn upper bound on the optimal value is %.4f.\n',max(min(avg_obj)),min(Mins));


%% Convex-concave procedure

% use average PCA to initialize
result_avg = solve_mmpca(Sigma, m, X_guess, tau_0, tau_inc, tau_max);

fprintf('Using the average PCA initialization found a value of %.4f.\n',result_avg);

% use random initializaiton
result_best = 0;
for trial = 1:trials
    X_c = rand(n,m);
    result = solve_mmpca(Sigma, m, X_c, tau_0, tau_inc, tau_max);
    if(result > result_best)
        result_best = result;
    end
end

fprintf('After %d random initializations, a value of %.4f was found.\n',trial, result_best);