view 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]