octave_packages/linear-algebra-2.2.0/nmf_pg.m

   1 ## Copyright (C) 2005-2006 Chih-Jen Lin <cjlin@csie.ntu.edu.tw>
   2 ## All rights reserved.
   3 ##
   4 ## Redistribution and use in source and binary forms, with or without
   5 ## modification, are permitted provided that the following conditions are met:
   6 ##
   7 ##     1 Redistributions of source code must retain the above copyright notice,
   8 ##       this list of conditions and the following disclaimer.
   9 ##     2 Redistributions in binary form must reproduce the above copyright
  10 ##       notice, this list of conditions and the following disclaimer in the
  11 ##       documentation and/or other materials provided with the distribution.
  12 ##
  13 ## THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ''AS IS''
  14 ## AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
  15 ## IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
  16 ## ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHOR OR CONTRIBUTORS BE LIABLE FOR
  17 ## ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  18 ## DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  19 ## SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  20 ## CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  21 ## OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  22 ## OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  23 ##
  24 ## The views and conclusions contained in the software and documentation are
  25 ## those of the authors and should not be interpreted as representing official
  26 ## policies, either expressed or implied, of the copyright holders.
  27
  28 ## -*- texinfo -*-
  29 ## @deftypefn {Function File} {[@var{W}, @var{H}] =} nmf_pg (@var{V}, @var{Winit}, @
  30 ## @var{Hinit}, @var{tol}, @var{timelimit}, @var{maxiter})
  31 ##
  32 ## Non-negative matrix factorization by alternative non-negative least squares
  33 ## using projected gradients.
  34 ##
  35 ## The matrix @var{V} is factorized into two possitive matrices @var{W} and
  36 ## @var{H} such that @code{V = W*H + U}. Where @var{U} is a matrix of residuals
  37 ## that can be negative or positive. When the matrix @var{V} is positive the order
  38 ## of the elements in @var{U} is bounded by the optional named argument @var{tol}
  39 ## (default value @code{1e-9}).
  40 ##
  41 ## The factorization is not unique and depends on the inital guess for the matrices
  42 ## @var{W} and @var{H}. You can pass this initalizations using the optional
  43 ## named arguments @var{Winit} and @var{Hinit}.
  44 ##
  45 ## timelimit, maxiter: limit of time and iterations
  46 ##
  47 ## Examples:
  48 ##
  49 ## @example
  50 ##   A     = rand(10,5);
  51 ##   [W H] = nmf_pg(A,tol=1e-3);
  52 ##   U     = W*H -A;
  53 ##   disp(max(abs(U)));
  54 ## @end example
  55 ##
  56 ## @end deftypefn
  57
  58 ## 2012 - Modified and adapted to Octave 3.6.1 by
  59 ## Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
  60
  61 function [W, H] = nmf_pg (V, varargin)
  62
  63 # JuanPi Fri 16 Mar 2012 10:49:11 AM CET
  64 # TODO:
  65 # - finish docstring
  66 # - avoid multiple transpositions
  67
  68   # --- Parse arguments --- #
  69   parser = inputParser ();
  70   parser.FunctionName = "nmf_pg";
  71   parser = addParamValue (parser,'Winit', [], @ismatrix);
  72   parser = addParamValue (parser,'Hinit', [], @ismatrix);
  73   parser = addParamValue (parser,'Tol', 1e-6, @(x)x>0);
  74   parser = addParamValue (parser,'TimeLimit', 10, @(x)x>0);
  75   parser = addParamValue (parser,'MaxIter', 100, @(x)x>0);
  76   parser = addParamValue (parser,'MaxSubIter', 1e3, @(x)x>0);
  77   parser = addParamValue (parser,'Verbose', true);
  78   parser = parse(parser,varargin{:});
  79
  80   Winit      = parser.Results.Winit;
  81   Hinit      = parser.Results.Hinit;
  82   tol        = parser.Results.Tol;
  83   timelimit  = parser.Results.TimeLimit;
  84   maxiter    = parser.Results.MaxIter;
  85   maxsubiter = parser.Results.MaxSubIter;
  86   verbose    = parser.Results.Verbose;
  87
  88   clear parser
  89   # ------ #
  90
  91   # --- Initialize matrices --- #
  92   [r c] = size (V);
  93   Hgiven = !isempty (Hinit);
  94   Wgiven = !isempty (Winit);
  95   if Wgiven && !Hgiven
  96
  97     W = Winit;
  98     H = ones (size (W,2),c);
  99
 100   elseif !Wgiven && Hgiven
 101
 102     H = Hinit;
 103     W = ones (r, size(H,2));
 104
 105   elseif !Wgiven && !Hgiven
 106
 107     if r == c
 108       W = ones (r)
 109       H = W
 110     else
 111       W = ones (r);
 112       H = ones (r,c);
 113     end
 114
 115   else
 116
 117     W = Winit;
 118     H = Hinit;
 119
 120   end
 121   [Hr,Hc] = size(H);
 122   [Wr,Wc] = size(W);
 123
 124   # start tracking time
 125   initt = cputime ();
 126
 127   gradW = W*(H*H') - V*H';
 128   gradH = (W'*W)*H - W'*V;
 129
 130   initgrad = norm([gradW; gradH'],'fro');
 131
 132   # Tolerances for matrices
 133   tolW = max(0.001,tol)*initgrad;
 134   tolH = tolW;
 135
 136   # ------ #
 137
 138   # --- Main Loop --- #
 139   if verbose
 140     fprintf ('--- Factorizing %d-by-%d matrix into %d-by-%d times %d-by-%d\n',...
 141          r,c,Wr,Wc,Hr,Hc);
 142     fprintf ("Initial gradient norm = %f\n", initgrad);
 143     fflush (stdout);
 144     text_waitbar(0,'Please wait ...');
 145   end
 146
 147   for iter = 1:maxiter
 148
 149     # stopping condition
 150     projnorm = norm ( [ gradW(gradW<0 | W>0); gradH(gradH<0 | H>0) ] );
 151     stop_cond = [projnorm < tol*initgrad , cputime-initt > timelimit];
 152     if any (stop_cond)
 153       if stop_cond(2)
 154         warning('mnf_pg:MaxIter',["Time limit exceeded.\n" ...
 155                                   "Could be solved increasing TimeLimit.\n"]);
 156       end
 157       break
 158     end
 159
 160
 161     # FIXME: avoid multiple transpositions
 162     [W, gradW, iterW] = nlssubprob(V', H', W', tolW, maxsubiter, verbose);
 163     W = W';
 164     gradW = gradW';
 165
 166     if iterW == 1,
 167       tolW = 0.1 * tolW;
 168     end
 169
 170     [H, gradH, iterH] = nlssubprob(V, W, H, tolH, maxsubiter, verbose);
 171     if iterH == 1,
 172       tolH = 0.1 * tolH;
 173     end
 174
 175     if (iterW == 1 && iterH == 1 && tolH + tolW < tol*initgrad),
 176       warning ('nmf_pg:InvalidArgument','Failed to move');
 177       break
 178     end
 179
 180      if verbose
 181       text_waitbar (iter/maxiter);
 182      end
 183   end
 184
 185   if iter == maxiter
 186     warning('mnf_pg:MaxIter',["Reached maximum iterations in main loop.\n" ...
 187                               "Could be solved increasing MaxIter.\n"]);
 188   end
 189
 190   if verbose
 191     fprintf ('\nIterations = %d\nFinal proj-grad norm = %f\n', iter, projnorm);
 192     fflush (stdout);
 193   end
 194 endfunction
 195
 196 function [H, grad,iter] = nlssubprob(V,W,Hinit,tol,maxiter,verbose)
 197 % H, grad: output solution and gradient
 198 % iter: #iterations used
 199 % V, W: constant matrices
 200 % Hinit: initial solution
 201 % tol: stopping tolerance
 202 % maxiter: limit of iterations
 203
 204   H = Hinit;
 205   WtV = W'*V;
 206   WtW = W'*W;
 207
 208   alpha = 1;
 209   beta = 0.1;
 210
 211   for iter=1:maxiter
 212     grad = WtW*H - WtV;
 213     projgrad = norm ( grad(grad < 0 | H >0) );
 214
 215     if projgrad < tol,
 216       break
 217     end
 218
 219     % search step size
 220     Hn = max(H - alpha*grad, 0);
 221     d = Hn-H;
 222     gradd = sum ( sum (grad.*d) );
 223     dQd = sum ( sum ((WtW*d).*d) );
 224
 225     if gradd + 0.5*dQd > 0.01*gradd,
 226       % decrease alpha
 227       while 1,
 228         alpha *= beta;
 229         Hn = max (H - alpha*grad, 0);
 230         d = Hn-H;
 231         gradd = sum (sum (grad.*d) );
 232         dQd = sum (sum ((WtW*d).*d));
 233
 234         if gradd + 0.5*dQd <= 0.01*gradd || alpha < 1e-20
 235           H = Hn;
 236           break
 237         end
 238
 239       endwhile
 240
 241     else
 242       % increase alpha
 243       while 1,
 244         Hp = Hn;
 245         alpha /= beta;
 246         Hn = max (H - alpha*grad, 0);
 247         d = Hn-H;
 248         gradd = sum ( sum (grad.*d) );
 249         dQd = sum (sum ( (WtW*d).*d ) );
 250
 251         if gradd + 0.5*dQd > 0.01*gradd || Hn == Hp || alpha > 1e10
 252           H = Hp;
 253           alpha *= beta;
 254           break
 255         end
 256
 257       endwhile
 258     end
 259
 260   endfor
 261
 262   if iter == maxiter
 263     warning('mnf_pg:MaxIter',["Reached maximum iterations in nlssubprob\n" ...
 264                               "Could be solved increasing MaxSubIter.\n"]);
 265   end
 266
 267 endfunction
 268
 269 %!demo
 270 %! t = linspace (0,1,100)';
 271 %!
 272 %! ## --- Build hump functions of different frequency
 273 %! W_true = arrayfun ( @(f)sin(2*pi*f*t).^2, linspace (0.5,2,4), ...
 274 %!                                                     'uniformoutput', false );
 275 %! W_true = cell2mat (W_true);
 276 %! ## --- Build combinator vectors
 277 %! c      = (1:4)';
 278 %! H_true = arrayfun ( @(f)circshift(c,f), linspace (0,3,4), ...
 279 %!                                                     'uniformoutput', false );
 280 %! H_true = cell2mat (H_true);
 281 %! ## --- Mix them
 282 %! V = W_true*H_true;
 283 %! ## --- Give good inital guesses
 284 %! Winit = W_true + 0.4*randn(size(W_true));
 285 %! Hinit = H_true + 0.2*randn(size(H_true));
 286 %! ## --- Factorize
 287 %! [W H] = nmf_pg(V,'Winit',Winit,'Hinit',Hinit,'Tol',1e-6,'MaxIter',1e3);
 288 %! disp('True mixer')
 289 %! disp(H_true)
 290 %! disp('Rounded factorized mixer')
 291 %! disp(round(H))
 292 %! ## --- Plot results
 293 %! plot(t,W,'o;factorized;')
 294 %! hold on
 295 %! plot(t,W_true,'-;True;')
 296 %! hold off
 297 %! axis tight