Code covered by the BSD License

# PCA and ICA Package

### Brian Moore (view profile)

24 Sep 2012 (Updated )

Implements Principal Component Analysis (PCA) and Independent Component Analysis (ICA).

myICA(z,NUM,varargin)
```function [z_ic A T mean_z] = myICA(z,NUM,varargin)
%--------------------------------------------------------------------------
% Syntax:       z_ic = myICA(z,NUM);
%               z_ic = myICA(z,NUM,'true');
%               z_ic = myICA(z,NUM,'false');
%               [z_ic A T mean_z] = myICA(z,NUM);
%               [z_ic A T mean_z] = myICA(z,NUM,'true');
%               [z_ic A T mean_z] = myICA(z,NUM,'false');
%
% Inputs:       z is an M x N matrix containing N samples of an
%               M-dimensional multivariate random variable
%
%               NUM is the desired number of independent components.
%
%               display can be {'true','false'}. The default is 'true'.
%
% Outputs:      z_ic is a NUM x N matrix containing the NUM independent
%               components (scaled to have variance 1) of each of the N
%               samples in z.
%
%               A and T are the ICA transformation matrices such that:
%               z_LD = T \ pinv(A) * z_ic + repmat(mean_z,1,size(z,2));
%               is the NUM-dimensional ICA approximation of z
%
%               mean_z is the M x 1 sample mean vector of z.
%
% Description:  This function performs independent component analysis (ICA)
%               on the input samples of a multivariate random variable and
%               returns the NUM independent components of each sample.
%
% Author:       Brian Moore
%               brimoor@umich.edu
%
% Date:         July 20, 2012
%--------------------------------------------------------------------------

eps = 1e-4; % Convergence criteria
maxSamples = 1000; % Max number of data points in sample mean calculations
maxIters = 100; % Maximum number of iterations allowed

% Parse user input
if (nargin == 3)
display = varargin{1};
else
display = 'true';
end

% Center and whiten the input data
[z_cw T mean_z] = myCenterAndWhiten(z);

% Get data dimension
[m,n] = size(z_cw);

if (n > maxSamples)
% Draw maxSamples data points at random for sample mean calculations
z_cw_trimmed = z_cw(:,randperm(n,maxSamples));
else
% Use all data points for sample mean calculations
z_cw_trimmed = z_cw;
end

% Choose random weights initially
w = rand(NUM,m);
for i = 1:NUM
w(i,:) = w(i,:) / norm(w(i,:));
end

% Initialize loop variables
err = ones(NUM,1);
its = 0;

while ((max(err) > eps) && (its < maxIters))
% Increment iteration counter
its = its + 1;

% Save last weight matrix
w_old = w;

% for each weight vector
for i = 1:NUM
% Last independent components
si = w_old(i,:) * z_cw_trimmed;

% Compute negentropy scores
% negentropy function : f(u) = -exp(-u^2/2)
g = si .* exp(-0.5 * (si.^2));
gp = -1.0 * ((si.^2) .* exp(-0.5 * (si.^2)));

% Update weights in the direction of maximum negentropy
w(i,:) = mean(z_cw_trimmed .* repmat(g,m,1),2)' - mean(gp) * w_old(i,:);

% Normalize weight vector
w(i,:) = w(i,:) / norm(w(i,:));
end

% Decorrelate weight vectors
[U,S,~] = svd(w,'econ');
Sinv = diag(1./diag(S));
w = U * Sinv * U' * w;

% Compute innovation
for i = 1:NUM
err(i) = 1 - w(i,:) * w_old(i,:)';
end

% Display innovation
if strcmpi(display,'true')
disp(['Iteration ' num2str(its) ': max(1 - <w' num2str(its) ',w' num2str(its-1) '>) = ' num2str(max(err))])
end
end

% Return transformation matrix
A = w;

% Compute independent components
z_ic = A * z_cw;
```