qp() in Octave 3.0.0 returns result egregiously violating input constraints

[Gabriele Pannocchia]

Hi Josh,

> I put in the expressions you suggested into qp.m right before the
> __qp__() call.  The output I get is:
> for  norm(A*x0-b):
>    ans =  5.7220e-16
> and for 
>       all((Ain*x0) > bin -
> 1e-15)                                                                                                                                
> I got 1.  (I got zero if I put a tolerance of 1e-18)

OK.

> I have a bit more debugging information.
> So when __qp__ runs, it starts off with a feasible solution.  Then it
> varies between thinking the problem is convex and not convex.
> That is, it varies between setting info to 0 and 1, and I think this
> depends on the set of active constraints.

Loosely speaking, yes. Depending on which inequality constraints are in
the current active set, the reduced Hessian could be positive definite
or not. Your Hessian is actually indefinite, so the reduced Hessian
could be itself indefinite. 

The reduced Hessian is checked in line 240 of __qp__.cc, and pR tells us
if it is positive definite or not. If it is, i.e. info=0, the Cholesky
factors are computed and the inverse of the reduced Hessian computed. 

If instead the reduced Hessian is not positive definite, i.e. info=1,
the "most negative curvature" is found from an eigenvalue/eigenvector
decomposition of rH.
The most negative eigenvalue is evaluated in line 277, the corresponding
index is found in the subsequent lines and finally eVrH is the
eigenvector associated to the minimum eigenvalue.

By definition (as far as I know), eigenvectors should be normalized to
length 1 (i.e. the 2-norm of eVrH should be 1), so that if you report
values in the vector eVrH that exceeds 1 (in magnitude), then the
problem may be with the library that makes the eigenvalue decomposition.

Can you confirm that the 2-norm of eVrH is huge?

Thanks,
	Gabriele