mgrad.m

gaussian procession 原程序,用于解决分类问题,有兴趣可以看看啊

function grad = mgrad(vectheta)

% gradient of the log probability
%
% grad = mgrad(vectheta)
%
% vectheta  : vector of hyperparameters

%            Matlab code for Gaussian Processes for Classification:
%                      GPCLASS version 0.2  10 Nov 97
%       Copyright (c) David Barber and Christopher K I Williams (1997)

%    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
%    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., 675 Mass Ave, Cambridge, MA 02139, USA.

global x_tr ntr ybin x_all jitter m	% general data and cov. mat parameters
global yeta invSxTVEC sinvSxTVEC	% global to start mode search at old points
global vec_mean_prior vec_var_prior;	% these are for the hyperprior distribution

th = 100;				% threshold value before exp
a_thresh = 100;

grad = zeros(1,length(vectheta));

%fprintf(1,'\n Gradient\n'); fprintf(1,'==========\n')

Theta = threshold(vec2mitheta(vectheta,m),th);	% covariance function parameters

noiseproc = getnoise_sp(full(yeta),m);	% initalise the noise process
pye = minvg(yeta);

% set up the covariance matrix:
mntr = m*ntr;
B=zeros(mntr,mntr);
Bno_bias=zeros(mntr,mntr);
o = size(Theta,2);		% 	# parameters in the covariance function
tnp = size(Theta,1)*o;		% 	dimension of gradient vector.

for ci = 1:m;
  D = zeros(ntr,ntr);
  a = threshold(exp(Theta(ci,1)),a_thresh);
  bias = exp(Theta(ci,o));	% bias defined by last component
  for in_mcomp = 2:o-1
    k = in_mcomp;	
    epsi = exp(Theta(ci,k));
    oo = ones(ntr,1); ot = ones(1,ntr);
    xv = x_tr(k-1,:)';      
    xxv = xv.*xv;
    C = xxv*ot + oo*xxv'-2.*xv*xv';
    D = D + epsi.*C;
  end  
  [D,tflag] = threshold(D,th);
  Bno_bias(1+(ci-1)*ntr:ci*ntr,1+(ci-1)*ntr:ci*ntr) = a.*exp(-0.5.*D);  
  B(1+(ci-1)*ntr:ci*ntr,1+(ci-1)*ntr:ci*ntr) =...
    Bno_bias(1+(ci-1)*ntr:ci*ntr,1+(ci-1)*ntr:ci*ntr) + bias.*ones(ntr,ntr);
end

B = B + eye(mntr).*jitter.^2;		% stabilize B
B = sparse(B);

% find the mode

noiseloop=0;
%fprintf(1,'Grad Pot noiseloop: \n');
pyeold = ones(mntr,1);

pye = minvg(yeta);
dim = length(pye);
toln = 0.00000001;

[yeta,pye,noiseproc,sinvSxTVEC,invSxTVEC]=getmode(sinvSxTVEC,yeta,ybin,toln,B,m,ntr);

invS4 = invIpSN(ntr,m,pye,B);
invS5 = noiseproc*invS4;
ybp = ybin-pye;

mpye = full(reshape(pye,ntr,m));

FB = full(Bno_bias);
SB = sparse(B);
invS4xB = invS4*SB;			% bit expensive
noiseprocxinvS4 = noiseproc*invS4;
Pmat = getmPmat(pye,ntr);
FP = invS4xB*Pmat;
diagF = diag(invS4xB);

for cm = 1:m
  sbnoiseprocxinvS4  = noiseprocxinvS4(1+(cm-1)*ntr:cm*ntr,1+(cm-1)*ntr:cm*ntr);
  for in_mcomp = 1:o
    k = in_mcomp;	% 		which component in m
    dB = zeros(mntr,mntr);
    dnB = zeros(ntr,ntr);        
    if in_mcomp == 1
      dnB = Bno_bias(1+(cm-1)*ntr:cm*ntr,1+(cm-1)*ntr:cm*ntr);
      dB(1+(cm-1)*ntr:cm*ntr,1+(cm-1)*ntr:cm*ntr)=dnB;
    elseif in_mcomp == o;
      bias = exp(Theta(cm,o));	% 	bias defined by last component
      dnB = ones(ntr,ntr).*bias;
      dB(1+(cm-1)*ntr:cm*ntr,1+(cm-1)*ntr:cm*ntr) =dnB;
    else      
      hw = -0.5.*exp(Theta(cm,k));
      oo = ones(ntr,1); ot = ones(1,ntr);
      xv = x_tr(k-1,:)'; xxv = xv.*xv;
      C = xxv*ot + oo*xxv'-2.*xv*xv';
      dnB = hw.*FB(1+(cm-1)*ntr:cm*ntr,1+(cm-1)*ntr:cm*ntr).*C;
      dB(1+(cm-1)*ntr:cm*ntr,1+(cm-1)*ntr:cm*ntr) = dnB;
    end    
    dB=sparse(dB);
    A=invS5*(dB*ybp);
    mA = full(reshape(A,ntr,m));
    Amat = getmPmat(A,ntr);
    dN = zeros(mntr,mntr); dN = diag(A.*(1-2*pye));
    for ci = 1:m;
      for cj = 1:ci;
	if ci~=cj
	  dN(1+ntr*(ci-1):ntr*ci,1+ntr*(cj-1):ntr*cj)=-diag(mpye(:,ci).*mA(:,cj) + mpye(:,cj).*mA(:,ci));
	  dN(1+ntr*(cj-1):ntr*cj,1+ntr*(ci-1):ntr*ci)=-diag(mpye(:,ci).*mA(:,cj) + mpye(:,cj).*mA(:,ci));     
	end
      end
    end  
   
    comp = in_mcomp + (cm-1)*o;

    dN = sparse(dN); noiseproc = sparse(noiseproc);
    grad(comp)=ybp'*dB*ybp;
    tr(comp)=diagF'*A - 2.*trFPA(FP,A,ntr,m)+trAB(sbnoiseprocxinvS4,dnB);
  end					% end of component loop
end

grad = 0.5.*(-grad+tr) - mgradlnptheta(vectheta);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function g = mgradlnptheta(theta)

global vec_mean_prior vec_var_prior

g = -(theta-vec_mean_prior)./vec_var_prior;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function t = trFPA(FP,A,nn,m)

t=0;
for c = 1:m
  d = diag(FP(1+(c-1)*nn:c*nn,1:nn));
  fe = A(1+(c-1)*nn:c*nn,1);
  t = t + d'*fe;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function P = getmPmat(pye,nn)

m = size(pye,1)./nn;

P = zeros(m*nn,nn);

for c = 1:m
  P(1+(c-1)*nn:c*nn,1:nn) = getPmat(c,pye,nn);
end

P = sparse(P);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function t= trAB(A,B)

A = full(A); B = full(B');

ts = size(A,1)*size(A,2);
a = reshape(A,ts,1);
b = reshape(B,ts,1);

t = a'*b;