--- /dev/null
+## Copyright (c) 2012 by Jingu Kim and Haesun Park <jingu@cc.gatech.edu>
+##
+## 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 3 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, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{W}, @var{H}, @var{iter}, @var{HIS}] = } nmf_bpas (@var{A}, @var{k})
+## Nonnegative Matrix Factorization by Alternating Nonnegativity Constrained Least Squares
+## using Block Principal Pivoting/Active Set method.
+##
+## This function solves one the following problems: given @var{A} and @var{k}, find @var{W} and @var{H} such that
+## (1) minimize 1/2 * || @var{A}-@var{W}@var{H} ||_F^2
+## (2) minimize 1/2 * ( || @var{A}-@var{W}@var{H} ||_F^2 + alpha * || @var{W} ||_F^2 + beta * || @var{H} ||_F^2 )
+## (3) minimize 1/2 * ( || @var{A}-@var{W}@var{H} ||_F^2 + alpha * || @var{W} ||_F^2 + beta * (sum_(i=1)^n || @var{H}(:,i) ||_1^2 ) )
+## where @var{W}>=0 and @var{H}>=0 elementwise.
+## The input arguments are @var{A} : Input data matrix (m x n) and @var{k} : Target low-rank.
+##
+##
+## @strong{Optional Inputs}
+## @table @samp
+## @item Type : Default is 'regularized', which is recommended for quick application testing unless 'sparse' or 'plain' is explicitly needed. If sparsity is needed for 'W' factor, then apply this function for the transpose of 'A' with formulation (3). Then, exchange 'W' and 'H' and obtain the transpose of them. Imposing sparsity for both factors is not recommended and thus not included in this software.
+## @table @asis
+## @item 'plain' to use formulation (1)
+## @item 'regularized' to use formulation (2)
+## @item 'sparse' to use formulation (3)
+## @end table
+##
+## @item NNLSSolver : Default is 'bp', which is in general faster.
+## @table @asis
+## item 'bp' to use the algorithm in [1]
+## item 'as' to use the algorithm in [2]
+## @end table
+##
+## @item Alpha : Parameter alpha in the formulation (2) or (3). Default is the average of all elements in A. No good justfication for this default value, and you might want to try other values.
+## @item Beta : Parameter beta in the formulation (2) or (3).
+## Default is the average of all elements in A. No good justfication for this default value, and you might want to try other values.
+## @item MaxIter : Maximum number of iterations. Default is 100.
+## @item MinIter : Minimum number of iterations. Default is 20.
+## @item MaxTime : Maximum amount of time in seconds. Default is 100,000.
+## @item Winit : (m x k) initial value for W.
+## @item Hinit : (k x n) initial value for H.
+## @item Tol : Stopping tolerance. Default is 1e-3. If you want to obtain a more accurate solution, decrease TOL and increase MAX_ITER at the same time.
+## @item Verbose :
+## @table @asis
+## @item 0 (default) - No debugging information is collected.@*
+## @item 1 (debugging purpose) - History of computation is returned by 'HIS' variable.
+## @item 2 (debugging purpose) - History of computation is additionally printed on screen.
+## @end table
+## @end table
+##
+## @strong{Outputs}
+## @table @samp
+## @item W : Obtained basis matrix (m x k)
+## @item H : Obtained coefficients matrix (k x n)
+## @item iter : Number of iterations
+## @item HIS : (debugging purpose) History of computation
+## @end table
+##
+## Usage Examples:
+## @example
+## nmf(A,10)
+## nmf(A,20,'verbose',2)
+## nmf(A,30,'verbose',2,'nnls_solver','as')
+## nmf(A,5,'verbose',2,'type','sparse')
+## nmf(A,60,'verbose',1,'type','plain','w_init',rand(m,k))
+## nmf(A,70,'verbose',2,'type','sparse','nnls_solver','bp','alpha',1.1,'beta',1.3)
+## @end example
+##
+## References:
+## [1] For using this software, please cite:@*
+## Jingu Kim and Haesun Park, Toward Faster Nonnegative Matrix Factorization: A New Algorithm and Comparisons,@*
+## In Proceedings of the 2008 Eighth IEEE International Conference on Data Mining (ICDM'08), 353-362, 2008@*
+## [2] If you use 'nnls_solver'='as' (see below), please cite:@*
+## Hyunsoo Kim and Haesun Park, Nonnegative Matrix Factorization Based @*
+## on Alternating Nonnegativity Constrained Least Squares and Active Set Method, @*
+## SIAM Journal on Matrix Analysis and Applications, 2008, 30, 713-730
+##
+## Check original code at @url{http://www.cc.gatech.edu/~jingu}
+##
+## @seealso{nmf_pg}
+## @end deftypefn
+
+## 2012 - Modified and adapted to Octave 3.6.1 by
+## Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
+
+# TODO
+# - Format code.
+# - Vectorize loops.
+
+function [W, H, iter, HIS] = nmf_bpas (A, k , varargin)
+ page_screen_output (0, "local");
+
+ [m,n] = size(A);
+ ST_RULE = 1;
+
+ # --- Parse arguments --- #
+ parser = inputParser ();
+ parser.FunctionName = "nmf_bpas";
+ parser = addParamValue (parser,'Winit', rand(m,k), @ismatrix);
+ parser = addParamValue (parser,'Hinit', rand(k,n), @ismatrix);
+ parser = addParamValue (parser,'Tol', 1e-3, @(x)x>0);
+ parser = addParamValue (parser,'Alpha', mean (A(:)), @(x)x>=0);
+ parser = addParamValue (parser,'Beta', mean (A(:)), @(x)x>=0);
+ parser = addParamValue (parser,'MaxIter', 100, @(x)x>0);
+ parser = addParamValue (parser,'MaxTime', 1e3, @(x)x>0);
+ parser = addParamValue (parser,'Verbose', false);
+
+ val_type = @(x,c) ischar (x) && any (strcmpi (x,c));
+ parser = addParamValue (parser,'Type', 'regularized', ...
+ @(x)val_type (x,{'regularized', 'sparse','plain'}));
+ parser = addParamValue (parser,'NNLSSolver', 'bp', ...
+ @(x)val_type (x,{'bp', 'as'}));
+
+ parser = parse(parser,varargin{:});
+
+ % Default configuration
+ par.m = m;
+ par.n = n;
+ par.type = parser.Results.Type;
+ par.nnls_solver = parser.Results.NNLSSolver;
+ par.alpha = parser.Results.Alpha;
+ par.beta = parser.Results.Beta;
+ par.max_iter = parser.Results.MaxIter;
+ par.min_iter = 20;
+ par.max_time = parser.Results.MaxTime;
+ par.tol = parser.Results.Tol;
+ par.verbose = parser.Results.Verbose;
+ W = parser.Results.Winit;
+ H = parser.Results.Hinit;
+
+ # TODO check if can be removed
+ argAlpha = par.alpha;
+ argBeta = par.beta;
+
+ clear parser val_type
+
+### PARSING TYPE
+# TODO add callbacks here to use during main loop. See [1]
+ % for regularized/sparse case
+ salphaI = sqrt (par.alpha) * eye (k);
+ zerokm = zeros (k,m);
+
+ if strcmpi (par.type, 'regularized')
+ sbetaI = sqrt (par.beta) * eye (k);
+ zerokn = zeros (k,n);
+
+ elseif strcmpi (par.type, 'sparse')
+ sbetaE = sqrt (par.beta) * ones (1,k);
+ betaI = par.beta * ones (k,k);
+ zero1n = zeros (1,n);
+
+ end
+###
+
+# Verbosity
+ display(par);
+### Done till here Sun Mar 25 19:00:26 2012
+
+ HIS = 0;
+ if par.verbose % collect information for analysis/debugging
+ [gradW,gradH] = getGradient(A,W,H,par.type,par.alpha,par.beta);
+ initGrNormW = norm(gradW,'fro');
+ initGrNormH = norm(gradH,'fro');
+ initNorm = norm(A,'fro');
+ numSC = 3;
+ initSCs = zeros(numSC,1);
+ for j=1:numSC
+ initSCs(j) = getInitCriterion(j,A,W,H,par.type,par.alpha,par.beta,gradW,gradH);
+ end
+%---(1)------(2)--------(3)--------(4)--------(5)---------(6)----------(7)------(8)-----(9)-------(10)--------------(11)-------
+% iter # | elapsed | totalTime | subIterW | subIterH | rel. obj.(%) | NM_GRAD | GRAD | DELTA | W density (%) | H density (%)
+%------------------------------------------------------------------------------------------------------------------------------
+ HIS = zeros(1,11);
+ HIS(1,[1:5])=0;
+ ver.initGrNormW = initGrNormW;
+ ver.initGrNormH = initGrNormH;
+ ver.initNorm = initNorm; HIS(1,6) = ver.initNorm;
+ ver.SC1 = initSCs(1); HIS(1,7) = ver.SC1;
+ ver.SC2 = initSCs(2); HIS(1,8) = ver.SC2;
+ ver.SC3 = initSCs(3); HIS(1,9) = ver.SC3;
+ ver.W_density = length(find(W>0))/(m*k); HIS(1,10) = ver.W_density;
+ ver.H_density = length(find(H>0))/(n*k); HIS(1,11) = ver.H_density;
+ if par.verbose == 2
+ disp (ver);
+ end
+ tPrev = cputime;
+ end
+
+ tStart = cputime;
+ tTotal = 0;
+ initSC = getInitCriterion(ST_RULE,A,W,H,par.type,par.alpha,par.beta);
+ SCconv = 0;
+ SC_COUNT = 3;
+
+#TODO: [1] Replace with callbacks avoid switching each time
+ for iter=1:par.max_iter
+ switch par.type
+ case 'plain'
+ [H,gradHX,subIterH] = nnlsm(W,A,H,par.nnls_solver);
+ [W,gradW,subIterW] = nnlsm(H',A',W',par.nnls_solver);, W=W';, gradW=gradW';
+ gradH = (W'*W)*H - W'*A;
+ case 'regularized'
+ [H,gradHX,subIterH] = nnlsm([W;sbetaI],[A;zerokn],H,par.nnls_solver);
+ [W,gradW,subIterW] = nnlsm([H';salphaI],[A';zerokm],W',par.nnls_solver);, W=W';, gradW=gradW';
+ gradH = (W'*W)*H - W'*A + par.beta*H;
+ case 'sparse'
+ [H,gradHX,subIterH] = nnlsm([W;sbetaE],[A;zero1n],H,par.nnls_solver);
+ [W,gradW,subIterW] = nnlsm([H';salphaI],[A';zerokm],W',par.nnls_solver);, W=W';, gradW=gradW';
+ gradH = (W'*W)*H - W'*A + betaI*H;
+ end
+
+ if par.verbose % collect information for analysis/debugging
+ elapsed = cputime-tPrev;
+ tTotal = tTotal + elapsed;
+ ver = [];
+ idx = iter+1;
+%---(1)------(2)--------(3)--------(4)--------(5)---------(6)----------(7)------(8)-----(9)-------(10)--------------(11)-------
+% iter # | elapsed | totalTime | subIterW | subIterH | rel. obj.(%) | NM_GRAD | GRAD | DELTA | W density (%) | H density (%)
+%------------------------------------------------------------------------------------------------------------------------------
+ ver.iter = iter; HIS(idx,1)=iter;
+ ver.elapsed = elapsed; HIS(idx,2)=elapsed;
+ ver.tTotal = tTotal; HIS(idx,3)=tTotal;
+ ver.subIterW = subIterW; HIS(idx,4)=subIterW;
+ ver.subIterH = subIterH; HIS(idx,5)=subIterH;
+ ver.relError = norm(A-W*H,'fro')/initNorm; HIS(idx,6)=ver.relError;
+ ver.SC1 = getStopCriterion(1,A,W,H,par.type,par.alpha,par.beta,gradW,gradH)/initSCs(1); HIS(idx,7)=ver.SC1;
+ ver.SC2 = getStopCriterion(2,A,W,H,par.type,par.alpha,par.beta,gradW,gradH)/initSCs(2); HIS(idx,8)=ver.SC2;
+ ver.SC3 = getStopCriterion(3,A,W,H,par.type,par.alpha,par.beta,gradW,gradH)/initSCs(3); HIS(idx,9)=ver.SC3;
+ ver.W_density = length(find(W>0))/(m*k); HIS(idx,10)=ver.W_density;
+ ver.H_density = length(find(H>0))/(n*k); HIS(idx,11)=ver.H_density;
+ if par.verbose == 2, display(ver);, end
+ tPrev = cputime;
+ end
+
+ if (iter > par.min_iter)
+ SC = getStopCriterion(ST_RULE,A,W,H,par.type,par.alpha,par.beta,gradW,gradH);
+ if (par.verbose && (tTotal > par.max_time)) || (~par.verbose && ((cputime-tStart)>par.max_time))
+ break;
+ elseif (SC/initSC <= par.tol)
+ SCconv = SCconv + 1;
+ if (SCconv >= SC_COUNT)
+ break;
+ end
+ else
+ SCconv = 0;
+ end
+ end
+ end
+ [m,n]=size(A);
+ norm2=sqrt(sum(W.^2,1));
+ toNormalize = norm2>0;
+ W(:,toNormalize) = W(:,toNormalize)./repmat(norm2(toNormalize),m,1);
+ H(toNormalize,:) = H(toNormalize,:).*repmat(norm2(toNormalize)',1,n);
+
+ final.iterations = iter;
+ if par.verbose
+ final.elapsed_total = tTotal;
+ else
+ final.elapsed_total = cputime-tStart;
+ end
+ final.relative_error = norm(A-W*H,'fro')/norm(A,'fro');
+ final.W_density = length(find(W>0))/(m*k);
+ final.H_density = length(find(H>0))/(n*k);
+ display(final);
+
+endfunction
+
+%------------------------------------------------------------------------------------------------------------------------
+% Utility Functions
+%-------------------------------------------------------------------------------
+function retVal = getInitCriterion(stopRule,A,W,H,type,alpha,beta,gradW,gradH)
+% STOPPING_RULE : 1 - Normalized proj. gradient
+% 2 - Proj. gradient
+% 3 - Delta by H. Kim
+% 0 - None (want to stop by MAX_ITER or MAX_TIME)
+ if nargin~=9
+ [gradW,gradH] = getGradient(A,W,H,type,alpha,beta);
+ end
+ [m,k]=size(W);, [k,n]=size(H);, numAll=(m*k)+(k*n);
+ switch stopRule
+ case 1
+ retVal = norm([gradW; gradH'],'fro')/numAll;
+ case 2
+ retVal = norm([gradW; gradH'],'fro');
+ case 3
+ retVal = getStopCriterion(3,A,W,H,type,alpha,beta,gradW,gradH);
+ case 0
+ retVal = 1;
+ end
+
+endfunction
+%-------------------------------------------------------------------------------
+function retVal = getStopCriterion(stopRule,A,W,H,type,alpha,beta,gradW,gradH)
+% STOPPING_RULE : 1 - Normalized proj. gradient
+% 2 - Proj. gradient
+% 3 - Delta by H. Kim
+% 0 - None (want to stop by MAX_ITER or MAX_TIME)
+ if nargin~=9
+ [gradW,gradH] = getGradient(A,W,H,type,alpha,beta);
+ end
+
+ switch stopRule
+ case 1
+ pGradW = gradW(gradW<0|W>0);
+ pGradH = gradH(gradH<0|H>0);
+ pGrad = [gradW(gradW<0|W>0); gradH(gradH<0|H>0)];
+ pGradNorm = norm(pGrad);
+ retVal = pGradNorm/length(pGrad);
+ case 2
+ pGradW = gradW(gradW<0|W>0);
+ pGradH = gradH(gradH<0|H>0);
+ pGrad = [gradW(gradW<0|W>0); gradH(gradH<0|H>0)];
+ retVal = norm(pGrad);
+ case 3
+ resmat=min(H,gradH); resvec=resmat(:);
+ resmat=min(W,gradW); resvec=[resvec; resmat(:)];
+ deltao=norm(resvec,1); %L1-norm
+ num_notconv=length(find(abs(resvec)>0));
+ retVal=deltao/num_notconv;
+ case 0
+ retVal = 1e100;
+ end
+
+endfunction
+%-------------------------------------------------------------------------------
+function [gradW,gradH] = getGradient(A,W,H,type,alpha,beta)
+ switch type
+ case 'plain'
+ gradW = W*(H*H') - A*H';
+ gradH = (W'*W)*H - W'*A;
+ case 'regularized'
+ gradW = W*(H*H') - A*H' + alpha*W;
+ gradH = (W'*W)*H - W'*A + beta*H;
+ case 'sparse'
+ k=size(W,2);
+ betaI = beta*ones(k,k);
+ gradW = W*(H*H') - A*H' + alpha*W;
+ gradH = (W'*W)*H - W'*A + betaI*H;
+ end
+
+endfunction
+
+%------------------------------------------------------------------------------------------------------------------------
+function [X,grad,iter] = nnlsm(A,B,init,solver)
+ switch solver
+ case 'bp'
+ [X,grad,iter] = nnlsm_blockpivot(A,B,0,init);
+ case 'as'
+ [X,grad,iter] = nnlsm_activeset(A,B,1,0,init);
+ end
+
+endfunction
+%------------------------------------------------------------------------------------------------------------------------
+function [ X,Y,iter,success ] = nnlsm_activeset( A, B, overwrite, isInputProd, init)
+% Nonnegativity Constrained Least Squares with Multiple Righthand Sides
+% using Active Set method
+%
+% This software solves the following problem: given A and B, find X such that
+% minimize || AX-B ||_F^2 where X>=0 elementwise.
+%
+% Reference:
+% Charles L. Lawson and Richard J. Hanson, Solving Least Squares Problems,
+% Society for Industrial and Applied Mathematics, 1995
+% M. H. Van Benthem and M. R. Keenan,
+% Fast Algorithm for the Solution of Large-scale Non-negativity-constrained Least Squares Problems,
+% J. Chemometrics 2004; 18: 441-450
+%
+% Written by Jingu Kim (jingu@cc.gatech.edu)
+% School of Computational Science and Engineering,
+% Georgia Institute of Technology
+%
+% Last updated Feb-20-2010
+%
+% <Inputs>
+% A : input matrix (m x n) (by default), or A'*A (n x n) if isInputProd==1
+% B : input matrix (m x k) (by default), or A'*B (n x k) if isInputProd==1
+% overwrite : (optional, default:0) if turned on, unconstrained least squares solution is computed in the beginning
+% isInputProd : (optional, default:0) if turned on, use (A'*A,A'*B) as input instead of (A,B)
+% init : (optional) initial value for X
+% <Outputs>
+% X : the solution (n x k)
+% Y : A'*A*X - A'*B where X is the solution (n x k)
+% iter : number of iterations
+% success : 1 for success, 0 for failure.
+% Failure could only happen on a numericall very ill-conditioned problem.
+
+ if nargin<3, overwrite=0;, end
+ if nargin<4, isInputProd=0;, end
+
+ if isInputProd
+ AtA=A;,AtB=B;
+ else
+ AtA=A'*A;, AtB=A'*B;
+ end
+
+ [n,k]=size(AtB);
+ MAX_ITER = n*5;
+ % set initial feasible solution
+ if overwrite
+ [X,iter] = solveNormalEqComb(AtA,AtB);
+ PassSet = (X > 0);
+ NotOptSet = any(X<0);
+ else
+ if nargin<5
+ X = zeros(n,k);
+ PassSet = false(n,k);
+ NotOptSet = true(1,k);
+ else
+ X = init;
+ PassSet = (X > 0);
+ NotOptSet = any(X<0);
+ end
+ iter = 0;
+ end
+
+ Y = zeros(n,k);
+ Y(:,~NotOptSet)=AtA*X(:,~NotOptSet) - AtB(:,~NotOptSet);
+ NotOptCols = find(NotOptSet);
+
+ bigIter = 0;, success=1;
+ while(~isempty(NotOptCols))
+ bigIter = bigIter+1;
+ if ((MAX_ITER >0) && (bigIter > MAX_ITER)) % set max_iter for ill-conditioned (numerically unstable) case
+ success = 0;, bigIter, break
+ end
+
+ % find unconstrained LS solution for the passive set
+ Z = zeros(n,length(NotOptCols));
+ [ Z,subiter ] = solveNormalEqComb(AtA,AtB(:,NotOptCols),PassSet(:,NotOptCols));
+ iter = iter + subiter;
+ %Z(abs(Z)<1e-12) = 0; % One can uncomment this line for numerical stability.
+ InfeaSubSet = Z < 0;
+ InfeaSubCols = find(any(InfeaSubSet));
+ FeaSubCols = find(all(~InfeaSubSet));
+
+ if ~isempty(InfeaSubCols) % for infeasible cols
+ ZInfea = Z(:,InfeaSubCols);
+ InfeaCols = NotOptCols(InfeaSubCols);
+ Alpha = zeros(n,length(InfeaSubCols));, Alpha(:,:) = Inf;
+ InfeaSubSet(:,InfeaSubCols);
+ [i,j] = find(InfeaSubSet(:,InfeaSubCols));
+ InfeaSubIx = sub2ind(size(Alpha),i,j);
+ if length(InfeaCols) == 1
+ InfeaIx = sub2ind([n,k],i,InfeaCols * ones(length(j),1));
+ else
+ InfeaIx = sub2ind([n,k],i,InfeaCols(j)');
+ end
+ Alpha(InfeaSubIx) = X(InfeaIx)./(X(InfeaIx)-ZInfea(InfeaSubIx));
+
+ [minVal,minIx] = min(Alpha);
+ Alpha(:,:) = repmat(minVal,n,1);
+ X(:,InfeaCols) = X(:,InfeaCols)+Alpha.*(ZInfea-X(:,InfeaCols));
+ IxToActive = sub2ind([n,k],minIx,InfeaCols);
+ X(IxToActive) = 0;
+ PassSet(IxToActive) = false;
+ end
+ if ~isempty(FeaSubCols) % for feasible cols
+ FeaCols = NotOptCols(FeaSubCols);
+ X(:,FeaCols) = Z(:,FeaSubCols);
+ Y(:,FeaCols) = AtA * X(:,FeaCols) - AtB(:,FeaCols);
+ %Y( abs(Y)<1e-12 ) = 0; % One can uncomment this line for numerical stability.
+
+ NotOptSubSet = (Y(:,FeaCols) < 0) & ~PassSet(:,FeaCols);
+ NewOptCols = FeaCols(all(~NotOptSubSet));
+ UpdateNotOptCols = FeaCols(any(NotOptSubSet));
+ if ~isempty(UpdateNotOptCols)
+ [minVal,minIx] = min(Y(:,UpdateNotOptCols).*~PassSet(:,UpdateNotOptCols));
+ PassSet(sub2ind([n,k],minIx,UpdateNotOptCols)) = true;
+ end
+ NotOptSet(NewOptCols) = false;
+ NotOptCols = find(NotOptSet);
+ end
+ end
+
+endfunction
+%------------------------------------------------------------------------------------------------------------------------
+function [ X,Y,iter,success ] = nnlsm_blockpivot( A, B, isInputProd, init )
+% Nonnegativity Constrained Least Squares with Multiple Righthand Sides
+% using Block Principal Pivoting method
+%
+% This software solves the following problem: given A and B, find X such that
+% minimize || AX-B ||_F^2 where X>=0 elementwise.
+%
+% Reference:
+% Jingu Kim and Haesun Park, Toward Faster Nonnegative Matrix Factorization: A New Algorithm and Comparisons,
+% In Proceedings of the 2008 Eighth IEEE International Conference on Data Mining (ICDM'08), 353-362, 2008
+%
+% Written by Jingu Kim (jingu@cc.gatech.edu)
+% Copyright 2008-2009 by Jingu Kim and Haesun Park,
+% School of Computational Science and Engineering,
+% Georgia Institute of Technology
+%
+% Check updated code at http://www.cc.gatech.edu/~jingu
+% Please send bug reports, comments, or questions to Jingu Kim.
+% This code comes with no guarantee or warranty of any kind. Note that this algorithm assumes that the
+% input matrix A has full column rank.
+%
+% Last modified Feb-20-2009
+%
+% <Inputs>
+% A : input matrix (m x n) (by default), or A'*A (n x n) if isInputProd==1
+% B : input matrix (m x k) (by default), or A'*B (n x k) if isInputProd==1
+% isInputProd : (optional, default:0) if turned on, use (A'*A,A'*B) as input instead of (A,B)
+% init : (optional) initial value for X
+% <Outputs>
+% X : the solution (n x k)
+% Y : A'*A*X - A'*B where X is the solution (n x k)
+% iter : number of iterations
+% success : 1 for success, 0 for failure.
+% Failure could only happen on a numericall very ill-conditioned problem.
+
+ if nargin<3, isInputProd=0;, end
+ if isInputProd
+ AtA = A;, AtB = B;
+ else
+ AtA = A'*A;, AtB = A'*B;
+ end
+
+ [n,k]=size(AtB);
+ MAX_ITER = n*5;
+ % set initial feasible solution
+ X = zeros(n,k);
+ if nargin<4
+ Y = - AtB;
+ PassiveSet = false(n,k);
+ iter = 0;
+ else
+ PassiveSet = (init > 0);
+ [ X,iter ] = solveNormalEqComb(AtA,AtB,PassiveSet);
+ Y = AtA * X - AtB;
+ end
+ % parameters
+ pbar = 3;
+ P = zeros(1,k);, P(:) = pbar;
+ Ninf = zeros(1,k);, Ninf(:) = n+1;
+ iter = 0;
+
+ NonOptSet = (Y < 0) & ~PassiveSet;
+ InfeaSet = (X < 0) & PassiveSet;
+ NotGood = sum(NonOptSet)+sum(InfeaSet);
+ NotOptCols = NotGood > 0;
+
+ bigIter = 0;, success=1;
+ while(~isempty(find(NotOptCols)))
+ bigIter = bigIter+1;
+ if ((MAX_ITER >0) && (bigIter > MAX_ITER)) % set max_iter for ill-conditioned (numerically unstable) case
+ success = 0;, break
+ end
+
+ Cols1 = NotOptCols & (NotGood < Ninf);
+ Cols2 = NotOptCols & (NotGood >= Ninf) & (P >= 1);
+ Cols3Ix = find(NotOptCols & ~Cols1 & ~Cols2);
+ if ~isempty(find(Cols1))
+ P(Cols1) = pbar;,Ninf(Cols1) = NotGood(Cols1);
+ PassiveSet(NonOptSet & repmat(Cols1,n,1)) = true;
+ PassiveSet(InfeaSet & repmat(Cols1,n,1)) = false;
+ end
+ if ~isempty(find(Cols2))
+ P(Cols2) = P(Cols2)-1;
+ PassiveSet(NonOptSet & repmat(Cols2,n,1)) = true;
+ PassiveSet(InfeaSet & repmat(Cols2,n,1)) = false;
+ end
+ if ~isempty(Cols3Ix)
+ for i=1:length(Cols3Ix)
+ Ix = Cols3Ix(i);
+ toChange = max(find( NonOptSet(:,Ix)|InfeaSet(:,Ix) ));
+ if PassiveSet(toChange,Ix)
+ PassiveSet(toChange,Ix)=false;
+ else
+ PassiveSet(toChange,Ix)=true;
+ end
+ end
+ end
+ NotOptMask = repmat(NotOptCols,n,1);
+ [ X(:,NotOptCols),subiter ] = solveNormalEqComb(AtA,AtB(:,NotOptCols),PassiveSet(:,NotOptCols));
+ iter = iter + subiter;
+ X(abs(X)<1e-12) = 0; % for numerical stability
+ Y(:,NotOptCols) = AtA * X(:,NotOptCols) - AtB(:,NotOptCols);
+ Y(abs(Y)<1e-12) = 0; % for numerical stability
+
+ % check optimality
+ NonOptSet = NotOptMask & (Y < 0) & ~PassiveSet;
+ InfeaSet = NotOptMask & (X < 0) & PassiveSet;
+ NotGood = sum(NonOptSet)+sum(InfeaSet);
+ NotOptCols = NotGood > 0;
+ end
+endfunction
+%------------------------------------------------------------------------------------------------------------------------
+function [ Z,iter ] = solveNormalEqComb( AtA,AtB,PassSet )
+% Solve normal equations using combinatorial grouping.
+% Although this function was originally adopted from the code of
+% "M. H. Van Benthem and M. R. Keenan, J. Chemometrics 2004; 18: 441-450",
+% important modifications were made to fix bugs.
+%
+% Modified by Jingu Kim (jingu@cc.gatech.edu)
+% School of Computational Science and Engineering,
+% Georgia Institute of Technology
+%
+% Last updated Aug-12-2009
+
+ iter = 0;
+ if (nargin ==2) || isempty(PassSet) || all(PassSet(:))
+ Z = AtA\AtB;
+ iter = iter + 1;
+ else
+ Z = zeros(size(AtB));
+ [n,k1] = size(PassSet);
+
+ ## Fixed on Aug-12-2009
+ if k1==1
+ Z(PassSet)=AtA(PassSet,PassSet)\AtB(PassSet);
+ else
+ ## Fixed on Aug-12-2009
+ % The following bug was identified by investigating a bug report by Hanseung Lee.
+ [sortedPassSet,sortIx] = sortrows(PassSet');
+ breaks = any(diff(sortedPassSet)');
+ breakIx = [0 find(breaks) k1];
+ % codedPassSet = 2.^(n-1:-1:0)*PassSet;
+ % [sortedPassSet,sortIx] = sort(codedPassSet);
+ % breaks = diff(sortedPassSet);
+ % breakIx = [0 find(breaks) k1];
+
+ for k=1:length(breakIx)-1
+ cols = sortIx(breakIx(k)+1:breakIx(k+1));
+ vars = PassSet(:,sortIx(breakIx(k)+1));
+ Z(vars,cols) = AtA(vars,vars)\AtB(vars,cols);
+ iter = iter + 1;
+ end
+ end
+ end
+endfunction
+
+%!shared m, n, k, A
+%! m = 30;
+%! n = 20;
+%! k = 10;
+%! A = rand(m,n);
+
+%!test
+%! [W,H,iter,HIS]=nmf_bpas(A,k);
+
+%!test
+%! [W,H,iter,HIS]=nmf_bpas(A,k,'verbose',2);
+
+%!test
+%! [W,H,iter,HIS]=nmf_bpas(A,k,'verbose',1,'nnlssolver','as');
+
+%!test
+%! [W,H,iter,HIS]=nmf_bpas(A,k,'verbose',1,'type','sparse');
+
+%!test
+%! [W,H,iter,HIS]=nmf_bpas(A,k,'verbose',1,'type','sparse','nnlssolver','bp','alpha',1.1,'beta',1.3);
+
+%!test
+%! [W,H,iter,HIS]=nmf_bpas(A,k,'verbose',2,'type','plain','winit',rand(m,k));
+
+%!demo
+%! m = 300;
+%! n = 200;
+%! k = 10;
+%!
+%! W_org = rand(m,k);, W_org(rand(m,k)>0.5)=0;
+%! H_org = rand(k,n);, H_org(rand(k,n)>0.5)=0;
+%!
+%! % normalize W, since 'nmf' normalizes W before return
+%! norm2=sqrt(sum(W_org.^2,1));
+%! toNormalize = norm2>0;
+%! W_org(:,toNormalize) = W_org(:,toNormalize)./repmat(norm2(toNormalize),m,1);
+%!
+%! A = W_org * H_org;
+%!
+%! [W,H,iter,HIS]=nmf_bpas (A,k,'type','plain','tol',1e-4);
+%!
+%! % -------------------- column reordering before computing difference
+%! reorder = zeros(k,1);
+%! selected = zeros(k,1);
+%! for i=1:k
+%! for j=1:k
+%! if ~selected(j), break, end
+%! end
+%! minIx = j;
+%!
+%! for j=minIx+1:k
+%! if ~selected(j)
+%! d1 = norm(W(:,i)-W_org(:,minIx));
+%! d2 = norm(W(:,i)-W_org(:,j));
+%! if (d2<d1)
+%! minIx = j;
+%! end
+%! end
+%! end
+%! reorder(i) = minIx;
+%! selected(minIx) = 1;
+%! end
+%!
+%! W_org = W_org(:,reorder);
+%! H_org = H_org(reorder,:);
+%! % ---------------------------------------------------------------------
+%!
+%! recovery_error_W = norm(W_org-W)/norm(W_org)
+%! recovery_error_H = norm(H_org-H)/norm(H_org)
git://git.creatis.insa-lyon.fr / CreaPhase.git / blobdiff