% 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. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%