Code covered by the BSD License  

Highlights from
Uncorrelated Multilinear Principal Component Analysis (UMPCA)

Uncorrelated Multilinear Principal Component Analysis (UMPCA)

by

 

The codes implement the Uncorrelated Multilinear Principal Component Analysis (UMPCA) algorithm.

UMPCA(TX,numP)
function [Us,TXmean,odrIdx]  = UMPCA(TX,numP)
% UMPCA: Uncorrelated Multilinear Principle Component Analysis
%
% %[Prototype]%
% function [Us,TXmean,odrIdx]  = UMPCA(TX,numP)
%
% %[Author Notes]%
% Author: Haiping LU
% Email : hplu@ieee.org   or   eehplu@gmail.com
% Release date: February 28, 2012 (Version 1.0)
% Please email me if you have any problem, question or suggestion
%
% %[Algorithm]%:
% This function implements the Uncorrelated Multilinear Principal Component
% Analysis (UMPCA) algorithm presented in the follwing paper:
%    Haiping Lu, K.N. Plataniotis, and A.N. Venetsanopoulos,
%    "Uncorrelated Multilinear Principal Component Analysis for Unsupervised Multilinear Subspace Learning",
%    IEEE Transactions on Neural Networks,
%    Vol. 20, No. 11, Page: 1820-1836, Nov. 2009.
% Please reference this paper when reporting work done using this code.
%
% %[Toolbox needed]%:
% Matlab Tensor Toolbox (included in this package)
% source: http://csmr.ca.sandia.gov/~tgkolda/TensorToolbox/
%
% %[Syntax]%: [Us,TXmean,odrIdx]  = UMPCA(TX,numP)
%
% %[Inputs]%:
%    TX: the input training data in tensorial representation, the last mode
%        is the sample mode. For Nth-order tensor data, TX is of 
%        (N+1)th-order with the (N+1)-mode to be the sample mode.
%        E.g., 30x20x10x100 for 100 samples of size 30x20x10
%        If your training data is too big, resulting in the "out of memory"
%        error, you could work around this problem by reading samples one 
%        by one from the harddisk, or you could email me for help.
%
%    numP: the dimension of the projected vector, denoted as P in the
%          paper. It is the number of elementary multilinear projections 
%          (EMPs) in tensor-to-vector projection.
%
%
% %[Outputs]%:
%    Us: the multilinear projection, consiting of numP (P in the paper) 
%        elementary multilinear projections (EMPs), each EMP is consisted
%        of N vectors, one in each mode 
%
%    TXmean: the mean of the input training samples TX
%
%    odrIdx: the ordering index of projected features in decreasing  
%            variance 
%
%
% %[Supported tensor order]%
% This function supports N=2,3,4, for other order N, please modify the
% codes accordingly or email hplu@ieee.org or eehplu@gmail.com for help
%
% %[Examples]%
%%%%%%%%%%%%%%%%%%%%%%%%%%Example on 2D face data%%%%%%%%%%%%%%%%%%%%%%%%%%
%       load FERETC70A15S8/FERETC70A15S8_80x80%each sample is a second-order tensor of size 80x80
%       N=ndims(fea2D)-1;%Order of the tensor sample
%       Is=size(fea2D);%80x80x721
%       numSpl=Is(3);%There are 721 face samples
%       numP=80;
%       load('FERETC70A15S8/3Train/1');%load partition for 3 samples per class
%       fea2D_Train = fea2D(:,:,trainIdx);
%       [Us,TXmean,odrIdx]  = UMPCA(fea2D,numP);
%       fea2D=fea2D-repmat(TXmean,[ones(1,N), numSpl]);%Centering
%       numP=length(odrIdx);
%       newfea = zeros(numSpl,numP);
%       for iP=1:numP
%           projFtr=ttv(tensor(fea2D),Us(:,iP),[1 2]);
%           newfea(:,iP)=projFtr.data;
%       end
%       newfea = newfea(:,odrIdx);%newfea is the final feature vector to be 
%       %fed into a standard classifier (e.g., nearest neighbor classifier)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% %[Notes]%:
% A. Developed using Matlab R2006a & Matlab Tensor Toolbox 2.1
% B. Revision history:
%       Version 1.0 released on February 28, 2012
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%TX: (N+1)-dimensional tensor Tensor Sample Dimension x NumSamples
N=ndims(TX)-1;%The order of samples.
IsTX=size(TX);
Is=IsTX(1:N);%The dimensions of the tensor

