lans_ppsgtm1.m

模式识别工具包

%	lans_ppsgtm1	- Generative Topographical Mapping approximation for PPS
%
%	[pps,Lcomp,Lavg]	= lans_ppsgtm1(pps,y,options)
%
%	_____OUTPUTS____________________________________________________________
%	pps	updated PPS					(structure)
%	Lcomp	updated true log likelihood			(scalar)
%	Lavg	updated expected log likelihood			(scalar)
%		
%	_____INPUTS_____________________________________________________________
%	pps	current PPS					(structure)
%	y	unormalized data in original dimension		(col vectors)
%	options	active options					(string)
%		-debug		{0,1}
%		-regularize	regularization factor		(real+)
%				0	default
%				overides pps.gamma
%				(also specified in pps.gamma by lans_ppsinit)
%
%	_____NOTES______________________________________________________________
%	- only R is computed (E-step) for alpha~=1, the M-step is approximated
%	  assuming spherical variance
%
%	_____SEE ALSO___________________________________________________________
%	lans_pps1	lans_ppsgem1	lans_ppspost	lans_ppsinit
%
%	(C) 1999.11.15 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/

function	[pps,Lcomp,Lavg]	= lans_ppsgtm1(pps,y,options)

debug		= paraget('-debug',options,0);
regularize	= paraget('-regularize',options,0);
pps.gamma	= regularize;

N		= size(y,2);
alpha		= pps.alpha;			% clampling factor
beta		= pps.beta;			% spherical variance
p		= pps.p;			% latent basis activations
D		= pps.D;			% data dimensionality
Q		= pps.Q;			% reduced dimensionality
M		= pps.M;			% # nodes
L		= pps.L;			% # latent basis

%__________	compute OLD posterior/responsibility matrix and likelihood
%__________	using old distribution parameters 
Rold	= pps.R;			% M x N

%__________	update pps.W
%	Use a sparse representation for
%	the weight regularizing matrix
LEFT	= zeros(L, L);
cholDcmp= zeros(L, L);
G	= sum(Rold,2);				% M x 1
if (regularize > 0)
	LAMBDA		= regularize*speye(L);
	LAMBDA(L,L)	= 0;
%	LEFT 	= full(p*spdiags(G, 0, M, M)*p' + LAMBDA/beta(1));
	LEFT 	= full(p*spdiags(G, 0, M, M)*p' + LAMBDA);
else
	LEFT	= full(p*spdiags(G, 0, M, M)*p');
end  

% A is a symmetric matrix likely to be positive definite, so try
% fast Cholesky decomposition to calculate WT, otherwise use SVD.
% ????(actT*(Rold*yT)) is computed rigth-to-left, as gtmGlobalR
% ???? and yT are normally (much) larger than actT.
[cholDcmp singular] = chol(LEFT);

if (singular)
	WT = pinv(LEFT)*(p*(Rold*y'));
else
	WT = cholDcmp \ (cholDcmp' \ (p*(Rold*y')));
end
pps.W	= WT';

%__________	update centers pps.f in data space
lf	= pps.W*p;
pps.f	= lans_lin2md(lf,pps.xdim);

%__________	update tangential manifold gradient pps.dfdx if alpha~=1
if (alpha~=1)|debug
	dfdx	= pps.W(:,1:end)*reshape(pps.dpdx,L,M*Q);
	pps.dfdx= reshape(dfdx,D,Q,M);	% NOT orthogonal but spans R^Q
end

%__________	update pps.beta
% assumptions used:
% 	beta is common (scalar)
% 	beta is spherical
pps.beta= N*D / sum(sum(Rold.*lans_dist(lf,y,'-metric Euclidean2')));

%__________	Compute new Responsibility matrix and Lcomp, Lavg
[pps.R,Lcomp,Lavg]	= lans_ppspost(pps,y,options,Rold);