Mercurial > hg > octave-lojdl
changeset 8406:0b7566aea627
fix & optimize lsqnonneg
author | Jaroslav Hajek <highegg@gmail.com> |
---|---|
date | Mon, 15 Dec 2008 12:59:34 +0100 |
parents | d79dfbff2f9d |
children | 096c22ce2b0b |
files | scripts/ChangeLog scripts/optimization/lsqnonneg.m |
diffstat | 2 files changed, 30 insertions(+), 17 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ChangeLog +++ b/scripts/ChangeLog @@ -1,3 +1,9 @@ +2008-12-15 Jaroslav Hajek <highegg@gmail.com> + + * optimization/lsqnonneg.m: Preprocess using QR for over-determined + systems. Simplify & fix indexing. Use left division for step problem. + Fix output args. + 2008-12-13 Francesco Potort́ <pot@gnu.org> * specfun/nchoosek.m: Check for input arguments, signal loss of
--- a/scripts/optimization/lsqnonneg.m +++ b/scripts/optimization/lsqnonneg.m @@ -1,4 +1,5 @@ ## Copyright (C) 2008 Bill Denney +## Copyright (C) 2008 Jaroslav Hajek ## ## This file is part of Octave. ## @@ -67,17 +68,28 @@ x = zeros (columns (c), 1); endif + MaxIter = optimget (options, "MaxIter", 1e5); ## Initialize the values. - p = []; - z = 1:numel (x); + p = false (1, numel (x)); + z = ~p; + ## If the problem is significantly over-determined, preprocess it using a + ## QR factorization first. + if (rows (c) >= 1.5 * columns (c)) + [q, r] = qr (c, 0); + d = q'*d; + c = r; + clear q + endif ## LH2: compute the gradient. w = c'*(d - c*x); + xtmp = zeros (columns (c), 1); + iter = 0; ## LH3: test for completion. - while (! isempty (z) && any (w(z) > 0) && iter < MaxIter) + while (any (z) && any (w(z) > 0) && iter < MaxIter) ## LH4: find the maximum gradient. idx = find (w == max (w)); if (numel (idx) > 1) @@ -86,8 +98,8 @@ idx = idx(1); endif ## LH5: move the index from Z to P. - z(z == idx) = []; - p(end+1) = idx; + z(idx) = false; + p(idx) = true; newx = false; while (! newx && iter < MaxIter) @@ -95,12 +107,9 @@ ## LH6: compute the positive matrix and find the min norm solution ## of the positive problem. - cpos = c; - cpos(:,z) = 0; ## Find min norm solution to the positive matrix. - [cpos_q, cpos_r] = qr (cpos, 0); - xtmp = (cpos_r\cpos_q')*d; - xtmp(z) = 0; + xtmp(:) = 0; + xtmp(p) = c(:,p) \ d; if (all (xtmp >= 0)) ## LH7: tmp solution found, iterate. newx = true; @@ -112,9 +121,8 @@ ## LH10: adjust X. x = x + alpha*(xtmp - x); ## LH11: move from P to Z all X == 0. - idx = find (x == 0); - p(ismember (p, idx)) = []; - z = [z idx]; + z |= (x == 0); + p = ~z; endif endwhile w = c'*(d - c*x); @@ -123,10 +131,10 @@ ## Generate the additional output arguments. if (nargout > 1) - resnorm = norm (C*x-d) ^ 2; + resnorm = norm (c*x - d) ^ 2; endif if (nargout > 2) - residual = d-C*x; + residual = d - c*x; endif exitflag = iter; if (nargout > 3 && iter >= MaxIter) @@ -136,8 +144,7 @@ output = struct ("algorithm", "nnls", "iterations", iter); endif if (nargout > 5) - ## FIXME - error ("lsqnonneg: lambda is not yet implemented"); + lambda = w; endif endfunction