% Matlab code for simulations and experiments in IJCNN2006 paper
% "Learning to segment any random vector"
% Aapo Hyvärinen and Jukka Perkiö
% 23 March 2006
% see end of file for release notes
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% STRUCTURAL PARAMETER SETTINGS %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%parameters for data creation and estimation
subspacesize=3; %how many components in a subspace, "k" in the paper
T=1000; %number of data points
n=50; %dimension of data
sigma2=0.01^2; % noise level, must be very small
%algorithmic parameter:
maxiter=200; %how many iterations in alternating optimization estimation
for randomdata=0:1
%when randomdata=1, create data with completely random segmentation
%when randomdata=0, create data with Boltzmann machine with intuitive structure
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% CREATE DATA ACCORDING TO GENERATIVE MODEL %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~randomdata
disp('creating Boltzmann data')
else
disp('creating random data')
end
%create model "mixing" matrix A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=orth(randn(n,subspacesize));
%create matrix B for boltzmann machine (BM) to be used in creating S
%(random data is generated by a fake boltzmann machine with
% parameter matrix B being all zeros )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if ~randomdata
%a simple beighbourhood structure
v=zeros(1,n);
v(1:floor(n/8))=1;
Bfull=toeplitz(v);
B=Bfull-diag(diag(Bfull));
else
%no structure at all
B=zeros(n);
end
%create segmentations s and observed data x
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
S=zeros(n,T);
X=zeros(n,T);
%create initial segmentations
for t=1:T
%discard segmentations that are just one segment
while sum(diff(S(:,t)))==0
%(actually, this condition is too stringent because it discards
%many other segmentations as well, but I realized this bug only
%after final paper submission so I leave it like this
%create segmentation by gibbs sampling in BM
s=sign(randn(n,1));
for biter=1:500*(1-randomdata)
i=ceil(rand*n);
y=B(i,:)*s;
s(i)=(rand<1./(1+exp(-2*y)))*2-1;
end
S(:,t)=s;
end
%create coefficients in each segment
y_plus=randn(subspacesize,1);
y_minus=randn(subspacesize,1);
%create observed data point
X(:,t)=((S(:,t)+1)/2).*(A*y_plus)+((1-S(:,t))/2).*(A*y_minus)+randn(n,1)*sqrt(sigma2);
end % of for t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ESTIMATE PARAMETERS %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%INPUT DATA TO THIS PART IS IN X which contains the x(t) as its columns
%you need to set also the structural parameters, see beginning of file
%Main results are output as
% Ahat, estimate of A
% and
% Shat, which contains estimates of s(t) as its columns
disp('estimating parameters')
%initialize parameters
%%%%%%%%%%%%%%%%%%%%%%
%uniform initial value for S
Shat=ones(n,T);
%initialize A using PCA
Cx=X*X'/T;
[E, D] = eig(Cx);
[dummy,order] = sort(diag(-D));
E = E(:,order);
Ahat=E(:,1:subspacesize);
%initialize variables
Serror=zeros(1,maxiter);
Y_plus=zeros(subspacesize,T);
Y_minus=zeros(subspacesize,T);
%START ITERATION
%%%%%%%%%%%%%%%%
for iter=1:maxiter %how many times through the data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%show progress of estimation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if rem(iter,10)==1
disp(maxiter-iter+1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimate latent variables given other parameters
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for t=1:T %for all data points
%shortcut notation
x=X(:,t);
%estimate y_plus and y_minus for a single point, given A and s
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
plusindices=find(Shat(:,t)==1);
minusindices=find(Shat(:,t)==-1);
A_plus=Ahat(plusindices,:);
A_minus=Ahat(minusindices,:);
Y_plus(:,t)=pinv(A_plus)*x(plusindices);
Y_minus(:,t)=pinv(A_minus)*x(minusindices);
if size(A_plus,1)==0, Y_plus(:,t)=zeros(subspacesize,1); end
if size(A_minus,1)==0, Y_minus(:,t)=zeros(subspacesize,1); end
%estimate s for a single point, given A and y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Shat(:,t)=sign(x.*(Ahat*(Y_plus(:,t)-Y_minus(:,t)))-.5*(Ahat*Y_plus(:,t)).^2+.5*(Ahat*Y_minus(:,t)).^2);
end %for t
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%estimate A given s and y
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Basic pseudoinverse solution:
for i=1:n;
plusindices=find(Shat(i,:)==1);
minusindices=find(Shat(i,:)==-1);
Ychosen=zeros(subspacesize,T);
Ychosen(:,plusindices)=Y_plus(:,plusindices);
Ychosen(:,minusindices)=Y_minus(:,minusindices);
Ahat(i,:)=X(i,:)*pinv(Ychosen);
end % of for i
%Ahat constrained to be orthogonal
Ahat=Ahat*sqrtm(inv(Ahat'*Ahat));
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end %of estimation iteration loop "for iter"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ANALYZE AND PLOT RESULTS %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%plot matches of estimated and real segmentations
clf
hist(abs(sum(S.*Shat/n)),20)
axis('on')
axis([0,1,0,1000])
pause
%show data for the last point (because we only know its ubest)
if ~randomdata
testpoint=1;
bar(S(:,testpoint))
axis([1,50,-2,2])
pause
bar(Shat(:,testpoint))
axis([1,50,-2,2])
pause
bar(X(:,testpoint))
axis([1,50,-2,2])
end
end %for the loop "for randomdata"
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% NOTES
%
% 1. This code does not exactly produce the simulations in the paper
% because I forgot to record all the random seeds needed.
% However, there is no qualitative difference in the results.
%
% 2. The main weakness of the method that can be seen in these
% simulations is that the noise level sigma has to be very low,
% otherwise the results are quite bad.
%
% 3. The experiments on real data in the paper were done
% with essentially the same code.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%