lans_ppsgem1.m

模式识别工具包

%	lans_ppsgem1	- Generalized EM algorithm iteration of PPS
%
%	[pps]	= lans_ppsgem1(pps,y,options)
%
%	_____OUTPUTS____________________________________________________________
%	pps	updated PPS					(structure)
%		
%	_____INPUTS_____________________________________________________________
%	pps	current PPS					(structure)
%	y	unormalized data in original dimension		(col vectors)
%	options	active options					(string)
%		-eta		learning rate for GEM algorithm	(scalar)
%		-iter		MAX # iterations		(scalar)
%		-regularize	regularization factor		(real+)
%				0	default
%		-tol		tolerance factor for GEM	(scalar)
%				0.01	default
%
%	_____NOTES______________________________________________________________
%	- learning rate eta must be SMALL <1e-4
%	- uses MATLAB's optimization toolbox's fsolve.m
%
%	_____SEE ALSO___________________________________________________________
%	ppsgem1	lans_ppsgtm1 lans_pps1	lans_ppspost	lans_ppsinit
%	fsolve
%
%	(C) 1999.10.14 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]	= lans_ppsgem1(pps,y,options)

eta		= paraget('-eta',options,0.0001);
iter		= paraget('-iter',options,200);
mom		= paraget('-mom',options,.1);		%momentum
regularize	= paraget('-regularize',options,0);
tol		= paraget('-tol',options,1e-2);

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 posterior/responsibility matrix and likelihood
[Rold,Lavg,Lcomp]	= lans_ppspost(pps,y,options);  % M x N

%__________	iteratively update pps.W	(unoptimize, computationally)
B	= (D-alpha*Q)./(beta*(D-Q));
S	= alpha./beta;

%	compute constants
coeff1	= 1/B;
coeff2	= (S-B)./(B.*S);
for m=1:M
	Psi(:,:,m)	= pps.dpdx(:,:,m)*pps.dpdx(:,:,m)';			
end

%	iterate
ioptions	= foptions;
ioptions(1)	= 1;
ioptions(2)	= tol;
ioptions(3)	= tol;
ioptions(5)	= 1;
ioptions(14)	= iter;

[W,roptions]	= fsolve('ppsgem',pps.W,ioptions,'',pps,y,Rold,regularize);
pps.W	= W;

%__________	update centers pps.f in data space
pps.f	= pps.W*p;

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

%__________	update pps.beta
coeff1	= (D-Q)/(D-alpha*Q);
coeff2	= (alpha-1)*D/(alpha*(D-alpha*Q));
thesum	= 0;
for n=1:N
	for m=1:M
		err	= y(:,n)-W*p(:,m);
%		t1	= W*Psi(:,:,m)*W';
		ot	= orth(W*pps.dpdx(:,:,m));
		t1	= ot*ot';
		term1	= coeff1*err'*err;
		term2	= coeff2*err'*t1*err;
		thesum	= thesum + Rold(m,n)*(term1-term2);
	end
end
%gpps	= lans_ppsgtm(pps,y,options);
%model_beta	= gpps.beta

beta		= N*D/thesum;
pps.beta	= beta;