Mercurial > hg > camir-aes2014

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [x, fval, status] = qplcprog(H, f, Ai, bi, Ae, be, lb, ub, varargin)
wolffd@0 2 %QPLCPROG Positive-semidefinite quadratic programming based on linear complementary programming.
wolffd@0 3 %
wolffd@0 4 % x = qplcprog(H, f);
wolffd@0 5 % [x, fval, status] = qplcprog(H, f, Ai, bi, Ae, be, lb, ub, varargin);
wolffd@0 6 %
wolffd@0 7 % Inputs:
wolffd@0 8 % H - Hessian for objective function (n x n matrix)
wolffd@0 9 %
wolffd@0 10 % Optional Inputs:
wolffd@0 11 % f - Vector for objective function, default is 0 (n x 1 vector)
wolffd@0 12 % Ai - Linear inequality constraint matrix, default is [] (ni x n matrix)
wolffd@0 13 % bi - Linear inequality constraint vector, default is [] (ni x 1 vector)
wolffd@0 14 % Ae - Linear equality constraint matrix, default is [] (ne x n matrix)
wolffd@0 15 % be - Linear equality constraint vector, default is [] (ne x 1 vector)
wolffd@0 16 % lb - Lower-bound constraint, default is 0 (n x 1 vector)
wolffd@0 17 % ub - Upper-bound constraint, default is [] (n x 1 vector)
wolffd@0 18 %
wolffd@0 19 % maxiter - LCP maximum iterations, default is 100000
wolffd@0 20 % tiebreak - Method to break ties for pivot selection (default is 'first'). Options are:
wolffd@0 21 % 'first' Select pivot with lowest index.
wolffd@0 22 % 'last' Select pivot with highest index.
wolffd@0 23 % 'random' Select a pivot at random.
wolffd@0 24 % tolpiv - LCP pivot selection tolerance, default is 1.0e-9
wolffd@0 25 % tolcon - QP constraint tolerance, default is 1.0e-6
wolffd@0 26 %
wolffd@0 27 % Outputs:
wolffd@0 28 % x - Solution vector (n x 1 vector)
wolffd@0 29 % fval - Value of objective function at solution (scalar)
wolffd@0 30 % status - Status of optimization with values (integer)
wolffd@0 31 % 6 Ray termination
wolffd@0 32 % 1 Convergence to a solution (normal termination)
wolffd@0 33 % 0 Termination prior to convergence
wolffd@0 34 % -2 Infeasible problem
wolffd@0 35 % -3 Bad pivot to an infeasible solution
wolffd@0 36 %
wolffd@0 37 % Notes:
wolffd@0 38 % Depending upon the inputs, various types of positive-semidefinite quadratic programs can be
wolffd@0 39 % solved. The general problem is
wolffd@0 40 %
wolffd@0 41 % min 0.5 * x' * H * x + f' * x
wolffd@0 42 % x
wolffd@0 43 %
wolffd@0 44 % subject to any or all of the following constraints
wolffd@0 45 %
wolffd@0 46 % Ai * x <= bi
wolffd@0 47 % Ae * x = be
wolffd@0 48 % lb <= x <= ub
wolffd@0 49 %
wolffd@0 50 % with the current requirement that lb must be specified and finite.
wolffd@0 51 %
wolffd@0 52 % It is assumed that H is positive-semidefinite. Note that if H = 0, qplcp solves a linear
wolffd@0 53 % program.
wolffd@0 54
wolffd@0 55 % Copyright 2008 The MathWorks, Inc.
wolffd@0 56 % $Revision: 1.1.6.3 $ $ Date: $
wolffd@0 57
wolffd@0 58 if nargin < 1
wolffd@0 59 error('Finance:qplcprog:MissingInputArgument', ...
wolffd@0 60 'Missing required Hessian matrix for QP.');
wolffd@0 61 end
wolffd@0 62 H = full(H);
wolffd@0 63 n = size(H,1);
wolffd@0 64
wolffd@0 65 if nargin < 2 || isempty(f)
wolffd@0 66 f = zeros(n,1);
wolffd@0 67 else
wolffd@0 68 f = full(f);
wolffd@0 69 f = f(:);
wolffd@0 70 end
wolffd@0 71
wolffd@0 72 if nargin == 3
wolffd@0 73 error('Finance:qplcprog:IncompleteConstraintSpecification', ...
wolffd@0 74 'Incomplete specification for linear inequality constraints.');
wolffd@0 75 end
wolffd@0 76 if nargin < 3
wolffd@0 77 Ai = [];
wolffd@0 78 bi = [];
wolffd@0 79 else
wolffd@0 80 Ai = full(Ai);
wolffd@0 81 bi = full(bi);
wolffd@0 82 bi = bi(:);
wolffd@0 83 end
wolffd@0 84
wolffd@0 85 if nargin == 5
wolffd@0 86 error('Finance:qplcprog:IncompleteConstraintSpecification', ...
wolffd@0 87 'Incomplete specification for linear equality constraints.');
wolffd@0 88 end
wolffd@0 89 if nargin < 5
wolffd@0 90 Ae = [];
wolffd@0 91 be = [];
wolffd@0 92 else
wolffd@0 93 Ae = full(Ae);
wolffd@0 94 be = full(be);
wolffd@0 95 be = be(:);
wolffd@0 96 end
wolffd@0 97
wolffd@0 98 if nargin < 7 || isempty(lb)
wolffd@0 99 error('Finance:qplcprog:MissingBoundConstraint', ...
wolffd@0 100 'Lower-bound constraint required.');
wolffd@0 101 else
wolffd@0 102 lb = full(lb);
wolffd@0 103 lb = lb(:);
wolffd@0 104 if any(~isfinite(lb))
wolffd@0 105 error('Finance:qplcprog:MissingBoundConstraint', ...
wolffd@0 106 'Finite lower-bound constraint required.');
wolffd@0 107 end
wolffd@0 108 end
wolffd@0 109
wolffd@0 110 if nargin < 8
wolffd@0 111 ub = [];
wolffd@0 112 else
wolffd@0 113 ub = full(ub);
wolffd@0 114 ub = ub(:);
wolffd@0 115 if all(ub == Inf)
wolffd@0 116 ub = [];
wolffd@0 117 end
wolffd@0 118 end
wolffd@0 119 if any(~isfinite(ub))
wolffd@0 120 error('Finance:qplcprog:InvalidBoundConstraint', ...
wolffd@0 121 'If specified, finite upper-bound constraint required.');
wolffd@0 122 end
wolffd@0 123
wolffd@0 124 % process name-value pairs (if any)
wolffd@0 125
wolffd@0 126 if nargin > 8
wolffd@0 127 if mod(nargin-8,2) ~= 0
wolffd@0 128 error('Finance:qplcprog:InvalidNameValuePair', ...
wolffd@0 129 'Invalid name-value pairs with either a missing name or value.');
wolffd@0 130 end
wolffd@0 131
wolffd@0 132 names = { 'maxiter', 'tiebreak', 'tolcon', 'tolpiv' }; % names
wolffd@0 133 values = { 10000, 'first', 1.0e-4, 1.0e-9 }; % default values
wolffd@0 134 try
wolffd@0 135 [maxiter, tiebreak, tolcon, tolpiv] = parsepvpairs(names, values, varargin{:});
wolffd@0 136
wolffd@0 137 catch E
wolffd@0 138 E.throw
wolffd@0 139 end
wolffd@0 140 else
wolffd@0 141 maxiter = 100000;
wolffd@0 142 tiebreak = 'first';
wolffd@0 143 tolpiv = 1.0e-9;
wolffd@0 144 tolcon = 1.0e-6;
wolffd@0 145 end
wolffd@0 146
wolffd@0 147 % check arguments
wolffd@0 148
wolffd@0 149 if maxiter <= 0
wolffd@0 150 error('Finance:qplcprog:NonPositiveInteger', ...
wolffd@0 151 'Maximum number of iterations (maxiter) must be a positive integer.');
wolffd@0 152 end
wolffd@0 153
wolffd@0 154 if tolcon < 2*eps
wolffd@0 155 error('Finance:qplcprog:InvalidTolerance', ...
wolffd@0 156 'Unrealistically small constraint tolerance (tolcon) specified.');
wolffd@0 157 end
wolffd@0 158
wolffd@0 159 if tolpiv < 2*eps
wolffd@0 160 error('Finance:qplcprog:InvalidTolerance', ...
wolffd@0 161 'Unrealistically small pivot tolerance (tolpiv) specified.');
wolffd@0 162 end
wolffd@0 163
wolffd@0 164 % set translation if necessary
wolffd@0 165
wolffd@0 166 if norm(lb) ~= 0
wolffd@0 167 fx = f + H*lb;
wolffd@0 168 if ~isempty(Ai) && ~isempty(bi)
wolffd@0 169 bxi = bi - Ai*lb;
wolffd@0 170 end
wolffd@0 171 if ~isempty(Ae) && ~isempty(be)
wolffd@0 172 bxe = be - Ae*lb;
wolffd@0 173 end
wolffd@0 174 if ~isempty(ub)
wolffd@0 175 uxb = ub - lb;
wolffd@0 176 end
wolffd@0 177 translate = true;
wolffd@0 178 else
wolffd@0 179 fx = f;
wolffd@0 180 if ~isempty(Ai) && ~isempty(bi)
wolffd@0 181 bxi = bi;
wolffd@0 182 end
wolffd@0 183 if ~isempty(Ae) && ~isempty(be)
wolffd@0 184 bxe = be;
wolffd@0 185 end
wolffd@0 186 if ~isempty(ub)
wolffd@0 187 uxb = ub;
wolffd@0 188 end
wolffd@0 189 translate = false;
wolffd@0 190 end
wolffd@0 191
wolffd@0 192 % set up constraints
wolffd@0 193
wolffd@0 194 if ~isempty(Ai) && ~isempty(bi)
wolffd@0 195 if ~isempty(Ae) && ~isempty(be)
wolffd@0 196 if ~isempty(ub)
wolffd@0 197 A = [ Ai; Ae; -Ae; eye(n,n) ];
wolffd@0 198 b = [ bxi; bxe; -bxe; uxb ];
wolffd@0 199 else
wolffd@0 200 A = [ Ai; Ae; -Ae ];
wolffd@0 201 b = [ bxi; bxe; -bxe ];
wolffd@0 202 end
wolffd@0 203 else
wolffd@0 204 if ~isempty(ub)
wolffd@0 205 A = [ Ai; eye(n,n) ];
wolffd@0 206 b = [ bxi; uxb ];
wolffd@0 207 else
wolffd@0 208 A = Ai;
wolffd@0 209 b = bxi;
wolffd@0 210 end
wolffd@0 211 end
wolffd@0 212 else
wolffd@0 213 if ~isempty(Ae) && ~isempty(be)
wolffd@0 214 if ~isempty(ub)
wolffd@0 215 A = [ Ae; -Ae; eye(n,n) ];
wolffd@0 216 b = [ bxe; -bxe; uxb ];
wolffd@0 217 else
wolffd@0 218 A = [ Ae; -Ae ];
wolffd@0 219 b = [ bxe; -bxe ];
wolffd@0 220 end
wolffd@0 221 else
wolffd@0 222 if ~isempty(ub)
wolffd@0 223 A = eye(n,n);
wolffd@0 224 b = uxb;
wolffd@0 225 else
wolffd@0 226 error('Finance:qplcprog:InvalidConstraints', ...
wolffd@0 227 'Insufficient or invalid constraints.');
wolffd@0 228 end
wolffd@0 229 end
wolffd@0 230 end
wolffd@0 231
wolffd@0 232 % set up lcp problem
wolffd@0 233
wolffd@0 234 ns = size(A,1);
wolffd@0 235
wolffd@0 236 M = [ zeros(ns,ns) -A; A' H ];
wolffd@0 237 q = [ b; fx ];
wolffd@0 238
wolffd@0 239 % solve lcp problem
wolffd@0 240
wolffd@0 241 if isempty(varargin)
wolffd@0 242 [z, w, status] = lcprog(M, q); %#ok
wolffd@0 243 else
wolffd@0 244 [z, w, status] = lcprog(M, q, 'maxiter', maxiter, 'tiebreak', tiebreak, 'tolpiv', tolpiv); %#ok
wolffd@0 245 end
wolffd@0 246
wolffd@0 247 % solve qp problem
wolffd@0 248
wolffd@0 249 if status > 0
wolffd@0 250 x = z(end-n+1:end);
wolffd@0 251 if translate
wolffd@0 252 x = x + lb;
wolffd@0 253 end
wolffd@0 254 fval = 0.5*x'*H*x + f'*x;
wolffd@0 255
wolffd@0 256 if qplcprogrealitycheck(Ai, bi, Ae, be, lb, ub, x)
wolffd@0 257 status = -2;
wolffd@0 258 x = NaN(n,1);
wolffd@0 259 fval = NaN;
wolffd@0 260 end
wolffd@0 261 else
wolffd@0 262 x = NaN(n,1);
wolffd@0 263 fval = NaN;
wolffd@0 264 end
wolffd@0 265
wolffd@0 266 function notok = qplcprogrealitycheck(Ai, bi, Ae, be, lb, ub, x, tolcon)
wolffd@0 267 %QPLCPROGREALITYCHECK - Make sure candidate solution x is feasible.
wolffd@0 268
wolffd@0 269 if nargin < 7
wolffd@0 270 error('Finance:qplcprog:MissingInputArgument', ...
wolffd@0 271 'Missing required input arguments Ai, bi, Ae, be, lb, ub, x.');
wolffd@0 272 end
wolffd@0 273 if nargin < 8 || isempty(tolcon)
wolffd@0 274 tolcon = 1.0e-6;
wolffd@0 275 end
wolffd@0 276
wolffd@0 277 notok = false;
wolffd@0 278
wolffd@0 279 if ~isempty(Ai)
wolffd@0 280 u = max(Ai*x - bi,0);
wolffd@0 281 if u > tolcon
wolffd@0 282 warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0 283 'Candidate solution violates inequality constraints.');
wolffd@0 284 notok = true;
wolffd@0 285 end
wolffd@0 286 end
wolffd@0 287
wolffd@0 288 if ~isempty(Ae)
wolffd@0 289 if norm(be - Ae*x,inf) > tolcon
wolffd@0 290 warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0 291 'Candidate solution violates equality constraints.');
wolffd@0 292 notok = true;
wolffd@0 293 end
wolffd@0 294 end
wolffd@0 295
wolffd@0 296 if ~isempty(lb)
wolffd@0 297 u = x - lb;
wolffd@0 298 if min(u) < -tolcon
wolffd@0 299 warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0 300 'Candidate solution violates lower-bound constraints.');
wolffd@0 301 notok = true;
wolffd@0 302 end
wolffd@0 303 end
wolffd@0 304
wolffd@0 305 if ~isempty(ub)
wolffd@0 306 u = ub - x;
wolffd@0 307 if min(u) < -tolcon
wolffd@0 308 warning('Finance:qplcprog:InfeasibleProblem', ...
wolffd@0 309 'Candidate solution violates upper-bound constraints.');
wolffd@0 310 notok = true;
wolffd@0 311 end
wolffd@0 312 end
wolffd@0 313
wolffd@0 314
wolffd@0 315 % [EOF]