1 %% Copyright (C) 2010, 2011 Olaf Till <i7tiol@t-online.de>
3 %% This program is free software; you can redistribute it and/or modify it under
4 %% the terms of the GNU General Public License as published by the Free Software
5 %% Foundation; either version 3 of the License, or (at your option) any later
8 %% This program is distributed in the hope that it will be useful, but WITHOUT
9 %% ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
10 %% FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
13 %% You should have received a copy of the GNU General Public License along with
14 %% this program; if not, see <http://www.gnu.org/licenses/>.
16 %% [lb, idx, ridx, mv] = cpiv_bard (v, m[, incl])
18 %% v: column vector; m: matrix; incl (optional): index. length (v)
19 %% must equal rows (m). Finds column vectors w and l with w == v + m *
20 %% l, w >= 0, l >= 0, l.' * w == 0. Chooses idx, w, and l so that
21 %% l(~idx) == 0, l(idx) == -inv (m(idx, idx)) * v(idx), w(idx) roughly
22 %% == 0, and w(~idx) == v(~idx) + m(idx, ~idx).' * l(idx). idx indexes
23 %% at least everything indexed by incl, but l(incl) may be < 0. lb:
24 %% l(idx) (column vector); idx: logical index, defined above; ridx:
25 %% ~idx & w roughly == 0; mv: [m, v] after performing a Gauss-Jordan
26 %% 'sweep' (with gjp.m) on each diagonal element indexed by idx.
27 %% Except the handling of incl (which enables handling of equality
28 %% constraints in the calling code), this is called solving the
29 %% 'complementary pivot problem' (Cottle, R. W. and Dantzig, G. B.,
30 %% 'Complementary pivot theory of mathematical programming', Linear
31 %% Algebra and Appl. 1, 102--125. References for the current
32 %% algorithm: Bard, Y.: Nonlinear Parameter Estimation, p. 147--149,
33 %% Academic Press, New York and London 1974; Bard, Y., 'An eclectic
34 %% approach to nonlinear programming', Proc. ANU Sem. Optimization,
35 %% Canberra, Austral. Nat. Univ.).
37 function [lb, idx, ridx, m] = cpiv_bard (v, m, incl)
41 error ('first argument is no column vector'); % the most typical mistake
45 elseif (islogical (incl))
51 if (length (incl) == n)
54 m = cat (2, m, m * v);
62 nz = eps; % This is arbitrary; components of w and -l are regarded as
63 % non-negative if >= -nz.
64 nl = 100 * n; % maximum number of loop repeats, after that give up
69 while (~ready && nl > 0)
70 [vm, idm] = min (sgn(nincl) .* m(nincl, end));
82 error ('not successful');
86 ridx = ~idx & abs (m(:, end)) <= nz;