Julien, I didn't say that there was no stopping criterion in the code.  I was just concerned that the criterion may never be reached because the numerical precision keeps stopping chol from succeeding.

Perhaps we can discuss this more after 4.2.0 is out.

Retagging from "Patch submitted" to "In progress".

Nice ideas, Julien.

Are these "robust" in the sense of numerical stability, or just in the sense of always getting an answer?  In the case of comment #4, the error was around 10^-15, but I don't see anything in the code that will stop when numerical precision is exhausted.

If they can be made numerically stable, perhaps they should be put in   chol   itself, rather than as a hack in _qp_.  What do you think?

For the numerical example of comment #4, it turns out that a decent Cholesky factorization can be obtained in two steps:

load qp_data.mat

% Find the largest eigenvalue (can be generalized to several)
[V, D] = eigs (H, 1, 'LM');

% Deflate
W = null (V');
H1 = W' * H * W;

% Now we can factorize using chol
R1 = chol (H1);

% Reconstruct a Cholesky factor for H
R = (sqrt (D)) * V * V' + W * R1 * W';

% Reconstruction
Hrec = R' * R;

% Error
err_rel = (H - Hrec) ./ (abs (H));
max (err_rel(:))

(I get a maximal relative error of approx 1.7e-15)

Octave 4.0.0 still displays this behavior with the qp test case from comment #4 yielding:

error: __qp__: operator *: nonconformant arguments (op1 is 6x6, op2 is 7x1)

Hello Michele,

All I can do right now is confirm that the problem still exists for me in Octave 3.8.2-RC2 (Ubuntu 13.04, 32bits, gcc 4.7.3).

@++
Julien

Here is a snippet that reproduces the _qp_ error:

clear all;
load qp_data;
maxit = 200;

[x, lambda, info, iter] = __qp__ ...
  (x0, H, q, [], [], Ain, bin, maxit);

(file #29590)

The same problem (or so it seems) is discussed here:

https://savannah.gnu.org/bugs/?32717

Hello,

I wrote a program grenzlinietest2.m which uses sqp and produces strange errors.

If I use octave version 3.2.4, everything works fine. After we updated to version 3.6.4 and installed the package optim, we received the following error

warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
warning: sqp: QP subproblem failed to converge in 200 iterations
error: _qp_: operator *: nonconformant arguments (op1 is 6x6, op2 is 7x1)
error: called from:
error:   /usr/local/share/octave/3.6.4/m/optimization/qp.m at line 393, column 26
error:   /usr/local/share/octave/3.6.4/m/optimization/sqp.m at line 414, column 32
error:   /home/hackstein/Dominantepole5/minimize6test2.m at line 9, column 5
error:   /home/hackstein/Dominantepole5/verschiebungize3test2.m at line 8, column 3
error:   /home/hackstein/Dominantepole5/grenzlinietest2.m at line 10, column 7

For the functions, see the attachment. As I can attach only four functions, I post three of them here:

function erg=funp6test2(x,p)
erg=(p(1)-real(dominant2(roots(polytest2(x(1),x(2),x(3),x(4),x(5),x(6),x(7)))))).^2+(p(2)-imag(dominant2(roots(polytest2(x(1),x(2),x(3),x(4),x(5),x(6),x(7)))))).^2
endfunction


function mini=minimize6test2(b)
funp6test2wrap = @(x)funp6test2(x,[real(b),imag(b)])
y0  = 2*ones(7,1);
lbound = [0.5,0.5,0.5,0.5,0.5,0.5,5]';
ubound = [20,20,20,20,20,20,200]';
mini=sqp (y0,funp6test2wrap,[],[],lbound,ubound)
endfunction

function poly=polytest2(la, lb, lc, Ia, Ib, Ic, RC)
poly=[1 la+lb*RC*45*la.^7 lb+lc-RC-3*Ib*Ic Ia*RC+Ib+87*RC Ib*Ic.^3+Ic/(Ia+Ic-la) RC-Ia*Ia*lc la*RC*Ib+Ia*RC koeff0(la,lb,lc,Ia,Ib,Ic,RC)]
endfunction

Minimize6test2.m minimizes the function
(real(b)-real(dominant2(roots(polytest2(x(1),x(2),x(3),x(4),x(5),x(6),x(7)))))).^2+(imag(b)-imag(dominant2(roots(polytest2(x(1),x(2),x(3),x(4),x(5),x(6),x(7)))))).^2

Here b is a complex number and polytest2(x(1),x(2),x(3),x(4),x(5),x(6),x(7)) is a polynomial of degree 7, parametrized by the seven real parameters x(1),x(2),x(3),x(4),x(5),x(6) and x(7).
dominant2 selects the root of the polynomial with the biggest real part. The term should be minimized with respect to 0.5<=x(1),x(2),x(3),x(4),x(5),x(6)<=20 and 5<=x(7)<=200.

The error above is produced by the function call grenzlinietest2(2,4,6,1,2,6,50). Note that more strange errors occur if one takes a polynomial more complicated than polytest2. See the discussion on help-octave for more details.
Thanks a lot in advance!