lans_ppsinit.m

模式识别工具包

%	lans_ppsinit	- Probabilistic Principal Surface (PPS) initialization
%
%	[pps,Lcomp,Lavg,mse]	= lans_ppsinit(y)
%		  		  lans_ppsinit(y,options)
%		  		  lans_ppsinit(y,options,state)
%
%	_____OUTPUTS____________________________________________________________
%	pps	probabilistic principal surface (PPS)		(structure)
%		.str	Description string			(string)
%		.manifold	type				(string)
%		.f	nodes f(x) in data space		( D x gx x gy)
%				.f{1}	contains grid minus poles
%				.f{2}	optionally contains [south north center]
%								(cell)
%		.x	nodes x in latent dimension		( Q x gx x gy )
%								(cell)
%			GRID information contained therein
%			GRID scaled to lie within [-1 1];
%			controlled by constants GRIDMIN, GRIDMAX
%			for 3sphere
%				.x{1}	contains grid minus poles
%				.x{2}	optionally contains [south north center]
%					
%			UNIVERSALLY use lans_md2lin to collapse!!
%		.xdim	latent grid structure of nodes		(cell)
%		.W	Weight matrix				( D x L )
%		.sigma	full covariance at each .f		( D x D x M)
%			used only if alpha~=1
%			can be derived directly from dfdx
%				sigma	= B*eye(D)+(S-B)dfdx
%			where
%				B	= (D-alpha*Q)*sig/(D-Q)
%				S	= alpha*sig;
%
%		.alpha	global attenuation factor		(scalar)
%			0 < alpha < 1	PPS (oriented GTM)
%			    alpha = 1	GTM/SOM
%			1 < alpha < D/Q manifold aligned GTM
%		.beta	global inverse spherical variance	(scalar)
%		.gamma	regularization parameter		(scalar)
%		.prior	prior probability for each ref vector	(vector)
%		.p_cen	latent basis centers			( Q x L-1 )
%			Each basis is Q-D; there are L-1 bases
%			(last one is bias)
%		.p_sig	latent basis width (1 sigma)		( 1 x L-1 )
%			max(average nearest basis distance,unexplained variance)
%			= sqrt(p_var);
%		.p_fac	factor for p_sig			(scalar)
%			1	default
%			(factor of p_var will be squared!)
%		.p	latent basis activations w.r.t. x	( L x M )
%			(l,m) represents activation of basis l
%			function of W and p
%			function of W and p
%			for latent node m
%			- Lth term denotes constant bias activation (=1)
%			- .p is constant for fixed x, p_cen, p_sig
%		.po	activation of latent Origin		( L x 1)
%			useful for "clustering"
%		.dpdx	latent basis gradient w.r.t. x		( L x Q x M )
%			row (:,:,m) represents latent basis
%			gradient at node m where the LxQ gradient matrix
%			dp/dx	= [
%				dp_1(x)/dx_1 ... dp_1(x)/dx_Q
%
%				dp_L(x)/dx_1 ... dp_L(x)/dx_Q
%				= [dP/dx_1 ... dP/dx_Q]
%			where each row
%			dp_l(X_m)/dX = (p_cen(:,l)-X_m)'*p(l,m)/(p_sigf^2)
%			- dpdx is only approximately orthogonal!!!!
%			- resulting in an approximately orthogonal dfdx!!!
%			- computed once
%			- the Lth row is 0, corresponding to constant bias
%		.dfdx	manifold gradient w.r.t. x		(D x Q x M)
%			Gradient vector in original dimension is 
%				df/dx_q	= W*dpdx(:,q,:);
%			which is a MxDxQ matrix
%			- computed only if alpha!=1
%		.D	dimensionality of data space		(integer+)
%		.Q	dimensionality of latent manifold	(integer+)
%		.M	# manifold nodes			(integer+)
%		.L	# latent basis functions (include bias)	(integer+)
%		.R	Responsibility matrix			(M x N)
%			
%	Lcomp	true log likelihood				(scalar)
%	Lavg	expected log likelihood				(scalar)
%	mse	Mean Squared Error using projection -method	(scalar)
%
%	_____INPUTS_____________________________________________________________
%	y	data in original dimension			(col vectors)
%	options	(see lanspara.m)				(string)
%		-algo	Update algorithm			(string)
%			"gtm"	approximated GTM		default
%			"gem"	Generalized EM, updating alpha
%		-alpha	attenuation factor			(scalar)
%			specifies starting alpha if "gem" algorithm used
%			0 < . < 1	PPS (oriented GTM)
%			1		GTM/SOM
%			1 < . < D/Q	manifold-aligned GTM	default
%		-center	Include center node			(scalar)
%			0	No	(default)
%			1	Yes		
%			(only works for 3sphere manifold)
%			(create a center latent node)
%		-debug	{0,1}
%		-eta	learning rate for GEM algorithm		(scalar)
%		-init	Initialization method for 3sphere	(string)
%			data*	based on data (assumed N=M)
%			manifold manifold in 3 Principal component dirs
%			gaussian about mean with radius = SIGFAC*s.d.
%			permute	randomly select M points
%
%		-L1	max # latent bases in each dimension	(integer+)
%			(excluding bias)
%			> 3 recommended
%			(not for 3sphere)
%		-L2	optional 2nd dimension
%			L1			default
%			(excluding bias)
%			(not for 3sphere)
%		-Lratio	basis follows nodes grid ratio?		(binary)
%			0 regular		default
%			1 (Yes)	
%			setting Lratio=1 uses fewer latent bases
%			(only for 3sphere & 3hemisphere)
%		-Lalpha	clamping factor for LATENT basis	(integer)
%			useful only for '3sphere' manifold
%			0 < Lalpha < 3/2
%			1			default
%		-Lfac	latent basis width factor		(scalar)
%			1			default
%		-M1	max # manifold nodes in each dimension	(integer+)
%		-M2	optional 2nd dimension
%			=M1			default
%		-manifold					(string)
%			'1line'
%				'1infiniteline'	(unimplemented)
%			'2square'
%				'2circle'	(unimplemented)
%			'3sphere'
%			'3hemisphere'
%		-method	projection method
%			see lans_ppsproj.m
%		-repeat	whether to repeat latent nodes at poles	(binary)
%			0	no (default)	
%			1	yes
%			(valid only for 3sphere manifold)
%			(latent basis ALWAYS non-repeating)
%		-gtm	GTM Toolbox compatibility mode		(scalar)
%			0	No		default
%			1	Yes		
%				adjust p_cen, p_sig
%		-eint	elevation interval			(degrees)
%			currently applicable only for '3sphere' manifold
%			10			default
%		-rint	rotation interval			(degrees)
%			currently applicable only for
%			'3sphere' and '3hemisphere' manifold
%			10			default
%	state	random state used for initialization		(col vector)
%		[]*
%
%		* default
%
%	_____NOTES______________________________________________________________
%	for demo, call function without parameters
%	- Curse of dimensionality may produce Lavg = 0 for high-D data
%	- dfdx is NOT orthogonal but span the Q-D space. It is artificially made
%	  orthonormal in the lans_pdf.m function
%	- if D=Q, alpha = 1 except for '3sphere' manifold
%	- (M) may contain additional dimensions e.g.
%				L+1 x 10 x 5 
%	- Lalpha~=1 reduces error
%
%	- For small values of L1, Lratio=1 results in a better fit, otherwise
%	  Lratio=0 is better (more regularly spaced latent basis)
%	- 3-sphere
%		- use lans_md2lin to linearize pps.x which puts poles last
%			[pps.x{1} pps.x{2}] 
%		- single latent basis @ poles
%		- # latent bases (for same x ratio) is
%			L = L1*(L1+1) - 2(L1-1) + 1 = L1^2 - L1 + 3
%	  	  or for regularly spaced bases
%	  		L = (L1-1)*L1 + 1
%	  	- latent basis width varies with location
%		- if no orientation desired, PPS will be initialized to small
%		  random Gaussian values about the mean of y
%	- actual latent basis variance is obtained as a product of p_sig,p_fac
%	  where p_sig is fixed. To alter variance, set p_fac
%	- after computation, pps.R should be set to [] in order to save memory
%	- latent basis always NON-repeating at poles
%	- if -center 1
%		- y is appended with mean(y')'
%		- pps.f{2}(:,3) denotes center node in data space
%		- an extra latent node at latent origin is created
%		- NO need to add a latent basis at latent origin
%
%	_____SEE ALSO___________________________________________________________
%	lans_pps	lans_basis	lans_pdf
%
%	(C) 2000.02.26 Kui-yu Chang
%	http://lans.ece.utexas.edu/~kuiyu

%	This program is free software; you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation; either version 2 of the License, or
%	(at your option) any later version.
%
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%
%	You should have received a copy of the GNU General Public License
%	along with this program; if not, write to the Free Software
%	Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA
%	or check
%			http://www.gnu.org/

%	_____TO DO______________________________________________________________
%	GTM compatibility implemented only for Q=1

function [pps,Lcomp,Lavg,mse]	= lans_ppsinit(y,options,state)
if nargin>0
%__________ REGULAR ____________________________________________________________
GRIDMIN	= -1;			% min value of unit grid
GRIDMAX	= 1;			% max value of unit grid
SIGFAC	= .5;

if nargin<3
	if nargin<2
		options	= [];
	end
else
	if ~isempty(state)
		if length(state>2)
			randn('state',state(1:2));% initialize random number generation
		else
			randn('state',state);
		end
		rand('state',state);		% initialize random number generation
	end
end

%__________	Get parameters
[D N]		= size(y);

debug		= paraget('-debug',options,0);
algo		= paraget('-algo',options,'gtm');
alpha		= paraget('-alpha',options,1);
center		= paraget('-center',options,0);
gamma		= paraget('-regularize',options,0);
gtm		= paraget('-gtm',options,0);
init		= paraget('-init',options,'data');
manifold	= paraget('-manifold',options);
Q		= str2num(manifold(1));
M1		= paraget('-M1',options,floor(N^(1/Q)));
M2		= paraget('-M2',options,M1);
L1		= paraget('-L1',options,round(M1/5));
L2		= paraget('-L2',options,L1);
p_fac		= paraget('-Lfac',options,1);
Lalpha		= paraget('-Lalpha',options,1);
Lratio		= paraget('-Lratio',options,0);

%__________	Initialize PPS structure
pps.str		= 'Insert remarks here';
pps.manifold	= manifold;
pps.f		= 0;
pps.x		= [];
pps.xdim	= [];
pps.W		= 0;

pps.sigma	= 0;
pps.alpha	= alpha;
pps.beta	= 1;
pps.gamma	= gamma;
pps.prior	= 1;

pps.p_cen	= 0;
pps.p_sig	= 1;
pps.p_fac	= p_fac;

pps.p		= 0;
pps.po		= 0;
pps.dpdx	= 0;
pps.dfdx	= 0;

pps.D		= D;
pps.Q		= Q;
pps.M		= M1;
pps.L		= L1;		% temp assignment, changes w.r.t. manifold
pps.R		= [];	

%__________	Check if D>Q
%if D<=Q
	%--------- different for '3sphere'
%	alpha	= 1;
%end

%__________	Initialize manifold
switch (manifold)
	case '1line'
		%__________	latent vectors x (arranged along a line) 
		pps.x{1}	= lans_scale(1:M1,GRIDMIN,GRIDMAX);
		pps.xdim	= lans_cellsize(pps.x,1);
		M		= M1;

		%__________	latent basis (also arranged along a line)