1 ## Copyright (C) 2005-2006 Chih-Jen Lin <cjlin@csie.ntu.edu.tw>
2 ## All rights reserved.
4 ## Redistribution and use in source and binary forms, with or without
5 ## modification, are permitted provided that the following conditions are met:
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.
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.
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.
29 ## @deftypefn {Function File} {[@var{W}, @var{H}] =} nmf_pg (@var{V}, @var{Winit}, @
30 ## @var{Hinit}, @var{tol}, @var{timelimit}, @var{maxiter})
32 ## Non-negative matrix factorization by alternative non-negative least squares
33 ## using projected gradients.
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}).
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}.
45 ## timelimit, maxiter: limit of time and iterations
51 ## [W H] = nmf_pg(A,tol=1e-3);
58 ## 2012 - Modified and adapted to Octave 3.6.1 by
59 ## Juan Pablo Carbajal <carbajal@ifi.uzh.ch>
61 function [W, H] = nmf_pg (V, varargin)
63 # JuanPi Fri 16 Mar 2012 10:49:11 AM CET
66 # - avoid multiple transpositions
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{:});
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;
91 # --- Initialize matrices --- #
93 Hgiven = !isempty (Hinit);
94 Wgiven = !isempty (Winit);
98 H = ones (size (W,2),c);
100 elseif !Wgiven && Hgiven
103 W = ones (r, size(H,2));
105 elseif !Wgiven && !Hgiven
124 # start tracking time
127 gradW = W*(H*H') - V*H';
128 gradH = (W'*W)*H - W'*V;
130 initgrad = norm([gradW; gradH'],'fro');
132 # Tolerances for matrices
133 tolW = max(0.001,tol)*initgrad;
138 # --- Main Loop --- #
140 fprintf ('--- Factorizing %d-by-%d matrix into %d-by-%d times %d-by-%d\n',...
142 fprintf ("Initial gradient norm = %f\n", initgrad);
144 text_waitbar(0,'Please wait ...');
150 projnorm = norm ( [ gradW(gradW<0 | W>0); gradH(gradH<0 | H>0) ] );
151 stop_cond = [projnorm < tol*initgrad , cputime-initt > timelimit];
154 warning('mnf_pg:MaxIter',["Time limit exceeded.\n" ...
155 "Could be solved increasing TimeLimit.\n"]);
161 # FIXME: avoid multiple transpositions
162 [W, gradW, iterW] = nlssubprob(V', H', W', tolW, maxsubiter, verbose);
170 [H, gradH, iterH] = nlssubprob(V, W, H, tolH, maxsubiter, verbose);
175 if (iterW == 1 && iterH == 1 && tolH + tolW < tol*initgrad),
176 warning ('nmf_pg:InvalidArgument','Failed to move');
181 text_waitbar (iter/maxiter);
186 warning('mnf_pg:MaxIter',["Reached maximum iterations in main loop.\n" ...
187 "Could be solved increasing MaxIter.\n"]);
191 fprintf ('\nIterations = %d\nFinal proj-grad norm = %f\n', iter, projnorm);
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
213 projgrad = norm ( grad(grad < 0 | H >0) );
220 Hn = max(H - alpha*grad, 0);
222 gradd = sum ( sum (grad.*d) );
223 dQd = sum ( sum ((WtW*d).*d) );
225 if gradd + 0.5*dQd > 0.01*gradd,
229 Hn = max (H - alpha*grad, 0);
231 gradd = sum (sum (grad.*d) );
232 dQd = sum (sum ((WtW*d).*d));
234 if gradd + 0.5*dQd <= 0.01*gradd || alpha < 1e-20
246 Hn = max (H - alpha*grad, 0);
248 gradd = sum ( sum (grad.*d) );
249 dQd = sum (sum ( (WtW*d).*d ) );
251 if gradd + 0.5*dQd > 0.01*gradd || Hn == Hp || alpha > 1e10
263 warning('mnf_pg:MaxIter',["Reached maximum iterations in nlssubprob\n" ...
264 "Could be solved increasing MaxSubIter.\n"]);
270 %! t = linspace (0,1,100)';
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
278 %! H_true = arrayfun ( @(f)circshift(c,f), linspace (0,3,4), ...
279 %! 'uniformoutput', false );
280 %! H_true = cell2mat (H_true);
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));
287 %! [W H] = nmf_pg(V,'Winit',Winit,'Hinit',Hinit,'Tol',1e-6,'MaxIter',1e3);
288 %! disp('True mixer')
290 %! disp('Rounded factorized mixer')
292 %! ## --- Plot results
293 %! plot(t,W,'o;factorized;')
295 %! plot(t,W_true,'-;True;')