view toolboxes/distance_learning/mlr/util/qplcprog.m @ 0:e9a9cd732c1e tip

first hg version after svn

function [x, fval, status] = qplcprog(H, f, Ai, bi, Ae, be, lb, ub, varargin)
%QPLCPROG Positive-semidefinite quadratic programming based on linear complementary programming.
%
%	x = qplcprog(H, f);
%	[x, fval, status] = qplcprog(H, f, Ai, bi, Ae, be, lb, ub, varargin);
%
% Inputs:
%	H - Hessian for objective function (n x n matrix)
%
% Optional Inputs:
%	f - Vector for objective function, default is 0 (n x 1 vector)
%	Ai - Linear inequality constraint matrix, default is [] (ni x n matrix)
%	bi - Linear inequality constraint vector, default is [] (ni x 1 vector)
%	Ae - Linear equality constraint matrix, default is [] (ne x n matrix)
%	be - Linear equality constraint vector, default is [] (ne x 1 vector)
%	lb - Lower-bound constraint, default is 0 (n x 1 vector)
%	ub - Upper-bound constraint, default is [] (n x 1 vector)
%
%	maxiter - LCP maximum iterations, default is 100000
%	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 - LCP pivot selection tolerance, default is 1.0e-9
%	tolcon - QP constraint tolerance, default is 1.0e-6
%
% Outputs:
%	x - Solution vector (n x 1 vector)
%	fval - Value of objective function at solution (scalar)
%	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
%
% Notes:
%	Depending upon the inputs, various types of positive-semidefinite quadratic programs can be
%	solved. The general problem is
%
%		min  0.5 * x' * H * x + f' * x
%		 x
%
%	subject to any or all of the following constraints
%
%		Ai * x <= bi
%		Ae * x = be
%		lb <= x <= ub
%
%	with the current requirement that lb must be specified and finite.
%
%	It is assumed that H is positive-semidefinite. Note that if H = 0, qplcp solves a linear
%	program.

% Copyright 2008 The MathWorks, Inc.
% $Revision: 1.1.6.3 $   $ Date: $

if nargin < 1
   error('Finance:qplcprog:MissingInputArgument', ...
      'Missing required Hessian matrix for QP.');
end
H = full(H);
n = size(H,1);

if nargin < 2 || isempty(f)
   f = zeros(n,1);
else
   f = full(f);
   f = f(:);
end

if nargin == 3
   error('Finance:qplcprog:IncompleteConstraintSpecification', ...
      'Incomplete specification for linear inequality constraints.');
end
if nargin < 3
   Ai = [];
   bi = [];
else
   Ai = full(Ai);
   bi = full(bi);
   bi = bi(:);
end

if nargin == 5
   error('Finance:qplcprog:IncompleteConstraintSpecification', ...
      'Incomplete specification for linear equality constraints.');
end
if nargin < 5
   Ae = [];
   be = [];
else
   Ae = full(Ae);
   be = full(be);
   be = be(:);
end

if nargin < 7 || isempty(lb)
   error('Finance:qplcprog:MissingBoundConstraint', ...
      'Lower-bound constraint required.');
else
   lb = full(lb);
   lb = lb(:);
   if any(~isfinite(lb))
      error('Finance:qplcprog:MissingBoundConstraint', ...
         'Finite lower-bound constraint required.');
   end
end

if nargin < 8
   ub = [];
else
   ub = full(ub);
   ub = ub(:);
   if all(ub == Inf)
      ub = [];
   end
end
if any(~isfinite(ub))
   error('Finance:qplcprog:InvalidBoundConstraint', ...
      'If specified, finite upper-bound constraint required.');
end

% process name-value pairs (if any)

if nargin > 8
   if mod(nargin-8,2) ~= 0
      error('Finance:qplcprog:InvalidNameValuePair', ...
         'Invalid name-value pairs with either a missing name or value.');
   end

   names = { 'maxiter', 'tiebreak', 'tolcon', 'tolpiv' };			% names
   values = { 10000, 'first', 1.0e-4, 1.0e-9 };					% default values
   try
      [maxiter, tiebreak, tolcon, tolpiv] = parsepvpairs(names, values, varargin{:});

   catch E
      E.throw
   end
else
   maxiter = 100000;
   tiebreak = 'first';
   tolpiv = 1.0e-9;
   tolcon = 1.0e-6;
end

% check arguments

