% 
% Written by:
% -- 
% John L. Weatherwax                2005-08-14
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; drawnow; 
clc; 
clear; 

% generate some data (by default we have a zero mean): 
n = 300; p = 5; 
n = 10000; p = 15; 
n = 50000; p = 50; 
X = randn(n,p); 

% get the covariance (and eigendecomposition): 
% here A = eigvec
%      L = eigval 
covm = cov( X ); 
[eigvec,eigval] = eig( covm ); 

% reconstruct covm from its eigendecomposition (covm = A L A^T = eigvec * eigval * (eigvec.')): 
Xorig = eigvec * eigval * ( eigvec.' ); 
norm( covm - Xorig ) 

% compute the eigenvalues from the eigenvectors and the original covariance matrix covm (L = A^T covm A) 
Lorig = ( eigvec.' ) * covm * ( eigvec );
norm( eigval - Lorig )