mvrvm.rar_MVRVM_relevance

% @file mvrvm.m % Author Arasanathan Thayananthan ( at315@cam.ac.uk) % (c) Copyright University of Cambridge % % This library is free software; you can redistribute it and/or % modify it under the terms of the GNU Lesser General Public % License as published by the Free Software Foundation; either % version 2 of the License, or (at your option) any later version. % % This library 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 % Lesser General Public License for more details. % % You should have received a copy of the GNU Lesser General Public % License along with this library; if not, write to the Free Software % Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA %This m file implements the bottom-up relevance vector machine with multivariate outputs, proposed %in % Multivariate Relevance Vector Machines for Tracking % Arasanathan Thayananthan et al. (University of Cambridge) % in Proc. 9th European Conference on Computer Vision 2006. %Detailed derivation can be found in %Arasanathan Thayananthan. Template-based pose estimation and tracking of 3 hand motion %PhD thesis, University of Cambridge, UK 2005. % output % weights trained weights % used index into column of PHI of relevant basis % alpha final set of hyperparamters % beta final values of noise precisions. % input % PHI basis function/design matrix % t targetmatrix (multivariate) % N number of training examples % M number of basis functions % [ N M] =size(PHI) % P the number of output dimensions % [ N P]= size(t) % maxIts number of iterations function [weights, used, alpha, beta, SIGMA_nz, N] = mvrvm(PHI,t,maxIts) ALPHA_MAX = 1e12; %initialise beta epsilon = std(t) * 10/100; beta=1./(epsilon.^2); [N,M] = size(PHI); [N,P] = size(t); w = zeros(M,P); PHIt = PHI'*t; alpha = ones(M,1); nonZero = logical(ones(M,1)); alpha(:)=ALPHA_MAX; %choosing the first basis function with non-zero weight max_change=-inf; max_k=-1; max_alpha=-inf; for k=1:M phi=PHI(:,k); for j=1:P S(j)=beta(j)*phi'*phi; Q(j)=beta(j)*phi'*t(:,j); end if(alpha(k)<ALPHA_MAX) s=alpha(k)*S./(alpha(k)-S); q=alpha(k)*Q./(alpha(k)-S); else s=S; q=Q; end [new_alpha,l_inc]=SolveForAlpha(s',q',S',Q',alpha(k)); l_inc_vec(k)=l_inc; alpha_vec(k)=new_alpha; if(max_change<l_inc) max_k=k; max_s=s; max_q=q; max_alpha=new_alpha; max_change=l_inc; end end alpha(max_k)=max_alpha; %start the iterations for i=1:maxIts i; nonZero = (alpha<ALPHA_MAX); alpha_nz = alpha(nonZero); w(~nonZero) = 0; nz=size(alpha_nz,1); PHI_nz = PHI(:,nonZero); SIGMA_nz=zeros(nz,nz,P); for j=1:P SIGMA_nz(:,:,j) = inv((PHI_nz'*PHI_nz)*beta(j) + diag(alpha_nz)); end %choose the basis function which minimises the marginal likelihood max_change=-inf; max_k=-1; max_alpha=-inf; for k=1:M phi=PHI(:,k); for j=1:P S(j)=beta(j)*phi'*phi-beta(j)*beta(j)*phi'*PHI_nz*SIGMA_nz(:,:,j)*PHI_nz'*phi; Q(j)=beta(j)*phi'*t(:,j) -beta(j)*beta(j)*phi'*PHI_nz*SIGMA_nz(:,:,j)*PHI_nz'*t(:,j); end if(alpha(k)<ALPHA_MAX) s=alpha(k)*S./(alpha(k)-S); q=alpha(k)*Q./(alpha(k)-S); else s=S; q=Q; end [new_alpha,l_inc]=SolveForAlpha(s',q',S',Q',alpha(k)); l_inc_vec(k)=l_inc; alpha_vec(k)=new_alpha; if(max_change<l_inc) max_k=k; max_s=s; max_q=q; max_alpha=new_alpha; max_change=l_inc; end end alpha(max_k)=max_alpha; if(nz~=0) d=zeros(nz,1); for j=1:P d=d+diag(SIGMA_nz(:,:,j)); end gamma=P-alpha_nz.*d; for j=1:P gamma=1-alpha_nz.*diag(SIGMA_nz(:,:,j)); weights=beta(j)*SIGMA_nz(:,:,j)*PHI_nz'*t(:,j); ED=sum((t(:,j)-PHI_nz*weights).^2); beta(j)=(N-sum(gamma))/ED; end end end weights=zeros(nz,P); for j=1:P weights(:,j)=beta(j)*SIGMA_nz(:,:,j)*PHI_nz'*t(:,j); end used = find(nonZero); function [new_alpha,l_inc]=SolveForAlpha(s,q,S,Q,old_alpha) ALPHA_MAX = 1e12; [n m]=size(s); index=[1:n]'; C=zeros(n,1); CC=zeros(2*n-1,1); for i=1:n SS=-s(find(index~=i)); P=poly(SS)'; PP=conv(P,P); CC=CC+q(i)*q(i)*PP; C=C+P; end P=poly(-s); PP=n*conv(P,P)'; CC0=conv(P,C)'; CC=[0 CC]; CCC=[CC+CC0 0]; r=roots(PP-CCC); r=r(find(abs(imag(r))==0)); r=r(find(r>0)); r=[r ALPHA_MAX]; [nn mm]=size(r); L=zeros(nn,1); for i=1:nn new_alpha=r(i); if((old_alpha>=ALPHA_MAX)& (r(i)<ALPHA_MAX)) tt(i)=1; for j=1:n L(i)=L(i)+ log(new_alpha/(new_alpha+s(j)))+((q(j)*q(j)./(new_alpha+s(j)))); end elseif((old_alpha<ALPHA_MAX)&(new_alpha>=ALPHA_MAX)) tt(i)=2; for j=1:n L(i)=L(i)+ (Q(j)*Q(j)/(S(j)-old_alpha))-log(1- (S(j)/old_alpha)); end elseif((old_alpha<ALPHA_MAX)&(new_alpha<ALPHA_MAX)) tt(i)=3; for j=1:n ad=(1/new_alpha)-(1/old_alpha); L(i)=L(i)+ (Q(j)*Q(j)/(S(j)+(1/ad))) -log(1+S(j)*ad); end else tt(i)=4; L(i)=0; end end [m ind]=max(L); new_alpha=r(ind); l_inc=L(ind);