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