bugGNU Octave - Bugs: bug #38875, strange errors using sqp

 
 

bug #38875: strange errors using sqp

Submitter:  Urs Hackstein <uhackstein>
Submitted:  Fri 03 May 2013 10:53:49 AM UTC
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Other
Status:  In Progress Assigned to:  None
Originator Name:  Open/Closed:  * Open
Release:  * dev Operating System:  * Any
* Mandatory Fields

Add a New Comment (Rich Markup)
   

Comment Type & Canned Response:
       

 

( Jump to the original submission )

Thu 30 Jun 2016 01:29:03 PM UTC, comment #14: 

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".

Lachlan Andrew <lachlan>
Tue 01 Mar 2016 06:16:07 AM UTC, comment #13: 

Robust in the sense of getting an answer more often that's fore sure

How "reasonable" is that answer? Can we claim that it will more robust in terms of numerical stability? I can't tell. To answer that theoretically you would have to know what's under the hood of eigs and null (and be much more knowledgeable in numerical analysis than I am).

Another approach would be to set up a large-scale numerical benchmark on Cholesky factorization and then (hopefully) demonstrate empirically that this trick improves the results.

I would recommend the following:  a) First see if this trick (insert somewhere in sqp or qp) solves the original problem identified in this bug report. I can prepare a patch for that. If it does solve the problem, apply it.  b) Create a new bug report clearly focused on trying to improve chol for ill-conditionned matrices.

Finally, there is a stopping condition in my code.  Recursion stops as soon as a plain chol succeeds.  For nice matrices, there is no iteration at all.

Julien Bect <jbect>
Tue 01 Mar 2016 01:13:03 AM UTC, comment #12: 

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?

Lachlan Andrew <lachlan>
Mon 29 Feb 2016 01:27:28 PM UTC, comment #11: 

We could using the trick that I just proposed in recursive manner, as follows:

function R = robust_chol (H)

n = size (H, 1);  fprintf ('n = %d  -->  ', n);

[R, p] = chol (H);

if p == 0

    fprintf ('OK\n');

else

    fprintf ('chol failed, deflating...\n');

    % 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 = robust_chol (H1);

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

end

end

Julien Bect <jbect>
Mon 29 Feb 2016 01:13:09 PM UTC, comment #10: 

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)

Julien Bect <jbect>
Sat 27 Feb 2016 10:37:19 AM UTC, comment #9: 

The problem that Julien identified in comment #4 is because the range of eigenvalues of the Hessian matrix is too large, and so the inverse computed by the Cholesky decomposition is of lower dimension than the Hessian itself.

The attached patch throws a more informative error ("Hessian must be positive definite") in that case.  It may be slightly misleading since H is actually positive definite, but I wanted to match the error message of chol.  Thoughts?

The 0x0 matrix in comment #1 is likely to have a different origin, like that in bug #36015.

(file #36489)

Lachlan Andrew <lachlan>
Thu 06 Aug 2015 09:33:51 AM UTC, comment #8: 

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)

Reinhard Prix <repr>
Fri 03 Oct 2014 12:50:23 PM UTC, comment #7: 

Same problem here both under 3.81 (ubuntu 14.04 LTS) and 3.6.4 (Ubuntu 13.04)

Anonymous
Wed 16 Jul 2014 08:09:05 PM UTC, comment #6: 

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

Julien Bect <jbect>
Sun 06 Jul 2014 08:53:26 AM UTC, comment #5: 

I am experiencing the same problem on Octave 3.8.1.

Any ideas? Thanks a lot!

Michele Dall'Arno <mdallarno>
Sun 10 Nov 2013 10:05:44 PM UTC, comment #4: 

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)

Julien Bect <jbect>
Sun 10 Nov 2013 09:50:53 PM UTC, comment #3: 

Mea culpa: this is not the same as bug #32717.

After applying the patch proposed here:

https://savannah.gnu.org/bugs/index.php?40536

the error that was originally described in this bug report re-appears. It seems to be a problem in _qp_.

Julien Bect <jbect>
Fri 08 Nov 2013 09:15:07 PM UTC, comment #2: 

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

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

Julien Bect <jbect>
Fri 08 Nov 2013 09:08:13 PM UTC, comment #1: 

The proposed example still fails in Octave 3.7.7+ (3851e5fde76d), but the error message has changed:

error: sqp: operator *: nonconformant arguments (op1 is 7x14, op2 is 0x0)
error: called from:
error:   /media/data/sources/octave/core/scripts/optimization/sqp.m at line 450, column 9
error:   /home/bect/Sandbox/matlabOctave/trackSQPbug/minimize6test2.m at line 9, column 5
error:   /home/bect/Sandbox/matlabOctave/trackSQPbug/verschiebungize3test2.m at line 8, column 3
error:   /home/bect/Sandbox/matlabOctave/trackSQPbug/grenzlinietest2.m at line 4, column 7

Julien Bect <jbect>
Fri 03 May 2013 10:53:49 AM UTC, original submission:  

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!

Urs Hackstein <uhackstein>

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #36489:  bug_38875_qp.cset added by lachlan (889B - text/x-diff)
file #29590:  qp_data.mat added by jbect (2KiB - application/octet-stream)
file #28007:  polytest2.m added by uhackstein (201B - text/x-matlab - Here are the three other functions)
file #28008:  funp6test2.m added by uhackstein (204B - text/x-matlab - Here are the three other functions)
file #28009:  minimize6test2.m added by uhackstein (308B - text/x-matlab - Here are the three other functions)
file #28003:  grenzlinietest2.m added by uhackstein (636B - text/x-matlab)
file #28004:  verschiebungize3test2.m added by uhackstein (2KiB - text/x-matlab)
file #28005:  dominant2.m added by uhackstein (2KiB - text/x-matlab)
file #28006:  koeff0.m added by uhackstein (1KiB - text/x-matlab)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by lachlan (Updated the item)
  • -email is unavailable- added by repr (Posted a comment)
  • -email is unavailable- added by mdallarno (Posted a comment)
  • -email is unavailable- added by jbect (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by uhackstein (Submitted the item)
  • -email is unavailable- added by uhackstein
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only project members can vote.

     

     

     

     

    Follow 15 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-06-30 lachlan StatusPatch Submitted In Progress
    2016-02-27 lachlan Attached File- Added bug_38875_qp.cset, #36489
        StatusNone Patch Submitted
        Release3.6.4 dev
        Operating SystemGNU/Linux Any
    2013-11-10 jbect Attached File- Added qp_data.mat, #29590
    2013-10-28 mtmiller CategoryNone Octave Function
    2013-05-03 uhackstein Attached File- Added polytest2.m, #28007
        Attached File- Added funp6test2.m, #28008
        Attached File- Added minimize6test2.m, #28009
    2013-05-03 uhackstein Attached File- Added grenzlinietest2.m, #28003
        Attached File- Added verschiebungize3test2.m, #28004
        Attached File- Added dominant2.m, #28005
        Attached File- Added koeff0.m, #28006
        Carbon-Copy- Added -email is unavailable-

    Back to the top


    Powered by Savane 3.9