Menu
▾
▴
[Octave-cvsupdate] SF.net SVN: octave:[9918] trunk/octave-forge/main/linear-algebra
|
From: <jpi...@us...> - 2012-03-16 11:24:15
|
Revision: 9918
http://octave.svn.sourceforge.net/octave/?rev=9918&view=rev
Author: jpicarbajal
Date: 2012-03-16 11:24:04 +0000 (Fri, 16 Mar 2012)
Log Message:
-----------
linear-algebra: Beginning to add nonnegative matrix factorization
Added Paths:
-----------
trunk/octave-forge/main/linear-algebra/devel/
trunk/octave-forge/main/linear-algebra/devel/nmf_pg.m
Added: trunk/octave-forge/main/linear-algebra/devel/nmf_pg.m
===================================================================
--- trunk/octave-forge/main/linear-algebra/devel/nmf_pg.m (rev 0)
+++ trunk/octave-forge/main/linear-algebra/devel/nmf_pg.m 2012-03-16 11:24:04 UTC (rev 9918)
@@ -0,0 +1,249 @@
+## Copyright (C) 2005-2006 Chih-Jen Lin <cj...@cs...>
+## All rights reserved.
+##
+## Redistribution and use in source and binary forms, with or without
+## modification, are permitted provided that the following conditions are met:
+##
+## 1 Redistributions of source code must retain the above copyright notice,
+## this list of conditions and the following disclaimer.
+## 2 Redistributions in binary form must reproduce the above copyright
+## notice, this list of conditions and the following disclaimer in the
+## documentation and/or other materials provided with the distribution.
+##
+## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
+## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
+## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
+## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
+## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
+## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
+## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+##
+## The views and conclusions contained in the software and documentation are
+## those of the authors and should not be interpreted as representing official
+## policies, either expressed or implied, of the copyright holders.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{W}, @var{H}] =} nmf_pg (@var{V}, @var{Winit}, @
+## @var{Hinit}, @var{tol}, @var{timelimit}, @var{maxiter})
+##
+## Non-negative matrix factorization by alternative non-negative least squares
+## using projected gradients.
+##
+## The matrix @var{V} is factorized into two possitive matrices @var{W} and
+## @var{H} such that @code{V = W*H + U}. Where @var{U} is a matrix of residuals
+## that can be negative or positive. When the matrix @var{V} is positive the order
+## of the elements in @var{U} is bounded by the optional named argument @var{tol}
+## (default value @code{1e-9}).
+##
+## The factorization is not unique and depends on the inital guess for the matrices
+## @var{W} and @var{H}. You can pass this initalizations using the optional
+## named arguments @var{Winit} and @var{Hinit}.
+##
+## timelimit, maxiter: limit of time and iterations
+##
+## Examples:
+##
+## @example
+## A = rand(10,5);
+## [W H] = nmf_pg(A,tol=1e-3);
+## U = W*H -A;
+## disp(max(abs(U)));
+## @end example
+##
+## @end deftypefn
+
+## 2012 - Modified and adapted to Octave 3.6.1 by
+## Juan Pablo Carbajal <car...@if...>
+
+function [W, H] = nmf_pg (V, Winit=[], Hinit=[], tol=1e-6, timelimit=1, maxiter=100)
+# JuanPi Fri 16 Mar 2012 10:49:11 AM CET
+# TODO:
+# - finish docstring
+# - verbose optional
+# - avoid multiple transpositions
+# - vectorize loops
+ [r c] = size (V);
+ Hgiven = !isempty (Hinit)
+ Wgiven = !isempty (Winit)
+ if Wgiven && !Hgiven
+
+ W = Winit;
+ H = ones (size (W,2),c);
+
+ elseif !Wgiven && Hgiven
+
+ H = Hinit;
+ W = ones (r, size(H,2));
+
+ elseif !Wgiven && !Hgiven
+
+ if r == c
+ W = ones (r)
+ H = W
+ else
+ W = ones (r);
+ H = ones (r,c);
+ end
+
+ else
+
+ W = Winit;
+ H = Hinit;
+
+ end
+
+ initt = cputime ();
+
+ gradW = W*(H*H') - V*H';
+ gradH = (W'*W)*H - W'*V;
+
+ initgrad = norm([gradW; gradH'],'fro');
+
+ fprintf('Init gradient norm %f\n', initgrad);
+
+ tolW = max(0.001,tol)*initgrad;
+ tolH = tolW;
+
+ for iter=1:maxiter,
+
+ % stopping condition
+ projnorm = norm ( [ gradW(gradW<0 | W>0); gradH(gradH<0 | H>0) ] );
+ if projnorm < tol*initgrad || cputime-initt > timelimit,
+ break;
+ end
+
+ [W, gradW, iterW] = nlssubprob(V', H', W', tolW, 1000);
+ W = W';
+ gradW = gradW';
+
+ if iterW == 1,
+ tolW = 0.1 * tolW;
+ end
+
+ [H, gradH, iterH] = nlssubprob(V, W, H, tolH, 1000);
+ if iterH == 1,
+ tolH = 0.1 * tolH;
+ end
+
+ if (iterW == 1 && iterH == 1 && tolH + tolW < tol*initgrad),
+ fprintf('Failed to move\n'); break;
+ end
+
+ if rem(iter,10)==0
+ fprintf('.');
+ end
+
+ end
+
+ fprintf('\nIter = %d Final proj-grad norm %f\n', iter, projnorm);
+
+endfunction
+
+function [H, grad,iter] = nlssubprob(V,W,Hinit,tol,maxiter)
+% H, grad: output solution and gradient
+% iter: #iterations used
+% V, W: constant matrices
+% Hinit: initial solution
+% tol: stopping tolerance
+% maxiter: limit of iterations
+
+ H = Hinit;
+ WtV = W'*V;
+ WtW = W'*W;
+
+ alpha = 1;
+ beta = 0.1;
+
+ for iter=1:maxiter
+ grad = WtW*H - WtV;
+ projgrad = norm ( grad(grad < 0 | H >0) );
+
+ if projgrad < tol,
+ break
+ end
+
+ % search step size
+ Hn = max(H - alpha*grad, 0);
+ d = Hn-H;
+ gradd = sum ( sum (grad.*d) );
+ dQd = sum ( sum ((WtW*d).*d) );
+
+ if gradd + 0.5*dQd > 0.01*gradd,
+ % decrease alpha
+ while 1,
+ alpha *= beta;
+ Hn = max (H - alpha*grad, 0);
+ d = Hn-H;
+ gradd = sum (sum (grad.*d) );
+ dQd = sum (sum ((WtW*d).*d));
+
+ if gradd + 0.5*dQd <= 0.01*gradd || alpha < 1e-20
+ H = Hn;
+ break
+ end
+
+ endwhile
+
+ else
+ % increase alpha
+ while 1,
+ Hp = Hn;
+ alpha /= beta;
+ Hn = max(H - alpha*grad, 0);
+ d = Hn-H;
+ gradd = sum ( sum (grad.*d) );
+ dQd = sum (sum ( (WtW*d).*d ) );
+
+ if gradd + 0.5*dQd > 0.01*gradd || Hn == Hp || alpha > 1e10
+ H = Hp;
+ alpha *= beta;
+ break
+ end
+
+ endwhile
+ end
+
+ endfor
+
+ if iter == maxiter
+ fprintf('Max iter in nlssubprob\n');
+ end
+
+endfunction
+
+%!demo
+%! t = linspace (0,1,100)';
+%!
+%! % Build hump functions of different frequency
+%! W_true = arrayfun ( @(f)sin(2*pi*f*t).^2, linspace (0.5,2,4), ...
+%! 'uniformoutput', false );
+%! W_true = cell2mat (W_true);
+%!
+%! % Build combinator vectors
+%! c = (1:4)';
+%! H_true = arrayfun ( @(f)circshift(c,f), linspace (0,3,4), ...
+%! 'uniformoutput', false );
+%! H_true = cell2mat (H_true);
+%!
+%! % Mix them
+%! V = W_true*H_true;
+%!
+%! % give good inital guesses
+%! Winit = W_true + 0.4*randn(size(W_true));
+%! Hinit = H_true + 0.2*randn(size(H_true));
+%!
+%! Factorize
+%! [W H] = nmf_pg(V,Winit,Hinit,1e-6,1,1e3);
+%! disp('True mixer')
+%! disp(H_true)
+%! disp('Rounded factorized mixer')
+%! disp(round(H))
+%!
+%! plot(t,W,'o;factorized;')
+%! hold on
+%! plot(t,W_true,'-;True;')
+%! hold off
+%! axis tight
This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.
|
