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

 
 

bug #38875: strange errors using sqp

Submitted by:  Urs Hackstein <uhackstein>
Submitted on:  Fri May 3 10:53:49 2013  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Other
Status: In ProgressAssigned to: None
Originator Name: Open/Closed: Open
Release: devOperating System: Any

Add a New Comment (Rich MarkupRich Markup):
   

You are not logged in

Please log in, so followups can be emailed to you.

 

(Jump to the original submission Jump to the original submission)

Thu Jun 30 13:29:03 2016, 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>
Project Member
Tue Mar 1 06:16:07 2016, 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 Mar 1 01:13:03 2016, 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>
Project Member
Mon Feb 29 13:27:28 2016, comment #11:

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

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

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

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

Julien Bect <jbect>
Sat Feb 27 10:37:19 2016, 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>
Project Member
Thu Aug 6 09:33:51 2015, comment #8:

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

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

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

Anonymous
Wed Jul 16 20:09:05 2014, 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 Jul 6 08:53:26 2014, comment #5:

I am experiencing the same problem on Octave 3.8.1.

Any ideas? Thanks a lot!

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

Here is a snippet that reproduces the _qp_ error:

(file #29590)

Julien Bect <jbect>
Sun Nov 10 21:50:53 2013, 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 Nov 8 21:15:07 2013, comment #2:

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

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

Julien Bect <jbect>
Fri Nov 8 21:08:13 2013, comment #1:

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

Julien Bect <jbect>
Fri May 3 10:53:49 2013, 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+lbRC45la.^7 lb+lc-RC-3IbIc IaRC+Ib+87RC IbIc.^3+Ic/(Ia+Ic-la) RC-IaIalc laRCIb+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 File(s):
   
   
Comment:
   

Attached Files
file #36489:  bug_38875_qp.cset added by lachlan (889B - text/x-diff)
file #29590:  qp_data.mat added by jbect (2kB - 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 (2kB - text/x-matlab)
file #28005:  dominant2.m added by uhackstein (2kB - text/x-matlab)
file #28006:  koeff0.m added by uhackstein (1kB - text/x-matlab)

 

Depends on the following items: None found

Items that depend on this one: None found

 

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

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 15 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Thu Jun 30 13:29:03 2016lachlanStatusPatch Submitted=>In Progress
    Sat Feb 27 10:37:19 2016lachlanAttached File-=>Added bug_38875_qp.cset, #36489
      StatusNone=>Patch Submitted
      Release3.6.4=>dev
      Operating SystemGNU/Linux=>Any
    Sun Nov 10 22:05:44 2013jbectAttached File-=>Added qp_data.mat, #29590
    Mon Oct 28 02:14:53 2013mtmillerCategoryNone=>Octave Function
    Fri May 3 10:57:55 2013uhacksteinAttached File-=>Added polytest2.m, #28007
      Attached File-=>Added funp6test2.m, #28008
      Attached File-=>Added minimize6test2.m, #28009
    Fri May 3 10:53:49 2013uhacksteinAttached 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 -unavailable-

    Back to the top


    Powered by Savane 3.1-cleanup