if maxiter <= 0
   error('Finance:qplcprog:NonPositiveInteger', ...
      'Maximum number of iterations (maxiter) must be a positive integer.');
end

if tolcon < 2*eps
   error('Finance:qplcprog:InvalidTolerance', ...
      'Unrealistically small constraint tolerance (tolcon) specified.');
end

if tolpiv < 2*eps
   error('Finance:qplcprog:InvalidTolerance', ...
      'Unrealistically small pivot tolerance (tolpiv) specified.');
end

% set translation if necessary

if norm(lb) ~= 0
   fx = f + H*lb;
   if ~isempty(Ai) && ~isempty(bi)
      bxi = bi - Ai*lb;
   end
   if ~isempty(Ae) && ~isempty(be)
      bxe = be - Ae*lb;
   end
   if ~isempty(ub)
      uxb = ub - lb;
   end
   translate = true;
else
   fx = f;
   if ~isempty(Ai) && ~isempty(bi)
      bxi = bi;
   end
   if ~isempty(Ae) && ~isempty(be)
      bxe = be;
   end
   if ~isempty(ub)
      uxb = ub;
   end
   translate = false;
end

% set up constraints

if ~isempty(Ai) && ~isempty(bi)
   if ~isempty(Ae) && ~isempty(be)
      if ~isempty(ub)
         A = [ Ai; Ae; -Ae; eye(n,n) ];
         b = [ bxi; bxe; -bxe; uxb ];
      else
         A = [ Ai; Ae; -Ae ];
         b = [ bxi; bxe; -bxe ];
      end
   else
      if ~isempty(ub)
         A = [ Ai; eye(n,n) ];
         b = [ bxi; uxb ];
      else
         A = Ai;
         b = bxi;
      end
   end
else
   if ~isempty(Ae) && ~isempty(be)
      if ~isempty(ub)
         A = [ Ae; -Ae; eye(n,n) ];
         b = [ bxe; -bxe; uxb ];
      else
         A = [ Ae; -Ae ];
         b = [ bxe; -bxe ];
      end
   else
      if ~isempty(ub)
         A = eye(n,n);
         b = uxb;
      else
         error('Finance:qplcprog:InvalidConstraints', ...
            'Insufficient or invalid constraints.');
      end
   end
end

% set up lcp problem

ns = size(A,1);

M = [ zeros(ns,ns) -A; A' H ];
q = [ b; fx ];

% solve lcp problem

if isempty(varargin)
   [z, w, status] = lcprog(M, q); %#ok
else
   [z, w, status] = lcprog(M, q, 'maxiter', maxiter, 'tiebreak', tiebreak, 'tolpiv', tolpiv); %#ok
end

% solve qp problem

if status > 0
   x = z(end-n+1:end);
   if translate
      x = x + lb;
   end
   fval = 0.5*x'*H*x + f'*x;

   if qplcprogrealitycheck(Ai, bi, Ae, be, lb, ub, x)
      status = -2;
      x = NaN(n,1);
      fval = NaN;
   end
else
   x = NaN(n,1);
   fval = NaN;
end

function notok = qplcprogrealitycheck(Ai, bi, Ae, be, lb, ub, x, tolcon)
%QPLCPROGREALITYCHECK - Make sure candidate solution x is feasible.

if nargin < 7
   error('Finance:qplcprog:MissingInputArgument', ...
      'Missing required input arguments Ai, bi, Ae, be, lb, ub, x.');
end
if nargin < 8 || isempty(tolcon)
   tolcon = 1.0e-6;
end

notok = false;

if ~isempty(Ai)
   u = max(Ai*x - bi,0);
   if u > tolcon
      warning('Finance:qplcprog:InfeasibleProblem', ...
         'Candidate solution violates inequality constraints.');
      notok = true;
   end
end

if ~isempty(Ae)
   if norm(be - Ae*x,inf) > tolcon
      warning('Finance:qplcprog:InfeasibleProblem', ...
         'Candidate solution violates equality constraints.');
      notok = true;
   end
end

if ~isempty(lb)
   u = x - lb;
   if min(u) < -tolcon
      warning('Finance:qplcprog:InfeasibleProblem', ...
         'Candidate solution violates lower-bound constraints.');
      notok = true;
   end
end

if ~isempty(ub)
   u = ub - x;
   if min(u) < -tolcon
      warning('Finance:qplcprog:InfeasibleProblem', ...
         'Candidate solution violates upper-bound constraints.');
      notok = true;
   end
end


% [EOF]