mfbox_snmf.m

toolbox for spm 5 for data, model free analysis

function [W,H,sobj,chg]=mfbox_snmf(V,numC,alpha,epsilon,maxiter)
% Copyright by Peter Gruber 
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

error(nargchk(2,5,nargin));

if (nargin<3) alpha = 0.5; end
if (nargin<4) epsilon = 0.00001; end
if (nargin<5) maxiter = 1000; end

[vdim,samples] = size(V);

mV = [];
if (min(V(:))<0)
    mV = min(V,[],2);
    V = V-repmat(mV,1,samples);
end
 
W = abs(randn(vdim,numC));
H = abs(randn(numC,samples));

W = W./(ones(vdim,1)*sum(W,1));

sobj = zeros(1,maxiter+1);
chg = zeros(1,maxiter);

sobj(1) = calcobj(V,W,H,alpha);

for i=1:maxiter
    Wold = W;
    Hold = H;
    
    WH = W*H;
    VWH = V./WH;
    H = (H.*(W'*VWH))/(1+alpha);
    W = W.*(VWH*H')./repmat(sum(H'),vdim,1);
    W = W./repmat(sum(W),vdim,1);
    
    if (nargout>2)
        sobj(i+1) = calcobj(V,W,H,alpha);
    end

    chg(i) = norm(eye(numC)-pinv(Wold)*W,'fro')+norm(eye(numC)-H*pinv(Hold),'fro');
    if (chg(i)<epsilon*numC)
        break
    end
end
if (i==maxiter)
    warning('maxiter reached before threshold %f',chg(i));
end
sobj = sobj(1:i+1);
chg = chg(1:i);

if (~isempty(mV))
    W = W+repmat(mV,1,numC);
end

return


function t=calcobj(V,W,H,alpha)
WH = W*H;
VWH = V./WH;
t = sum(sum((V.*log(VWH)) - V + WH));
t = t + alpha*sum(H(:));