comparison toolboxes/distance_learning/mlr/util/lcprog.m @ 0:e9a9cd732c1e tip

first hg version after svn

+function [z, w, status, z0] = lcprog(M, q, varargin)
+%LCPROG - Linear complementary programming with Lemke's algorithm.
+%
+%	[z, w] = lcprog(M, q);
+%	[z, w, status, z0] = lcprog(M, q, varargin);
+%
+% Inputs:
+%	M - Matrix (n x n matrix)
+%	q - Vector (n vector)
+%
+% Optional Inputs:
+%	Optional name-value pairs are:
+%	'maxiter'	Maximum number of iterations before termination of algorithm (default is (1+n)^3).
+%	'tiebreak'	Method to break ties for pivot selection (default is 'first'). Options are:
+%				'first'		Select pivot with lowest index.
+%				'last'		Select pivot with highest index.
+%				'random'	Select a pivot at random.
+%	'tolpiv'	Pivot tolerance below which a number is considered to be zero (default is 1.0e-9).
+%
+% Outputs:
+%	z - Solution vector (n vector)
+%	w - Complementary solution vector (n vector)
+%	status - Status of optimization with values (integer)
+%		 6		Ray termination
+%		 1		Convergence to a solution (normal termination)
+%		 0		Termination prior to convergence
+%		-2		Infeasible problem
+%		-3		Bad pivot to an infeasible solution
+%	z0 - Artificial variable which is 0 upon normal termination (scalar)
+%
+% Comments:
+%	Given a matrix M and a vector q, the linear complementary programming (LCP) problem seeks an
+%	orthogonal pair of non-negative vectors w and z such that
+%		w = M * z + q
+%	with
+%		w' * z = 0
+%		w >= 0
+%		z >= 0 .
+%
+%	The algorithm starts with the problem
+%		w = M * z + q + z0 * 1
+%	with an artificial variable z0 > 0 such that q + z0*1 >= 0. The algorithm goes through a series
+%	of pivot steps that seeks to drive z0 to 0. The algorithm can stop in two ways with either
+%	normal termination or ray termination.
+%
+%	With normal termination, the artificial variable z0 is zero so that the original problem is
+%	solved.
+%
+%	With ray termination, the artificial variable z0 is non-negative with "ray" solution in the form
+%		w = wR + lambda * dw
+%		z = zR + lambda * dz
+%		z0 = z0R + lambda * dz0
+%	that satisfies
+%		w = M * z + q + z0 * 1
+%		w' * z = 0
+%	for any lambda >= 0. If ray termination occurs, the returned variables w, z, and z0 are the
+%	points on the ray with lambda = 0.
+%
+%	If tolpiv is too small, the algorithm may identify a pivot that is effectively zero. If this
+%	situation occurs, the status code -3 is triggered. Try a larger pivot tolerance 'tolpiv' although
+%	a pivot that is too large can trigger false convergence due to failure to find a pivot.
+%
+
+%
+% Copyright 2007 The MathWorks, Inc.
+% $Revision: 1.1.6.4 $   $ Date: $
+
+% initialization
+
+if nargin < 2 || isempty(M) || isempty(q)
+   error('finance:lcprog:MissingInputArgument', ...
+      'One or more missing required input arguments M, q.');
+end
+
+q = q(:);
+n = size(q,1);
+
+% process name-value pairs (if any)
+
+if mod(nargin-2,2) ~= 0
+   error('finance:lcprog:InvalidNameValuePair', ...
+      'Invalid name-value pairs with either a missing name or value.');
+end
+
+names = { 'maxiter', 'tiebreak', 'tolpiv' };				% names
+values = { (1 + n)^3, 'first', 1.0e-9 };					% default values
+try
+   [maxiter, tiebreak, tolpiv] = parsepvpairs(names, values, varargin{:});
+
+catch E
+   E.throw
+end
+
+% check optional arguments
+
+if maxiter <= 0
+   error('finance:lcprog:NonPositiveInteger', ...
+      'Maximum number of iterations (maxiter) must be a positive integer.');
+end
+
+if ~isempty(strmatch(tiebreak,{'default','first'}))
+   tb = 1;
+elseif ~isempty(strmatch(tiebreak,'last'))
+   tb = 2;
+elseif ~isempty(strmatch(tiebreak,'random'))
+   tb = 3;
+else
+   tb = 1;
+end
+
+if tolpiv < 2*eps
+   error('finance:lcprog:InvalidTolerance', ...
+      'Unrealistically small tolerance (tolpiv) specified.');
+end
+
+% trivial solution
+
+if all(q >= 0)
+   z = zeros(n,1);
+   w = q;
+   status = 1;
+   z0 = 0;
+   return
+end
+
+% set up general problem
+
+w = zeros(n,1);
+z = zeros(n,1);
+u = ones(n,1);
+z0 = max(-q);
+
+tableau = [ eye(n), -M, -u, q];
+basic = (1:n)';
+
+% initial pivot and initial basis
+
+r = find(q == min(q),1);			% r are pivot rows
+cx = r;								% w(r) leaves the basis
+c = 2*n + 1;						% initial pivot column
+basic(r) = c;						% z0 enters the basis
+
+% pivot step
+
+pr = (1/tableau(r,c)) * tableau(r,:);
+pc = tableau(:,c);
+tableau = tableau - pc * pr;
+tableau(r,:) = pr;
+
+% main loop
+
+for iter = 2:maxiter
+
+   % get new basic variable
+
+   if cx <= n
+      c = n + cx;
+   elseif cx <= 2*n
+      c = cx - n;
+   end
+
+   % do ratio test
+
+   d = tableau(:,c);
+   mind = NaN(n,1);
+   for i = 1:n
+      if d(i) > tolpiv
+         mind(i) = tableau(i,end)/d(i);
+      end
+   end
+   iind = find(mind == nanmin(mind));
+
+   % ray termination
+
+   if isempty(iind)						% cannot find a new basic variable
+      for i = 1:n
+         if basic(i) <= n
+            ii = basic(i);
+            w(ii) = tableau(i,end);
+            z(ii) = 0;
+         elseif basic(i) <= 2*n
+            ii = basic(i) - n;
+            w(ii) = 0;
+            z(ii) = tableau(i,end);
+         else
+            z0 = tableau(i,end);
+         end
+      end
+      if lcprealitycheck(M, q, w, z)
+         status = -3;
+         z = NaN(n,1);
+         w = NaN(n,1);
+         z0 = NaN;
+      else
+         status = 6;
+      end
+      return
+   end
+
+   % identify next basic variable to be replaced using ratio test
+
+   if tb == 3
+      ii = ceil(rand()*numel(iind));
+      if ii < 1
+         ii = 1;
+      end
+      r = iind(ii);
+   elseif tb == 2
+      r = iind(end);
+   elseif tb == 1
+      r = iind(1);
+   end
+   for i = 1:numel(iind)					% always select z0 if it is in the index set
+      if basic(iind(i)) == (2*n + 1)
+         r = iind(i);
+      end
+   end
+
+   % new basis
+
+   cx = basic(r);							% get variable leaving the basis
+   basic(r) = c;							% move new basic variable into the basis
+
+   % pivot step
+
+   pr = (1/tableau(r,c)) * tableau(r,:);
+   pc = tableau(:,c);
+   tableau = tableau - pc * pr;
+   tableau(r,:) = pr;
+
+   if ~isfinite(tableau(1,end))			% trap overflow/underflow
+      warning('finance:lcprog:InfeasibleProblem', ...
+         'Tableau overflow/underflow due to bad pivot.');
+      status = -3;
+      z = NaN(n,1);
+      w = NaN(n,1);
+      z0 = NaN;
+      return
+   end
+
+   % convergence test
+
+   if all(basic <= 2*n)					% test if z0 left the basis -> convergence
+      for i = 1:n
+         if basic(i) <= n
+            ii = basic(i);
+            w(ii) = tableau(i,end);
+            z(ii) = 0;
+         else
+            ii = basic(i) - n;
+            w(ii) = 0;
+            z(ii) = tableau(i,end);
+         end
+      end
+      z0 = 0;
+      if lcprealitycheck(M, q, w, z)
+         status = -3;
+         z = NaN(n,1);
+         w = NaN(n,1);
+         z0 = NaN;
+      else
+         status = 1;
+      end
+      return
+   end
+end
+
+% too many iterations
+
+for i = 1:n
+   if basic(i) <= n
+      ii = basic(i);
+      w(ii) = tableau(i,end);
+      z(ii) = 0;
+   elseif basic(i) <= 2*n
+      ii = basic(i) - n;
+      w(ii) = 0;
+      z(ii) = tableau(i,end);
+   else
+      z0 = tableau(i,end);
+   end
+end
+if lcprealitycheck(M, q, w, z)
+   status = -3;
+   z = NaN(n,1);
+   w = NaN(n,1);
+   z0 = NaN;
+else
+   status = 0;
+end
+
+
+function notok = lcprealitycheck(M, q, w, z)
+%LCPREALITYCHECK - Make sure candidate solution z and w is complementary basic feasible.
+
+u = w - M*z - q;
+
+if (max(u) - min(u)) > 1
+   warning('finance:lcprog:InfeasibleProblem', ...
+      'Candidate solution is infeasible due to a bad pivot.');
+   notok = true;
+else
+   notok = false;
+end
+
+
+% [EOF]