Mercurial > octave-nkf

comparison scripts/optimization/qp.m @ 20298:5c088348fddb

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
qp.m: Overhaul function (fixes bug #45324). * qp.m: Put input validation first. Validate both x0 and q are vectors. Transform x0 and q to column vectors for rest of computation (bug #45324). Use double quotes instead of single to match Octave coding conventions. Re-wrap multi-line comments to break at meaningful boundaries. Use isargout to avoid calculating unnecessary outputs.
author Rik <rik@octave.org>
date Mon, 15 Jun 2015 21:56:34 +0200
parents 83792dd9bcc1
children
comparison
equal deleted inserted replaced
20293:530803d4f65f 20298:5c088348fddb
114 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup. 114 ## PKG_ADD: ## Discard result to avoid polluting workspace with ans at startup.
115 ## PKG_ADD: [~] = __all_opts__ ("qp"); 115 ## PKG_ADD: [~] = __all_opts__ ("qp");
116 116
117 function [x, obj, INFO, lambda] = qp (x0, H, varargin) 117 function [x, obj, INFO, lambda] = qp (x0, H, varargin)
118 118
119 nargs = nargin; 119 if (nargin == 1 && ischar (x0) && strcmp (x0, "defaults"))
120
121 if (nargin == 1 && ischar (x0) && strcmp (x0, 'defaults'))
122 x = optimset ("MaxIter", 200); 120 x = optimset ("MaxIter", 200);
123 return; 121 return;
124 endif 122 endif
125 123
124 nargs = nargin;
126 if (nargs > 2 && isstruct (varargin{end})) 125 if (nargs > 2 && isstruct (varargin{end}))
127 options = varargin{end}; 126 options = varargin{end};
128 nargs--; 127 nargs--;
129 else 128 else
130 options = struct (); 129 options = struct ();
130 endif
131
132 if (nargs != 2 && nargs != 3 && nargs != 5 && nargs != 7 && nargs != 10)
133 print_usage ();
131 endif 134 endif
132 135
133 if (nargs >= 3) 136 if (nargs >= 3)
134 q = varargin{1}; 137 q = varargin{1};
135 else 138 else
160 A_lb = []; 163 A_lb = [];
161 A_in = []; 164 A_in = [];
162 A_ub = []; 165 A_ub = [];
163 endif 166 endif
164 167
165 if (nargs == 2 || nargs == 3 || nargs == 5 || nargs == 7 || nargs == 10) 168 maxit = optimget (options, "MaxIter", 200);
166 169
167 maxit = optimget (options, "MaxIter", 200); 170 ## Validate the quadratic penalty.
168 171 if (! issquare (H))
169 ## Checking the quadratic penalty 172 error ("qp: quadratic penalty matrix must be square");
170 if (! issquare (H)) 173 elseif (! ishermitian (H))
171 error ("qp: quadratic penalty matrix not square"); 174 ## warning ("qp: quadratic penalty matrix not hermitian");
172 elseif (! ishermitian (H)) 175 H = (H + H')/2;
173 ## warning ("qp: quadratic penalty matrix not hermitian"); 176 endif
174 H = (H + H')/2; 177 n = rows (H);
175 endif 178
176 n = rows (H); 179 ## Validate the initial guess.
177 180 ## If empty it is resized to the right dimension and filled with 0.
178 ## Checking the initial guess (if empty it is resized to the 181 if (isempty (x0))
179 ## right dimension and filled with 0) 182 x0 = zeros (n, 1);
180 if (isempty (x0)) 183 else
181 x0 = zeros (n, 1); 184 if (! isvector (x0))
185 error ("qp: the initial guess X0 must be a vector");
182 elseif (numel (x0) != n) 186 elseif (numel (x0) != n)
183 error ("qp: the initial guess has incorrect length"); 187 error ("qp: the initial guess X0 has incorrect length");
184 endif 188 endif
185 189 x0 = x0(:); # always use column vector.
186 ## Linear penalty. 190 endif
187 if (isempty (q)) 191
188 q = zeros (n, 1); 192 ## Validate linear penalty.
193 if (isempty (q))
194 q = zeros (n, 1);
195 else
196 if (! isvector (q))
197 error ("qp: Q must be a vector");
189 elseif (numel (q) != n) 198 elseif (numel (q) != n)
190 error ("qp: the linear term has incorrect length"); 199 error ("qp: Q has incorrect length");
191 endif 200 endif
192 201 q = q(:); # always use column vector.
193 ## Equality constraint matrices 202 endif
194 if (isempty (A) || isempty (b)) 203
195 A = zeros (0, n); 204 ## Validate equality constraint matrices.
196 b = zeros (0, 1); 205 if (isempty (A) || isempty (b))
197 n_eq = 0; 206 A = zeros (0, n);
207 b = zeros (0, 1);
208 n_eq = 0;
209 else
210 [n_eq, n1] = size (A);
211 if (n1 != n)
212 error ("qp: equality constraint matrix has incorrect column dimension");
213 endif
214 if (numel (b) != n_eq)
215 error ("qp: equality constraint matrix and vector have inconsistent dimensions");
216 endif
217 endif
218
219 ## Validate bound constraints.
220 Ain = zeros (0, n);
221 bin = zeros (0, 1);
222 n_in = 0;
223 if (nargs > 5)
224 if (! isempty (lb))
225 if (numel (lb) != n)
226 error ("qp: lower bound LB has incorrect length");
227 elseif (isempty (ub))
228 Ain = [Ain; eye(n)];
229 bin = [bin; lb];
230 endif
231 endif
232
233 if (! isempty (ub))
234 if (numel (ub) != n)
235 error ("qp: upper bound UB has incorrect length");
236 elseif (isempty (lb))
237 Ain = [Ain; -eye(n)];
238 bin = [bin; -ub];
239 endif
240 endif
241
242 if (! isempty (lb) && ! isempty (ub))
243 rtol = sqrt (eps);
244 for i = 1:n
245 if (abs (lb (i) - ub(i)) < rtol*(1 + max (abs (lb(i) + ub(i)))))
246 ## These are actually an equality constraint
247 tmprow = zeros (1,n);
248 tmprow(i) = 1;
249 A = [A;tmprow];
250 b = [b; 0.5*(lb(i) + ub(i))];
251 n_eq += 1;
252 else
253 tmprow = zeros (1,n);
254 tmprow(i) = 1;
255 Ain = [Ain; tmprow; -tmprow];
256 bin = [bin; lb(i); -ub(i)];
257 n_in += 2;
258 endif
259 endfor
260 endif
261 endif
262
263 ## Validate inequality constraints.
264 if (nargs > 7)
265 [dimA_in, n1] = size (A_in);
266 if (n1 != n)
267 error ("qp: inequality constraint matrix has incorrect column dimension");
198 else 268 else
199 [n_eq, n1] = size (A); 269 if (! isempty (A_lb))
200 if (n1 != n) 270 if (numel (A_lb) != dimA_in)
201 error ("qp: equality constraint matrix has incorrect column dimension"); 271 error ("qp: inequality constraint matrix and lower bound vector are inconsistent");
202 endif 272 elseif (isempty (A_ub))
203 if (numel (b) != n_eq) 273 Ain = [Ain; A_in];
204 error ("qp: equality constraint matrix and vector have inconsistent dimension"); 274 bin = [bin; A_lb];
205 endif
206 endif
207
208 ## Bound constraints
209 Ain = zeros (0, n);
210 bin = zeros (0, 1);
211 n_in = 0;
212 if (nargs > 5)
213 if (! isempty (lb))
214 if (numel (lb) != n)
215 error ("qp: lower bound has incorrect length");
216 elseif (isempty (ub))
217 Ain = [Ain; eye(n)];
218 bin = [bin; lb];
219 endif 275 endif
220 endif 276 endif
221 277 if (! isempty (A_ub))
222 if (! isempty (ub)) 278 if (numel (A_ub) != dimA_in)
223 if (numel (ub) != n) 279 error ("qp: inequality constraint matrix and upper bound vector are inconsistent");
224 error ("qp: upper bound has incorrect length"); 280 elseif (isempty (A_lb))
225 elseif (isempty (lb)) 281 Ain = [Ain; -A_in];
226 Ain = [Ain; -eye(n)]; 282 bin = [bin; -A_ub];
227 bin = [bin; -ub];
228 endif 283 endif
229 endif 284 endif
230 285
231 if (! isempty (lb) && ! isempty (ub)) 286 if (! isempty (A_lb) && ! isempty (A_ub))
232 rtol = sqrt (eps); 287 rtol = sqrt (eps);
233 for i = 1:n 288 for i = 1:dimA_in
234 if (abs (lb (i) - ub(i)) < rtol*(1 + max (abs (lb(i) + ub(i))))) 289 if (abs (A_lb(i) - A_ub(i))
290 < rtol*(1 + max (abs (A_lb(i) + A_ub(i)))))
235 ## These are actually an equality constraint 291 ## These are actually an equality constraint
236 tmprow = zeros (1,n); 292 tmprow = A_in(i,:);
237 tmprow(i) = 1;
238 A = [A;tmprow]; 293 A = [A;tmprow];
239 b = [b; 0.5*(lb(i) + ub(i))]; 294 b = [b; 0.5*(A_lb(i) + A_ub(i))];
240 n_eq += 1; 295 n_eq += 1;
241 else 296 else
242 tmprow = zeros (1,n); 297 tmprow = A_in(i,:);
243 tmprow(i) = 1;
244 Ain = [Ain; tmprow; -tmprow]; 298 Ain = [Ain; tmprow; -tmprow];
245 bin = [bin; lb(i); -ub(i)]; 299 bin = [bin; A_lb(i); -A_ub(i)];
246 n_in += 2; 300 n_in += 2;
247 endif 301 endif
248 endfor 302 endfor
249 endif 303 endif
250 endif 304 endif
251 305 endif
252 ## Inequality constraints 306
253 if (nargs > 7) 307 ## Now we should have the following QP:
254 [dimA_in, n1] = size (A_in); 308 ##
255 if (n1 != n) 309 ## min_x 0.5*x'*H*x + x'*q
256 error ("qp: inequality constraint matrix has incorrect column dimension"); 310 ## s.t. A*x = b
257 else 311 ## Ain*x >= bin
258 if (! isempty (A_lb)) 312
259 if (numel (A_lb) != dimA_in) 313 ## Discard inequality constraints that have -Inf bounds since those
260 error ("qp: inequality constraint matrix and lower bound vector inconsistent"); 314 ## will never be active.
261 elseif (isempty (A_ub)) 315 idx = (bin == -Inf);
262 Ain = [Ain; A_in]; 316
263 bin = [bin; A_lb]; 317 bin(idx) = [];
318 Ain(idx,:) = [];
319
320 n_in = numel (bin);
321
322 ## Check if the initial guess is feasible.
323 if (isa (x0, "single") || isa (H, "single") || isa (q, "single")
324 || isa (A, "single") || isa (b, "single"))
325 rtol = sqrt (eps ("single"));
326 else
327 rtol = sqrt (eps);
328 endif
329
330 eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+abs (b)));
331 in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+abs (bin))));
332
333 info = 0;
334 if (eq_infeasible || in_infeasible)
335 ## The initial guess is not feasible.
336 ## First, define an xbar that is feasible with respect to the
337 ## equality constraints.
338 if (eq_infeasible)
339 if (rank (A) < n_eq)
340 error ("qp: equality constraint matrix must be full row rank");
341 endif
342 xbar = pinv (A) * b;
343 else
344 xbar = x0;
345 endif
346
347 ## Second, check that xbar is also feasible with respect to the
348 ## inequality constraints.
349 if (n_in > 0)
350 res = Ain * xbar - bin;
351 if (any (res < -rtol * (1 + abs (bin))))
352 ## xbar is not feasible with respect to the inequality constraints.
353 ## Compute a step in the null space of the equality constraints,
354 ## by solving a QP. If the slack is small, we have a feasible initial
355 ## guess. Otherwise, the problem is infeasible.
356 if (n_eq > 0)
357 Z = null (A);
358 if (isempty (Z))
359 ## The problem is infeasible because A is square and full rank,
360 ## but xbar is not feasible.
361 info = 6;
264 endif 362 endif
265 endif 363 endif
266 if (! isempty (A_ub)) 364
267 if (numel (A_ub) != dimA_in) 365 if (info != 6)
268 error ("qp: inequality constraint matrix and upper bound vector inconsistent"); 366 ## Solve an LP with additional slack variables
269 elseif (isempty (A_lb)) 367 ## to find a feasible starting point.
270 Ain = [Ain; -A_in]; 368 gamma = eye (n_in);
271 bin = [bin; -A_ub]; 369 if (n_eq > 0)
370 Atmp = [Ain*Z, gamma];
371 btmp = -res;
372 else
373 Atmp = [Ain, gamma];
374 btmp = bin;
272 endif 375 endif
273 endif 376 ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];
274 377 lb = [-Inf(n-n_eq,1); zeros(n_in,1)];
275 if (! isempty (A_lb) && ! isempty (A_ub)) 378 ub = [];
276 rtol = sqrt (eps); 379 ctype = repmat ("L", n_in, 1);
277 for i = 1:dimA_in 380 [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
278 if (abs (A_lb(i) - A_ub(i)) 381 if ((status == 0)
279 < rtol*(1 + max (abs (A_lb(i) + A_ub(i))))) 382 && all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp))))
280 ## These are actually an equality constraint 383 ## We found a feasible starting point
281 tmprow = A_in(i,:); 384 if (n_eq > 0)
282 A = [A;tmprow]; 385 x0 = xbar + Z*P(1:n-n_eq);
283 b = [b; 0.5*(A_lb(i) + A_ub(i))];
284 n_eq += 1;
285 else 386 else
286 tmprow = A_in(i,:); 387 x0 = P(1:n);
287 Ain = [Ain; tmprow; -tmprow];
288 bin = [bin; A_lb(i); -A_ub(i)];
289 n_in += 2;
290 endif 388 endif
291 endfor 389 else
292 endif 390 ## The problem is infeasible
293 endif 391 info = 6;
294 endif
295
296 ## Now we should have the following QP:
297 ##
298 ## min_x 0.5*x'*H*x + x'*q
299 ## s.t. A*x = b
300 ## Ain*x >= bin
301
302 ## Discard inequality constraints that have -Inf bounds since those
303 ## will never be active.
304 idx = isinf (bin) & bin < 0;
305
306 bin(idx) = [];
307 Ain(idx,:) = [];
308
309 n_in = numel (bin);
310
311 ## Check if the initial guess is feasible.
312 if (isa (x0, "single") || isa (H, "single") || isa (q, "single")
313 || isa (A, "single") || isa (b, "single"))
314 rtol = sqrt (eps ("single"));
315 else
316 rtol = sqrt (eps);
317 endif
318
319 eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+abs (b)));
320 in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+abs (bin))));
321
322 info = 0;
323 if (eq_infeasible || in_infeasible)
324 ## The initial guess is not feasible.
325 ## First define xbar that is feasible with respect to the equality
326 ## constraints.
327 if (eq_infeasible)
328 if (rank (A) < n_eq)
329 error ("qp: equality constraint matrix must be full row rank");
330 endif
331 xbar = pinv (A) * b;
332 else
333 xbar = x0;
334 endif
335
336 ## Check if xbar is feasible with respect to the inequality
337 ## constraints also.
338 if (n_in > 0)
339 res = Ain * xbar - bin;
340 if (any (res < -rtol * (1 + abs (bin))))
341 ## xbar is not feasible with respect to the inequality
342 ## constraints. Compute a step in the null space of the
343 ## equality constraints, by solving a QP. If the slack is
344 ## small, we have a feasible initial guess. Otherwise, the
345 ## problem is infeasible.
346 if (n_eq > 0)
347 Z = null (A);
348 if (isempty (Z))
349 ## The problem is infeasible because A is square and full
350 ## rank, but xbar is not feasible.
351 info = 6;
352 endif
353 endif 392 endif
354
355 if (info != 6)
356 ## Solve an LP with additional slack variables to find
357 ## a feasible starting point.
358 gamma = eye (n_in);
359 if (n_eq > 0)
360 Atmp = [Ain*Z, gamma];
361 btmp = -res;
362 else
363 Atmp = [Ain, gamma];
364 btmp = bin;
365 endif
366 ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];
367 lb = [-Inf(n-n_eq,1); zeros(n_in,1)];
368 ub = [];
369 ctype = repmat ("L", n_in, 1);
370 [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
371 if ((status == 0)
372 && all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp))))
373 ## We found a feasible starting point
374 if (n_eq > 0)
375 x0 = xbar + Z*P(1:n-n_eq);
376 else
377 x0 = P(1:n);
378 endif
379 else
380 ## The problem is infeasible
381 info = 6;
382 endif
383 endif
384 else
385 ## xbar is feasible. We use it a starting point.
386 x0 = xbar;
387 endif 393 endif
388 else 394 else
389 ## xbar is feasible. We use it a starting point. 395 ## xbar is feasible. We use it a starting point.
390 x0 = xbar; 396 x0 = xbar;
391 endif 397 endif
392 endif
393
394 if (info == 0)
395 ## The initial (or computed) guess is feasible.
396 ## We call the solver.
397 [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);
398 else 398 else
399 iter = 0; 399 ## xbar is feasible. We use it a starting point.
400 x = x0; 400 x0 = xbar;
401 lambda = []; 401 endif
402 endif 402 endif
403
404 if (info == 0)
405 ## The initial (or computed) guess is feasible. Call the solver.
406 [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);
407 else
408 iter = 0;
409 x = x0;
410 lambda = [];
411 endif
412 if (isargout (2))
403 obj = 0.5 * x' * H * x + q' * x; 413 obj = 0.5 * x' * H * x + q' * x;
414 endif
415 if (isargout (3))
404 INFO.solveiter = iter; 416 INFO.solveiter = iter;
405 INFO.info = info; 417 INFO.info = info;
406
407 else
408 print_usage ();
409 endif 418 endif
410 419
411 endfunction 420 endfunction
412 421
413 422