new function for quadratic programming: pqpnonneg

[Jaroslav Hajek]

hi all,

I just contributed a new function: pqpnonneg (Positive Quadratic
Programming in NONNEGative variables).
It can be used to solve problems of the form

min 1/2*x'*A*x + b'*x, all (x >= 0)
where A is a positive definite matrix. The implementation exploits a
duality between least squares and positive qp problems; it is very
similar to lsqnonneg except that it solves the dual system A \ -b and
works with a Cholesky factorization instead of QR (via qrinsert &
qrdelete).

If applicable, there are two advantages against qp:

1. it is usually significantly faster
2. it will always converge in a finite number of iterations (I think
this follows from the Lawson-Hanson proof and the duality)

both of these are of course very desirable. Here's a small benchmark
creating an artificial problem and solving it via both pqpnonneg and
qp:

n = 500;
## Generate SPD matrix
for i = 1:3
  A = rand (n);
  A = A'*A + 1e-3*eye (n);
  x = rand (n, 1) - 0.5;
  b = -(A*x);
  disp ("pqpnonneg")
  tic;
  [x1, obj1] = pqpnonneg (A, b);
  toc
  norm (x1 - x)

  disp ("qp")
  tic;
  [x1, obj1] = qp (zeros (n, 1), A, b, [], [], zeros (n, 1), []);
  toc
  norm (x1 - x)
  disp ("--------------------------")
endfor

On my Core 2 Duo @ 2.83 GHz, g++ -O3 -march=native + serial GotoBLAS, I got:

pqpnonneg
Elapsed time is 0.00449514 seconds.
ans =  6.2362
qp
Elapsed time is 1.162 seconds.
ans =  6.2362
--------------------------
pqpnonneg
Elapsed time is 0.001014 seconds.
ans =  6.7428
qp
Elapsed time is 1.144 seconds.
ans =  6.7428
--------------------------
pqpnonneg
Elapsed time is 0.02081 seconds.
ans =  6.2335
qp
Elapsed time is 32.09 seconds.
ans =  6.2335
--------------------------
pqpnonneg
Elapsed time is 0.0235 seconds.
ans =  6.0290
qp
Elapsed time is 35.54 seconds.
ans =  6.0290
--------------------------
pqpnonneg
Elapsed time is 0.02568 seconds.
ans =  6.1129
qp
Elapsed time is 37.88 seconds.
ans =  6.1129
--------------------------

so, as you can see, pqpnonneg can be many times faster, as much as
1500x in this benchmark.
Btw. the reason why this goes into Octave rather than
OctaveForge/optim is the duality with lsqnonneg. If ever lsqnonneg is
moved to optim, so should pqpnonneg.

enjoy!

-- 
RNDr. Jaroslav Hajek
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz