[Home]
MBLDA
A Demo Matlab code for MBLDA. [Code]
Reference
Zhen Wang, Yuan-Hai Shao*, Lan Bai, Chun-Na Li, Li-Ming Liu, Nai-Yang Deng. MBLDA: a novel multiple between-class linear discriminant analysis[J]. Submitted.
Main Function
function function [rX,W]= MBLDA(X,Y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% MBLDA: Multiple between-class LDA
%
% [rX,W]= MBLDA(X,Y)
%
% Input:
% X - Data matrix. Each row vector is a data point.
% Y - Data labels. Each element is a label of a point.
% Output:
% rX - Data matrix after reducing dimension. The upper bound of
% dimension is 0.5*m*(m-1), where m is the number of points.
% W - Projection matrix. rX=X*W.
% Examples:
% X = rand(50,10);
% Y = [ones(10,1);2*ones(10,1);3*ones(10,1);4*ones(10,1);5*ones(10,1)];
% [rX,W]= MBLDA(X,Y)
%
% Reference:
% Zhen Wang, Yuan-Hai Shao, Lan Bai, Chun-Na Li, Li-Ming Liu, and Nai-
% Yang Deng. MBLDA: a novel multiple between-class linear discriminant
% analysis. IEEE Transactions on Cybernetics, submitted, 2015
% Version 1.0 --Apr/2015
%
% Written by Zhen Wang (wangzhen1882@126.com)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Main
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=max(Y); % number of classes
[m,n]=size(X);
A=zeros(n,n);
meanA=zeros(t,n); % the mean of each class
for i=1:t
Xt=X(Y==i,:);
mt=size(Xt,1);
meanA(i,:)=sum(Xt,1)/mt;
A=A+Xt'*(eye(mt)-ones(mt,mt)/mt)*Xt;
end
A=A/m+eye(n)*1e-10;
nW=t*(t-1)/2;
tmpW=zeros(n,nW); % all projection
tmpL=zeros(1,nW); % all minimum lambda
k=1;
for i=1:t
for j=i+1:t
a=meanA(i,:)-meanA(j,:);
norma2=a*a';
a1_sign=find(a~=0,1);
if a1_sign~=1
a(1,1)=a(1,a1_sign);
a(1,a1_sign)=0;
tmp=A(:,1);
newA=[A(:,a1_sign),A(:,2:n)];
newA(:,a1_sign)=tmp;
tmp=newA(1,:);
newA(1,:)=newA(a1_sign,:);
newA(a1_sign,:)=tmp;
else
newA=A;
end
V=eye(n);
V(:,1)=a';
V(1,2:n)=-a(2:n)/a(1);
U=zeros(n-1,n);
for h=2:n
U(:,h)=-a(h)*a(2:n)';
U(h-1,h)=U(h-1,h)+norma2;
end
U(:,1)=-a(1)*a(2:n)';
U=U*newA*V;
for h=1:n
f=zeros(n,1);
f(h)=1;
tmp=(U*U'+1e-10*eye(n-1))\U(:,h);
b=f-U'*tmp;
if b'*b>1e-8
break;
end
end
if b(1)==0
tmpL(1,k)=inf;
tmpW(:,k)=zeros(n,1);
else
tmpL(1,k)=a*newA*V*b/(norma2^2*b(1));
tmpW(:,k)=V*b;
if a1_sign~=1
tmp=tmpW(1,k);
tmpW(1,k)=tmpW(a1_sign,k);
tmpW(a1_sign,k)=tmp;
end
end
k=k+1;
end
end
ind=find(tmpL~=inf);
tW=tmpW(:,ind);
meanAp=meanA*tW;
Dis=zeros(nW,length(ind));
for i=1:length(ind)
k=1;
for j=1:t
for h=j+1:t
Dis(k,i)=abs(meanAp(j,i)-meanAp(h,i));
k=k+1;
end
end
end
tL=min(Dis,[],1);
W=SelectW(tW,1./tL);
rX=X*W;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Select W
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fW=SelectW(tW,V)
tol=1e-10;
% (1) the best w is the smallest value in V;
% (2) v in W who abs(v'w)=1 is discarded
% (3) one direction is selected if the cosine^2 with other selected
% directions is the smallest.
[m,n]=size(tW);
W=zeros(m,n);
for i=1:n
W(:,i)=tW(:,i)/norm(tW(:,i));
end
Monitor=1:n;
[tmp,ind]=min(V);
Monitor(ind)=0; % 0 stands for selected, -1 stands for ignored
fW=[W(:,ind),zeros(m,n-1)];
for i=2:n
sM=0; % special monitor
ind=Monitor(Monitor>0);
l=length(ind);
if l<1
break;
end
val=zeros(l,1);
tmpW=W(:,Monitor==0); % selected w
len=size(tmpW,2);
for j=1:l
for r=1:len
tmp=abs(tmpW(:,r)'*W(:,ind(j)));
if abs(tmp-1)1
tlen=ind(index2);
newV=V(tlen);
[tmp,index2]=min(newV);
fW(:,i)=W(:,tlen(index2));
Monitor(tlen(index2))=0;
else
fW(:,i)=W(:,ind(index));
Monitor(ind(index))=0;
end
end
end
fW(:,n-length(find(Monitor==-1))+1:n)=[];
end
Any question or advice please email to shaoyuanhai21@163.com and wangzhen1882@126.com.
- Last updated: Apri 15, 2015