%Please see Corollary 1 in the TNN paper: numP<=min(min(Is),M}
if numP>min(Is), numP=min(Is);end
numSpl=IsTX(N+1);%Number of samples
if numSpl<numP, numP=numSpl-1;end%Centering makes the samples to be dependent, hence rank lowered by 1

%%%%%%%%%%%%%Zero-Mean%%%%%%%%%%
TXmean=mean(TX,N+1);%The mean
TX=TX-repmat(TXmean,[ones(1,N), numSpl]);%Centering

%%%%%%%%%%%%%%%UMPCA parameters%%%%%%%%%%%%%%%%
maxK=10; %maximum number of iterations, you can change this number
Us=cell(N,numP);
Us0=cell(N,1);
for iP=1:numP %Get each EMP one by one
    %Initialization
    for n=1:N
        if iP==1
            Un=ones(Is(n),1);
            Un=Un/norm(Un);
            Us0{n}=Un;
        end
        Us{n,iP}=Us0{n};
    end
    %End Initialization
    
    %Start iterations
    for k=1:maxK
        for n=1:N
            switch N
                case 2
                    switch n
                        case 1
                            Ypn=ttv(tensor(TX),Us(2,iP),2);
                        case 2
                            Ypn=ttv(tensor(TX),Us(1,iP),1);
                    end
                case 3
                    switch n
                        case 1
                            Ypn=ttv(tensor(TX),Us(2:3,iP),[2 3]);
                        case 2
                            Ypn=ttv(tensor(TX),Us([1,3],iP),[1 3]);
                        case 3
                            Ypn=ttv(tensor(TX),Us(1:2,iP),[1 2]);
                    end
                case 4
                    switch n
                        case 1
                            Ypn=ttv(tensor(TX),Us(2:4,iP),[2 3 4]);
                        case 2
                            Ypn=ttv(tensor(TX),Us([1,3,4],iP),[1 3 4]);
                        case 3
                            Ypn=ttv(tensor(TX),Us([1,2,4],iP),[1 2 4]);
                        case 3
                            Ypn=ttv(tensor(TX),Us(1:3,iP),[1 2 3]);
                    end
                otherwise
                    error('Order N not supported. Please modify the code here or email hplu@ieee.org for help.')
            end
            Ypn=Ypn.data;
            ST=zeros(Is(n));
            for i=1:numSpl
                YDiff=Ypn(:,i);
                ST=ST+YDiff*YDiff'; %Within-class Scatter
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            if iP>1
                GYYG=Gps'*Ypn'*Ypn*Gps; %equation (16) in the paper
                ThetaST=(eye(Is(n))-Ypn*Gps*inv(GYYG)*Gps'*Ypn'); %equation (15) in the paper
                ST=ThetaST*ST;
                [Lmdn,Un]=maxeig(ST);
            else
                [Lmdn,Un]=maxeig(ST);
            end
            Un=Un/norm(Un);
            Us{n,iP}=Un;
        end
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
    %%%%%%%%%%%%%%%%%%%%%%%Projection%%%%%%%%%%%%%%%%%%%%%%%%%%%
    gp=ttv(tensor(TX),Us(:,iP),1:N);
    gp=gp.data;
    if iP==1
        Gps=gp;
    else
        Gps=[Gps gp];
    end
    %%%%%%%%%%%%%%%%%%%%%%%End Projection%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

%%%%%%%%%%%%%%%%%%%%%%%%Sort by Variance%%%%%%%%%%%%%%%%%%%%%%%%%
vecYps=Gps';%vecYps contains the feature vectors for training data
Ymean=mean(vecYps,2);%Should be zero
TVars=diag(vecYps*vecYps');%Calculate variance
[stTVars,odrIdx]=sort(TVars,'descend');
odrIdx=odrIdx(1:numP); %Take the first numP
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Contact us