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 Variations and Extension of the Convex-Concave Procedure.

This example uses CVX, which is available here. You can download the code mmpca_example.m along with the subroutine solve_mmpca.m.

Contents

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));

The best known value is 2.7870.
An upper bound on the optimal value is 3.1881.

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);

Using the average PCA initialization found a value of 2.9600.
After 5 random initializations, a value of 2.9428 was found.

Published with MATLAB® 7.10