bugGNU Octave - Bugs: bug #53991, Octave crashes using fsolve when...

 
 

bug #53991: Octave crashes using fsolve when x0 is a stationary point

Submitter:  None
Submitted:  Sun 27 May 2018 11:23:35 AM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Unexpected Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Stefano Toro Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 4.4.0
Operating System:  * Microsoft Windows Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Mon 25 Jun 2018 08:53:36 PM UTC, comment #10: 

I backported the patch to exit with info=-2 when a singularity is reached (http://hg.savannah.gnu.org/hgweb/octave/rev/68e2daa29222).

I added a BIST test on stable to verify that the patch above works correctly (http://hg.savannah.gnu.org/hgweb/octave/rev/c29a2107c559).

Finally, I checked in a third patch to return a success code (info=1) if the initial guess for the solution is already accurate enough (http://hg.savannah.gnu.org/hgweb/octave/rev/b3c35a130f94).  I added a BIST test for that too.

Marking bug as fixed and closing report.

Rik <rik5>
Group administrator
Thu 21 Jun 2018 02:00:16 PM UTC, comment #9: 

I checked in a changeset to address the problem of bug #39000, mentioned here in comment #4.  See that report for details.

John W. Eaton <jwe>
Group administrator
Thu 21 Jun 2018 08:55:34 AM UTC, comment #8: 

What about something like this


-        if (delta <= 1e1*macheps*xn)
+        if (fn <= tolf*n*xn)
+          info = 1;
+        elseif (delta <= 1e1*macheps*xn)


Marco Caliari <caliari>
Group Member
Wed 20 Jun 2018 06:35:33 PM UTC, comment #7: 

I checked in a change to avoid the crash at a stationary point (http://hg.savannah.gnu.org/hgweb/octave/rev/dffd9f6ee85c).  I used this test


+        if (norm (fjac, 1) < macheps * rows (fjac))
+          info = -2;
+          break;
+        endif


Since norm (x,1) is an absolute column sum I needed to scale macheps by the number of rows in the matrix.

I will probably backport this to stable, but not before the other issue is solved.  You are correct that if the user happens to specify a solution as the initial starting point then fsolve returns info=-3.  A simple test is just a line,


[x, fvec, info] = fsolve (@(x) 5*x, 0)
x = 0
fvec = 0
info = -3


The simplest thing would just be to copy the test for success to be at the top of the file just outside the start of the outer loop.


      ## FIXME: Why tolf*n*xn? If abs (e) ~ abs(x) * eps is a vector
      ## of perturbations of x, then norm (fjac*e) <= eps*n*xn, i.e., by
      ## tolf ~ eps we demand as much accuracy as we can expect.
      if (fn <= tolf*n*xn)
        info = 1;


However, we don't have xn available at that time.  Maybe someone can work out what xn should be for the first case of iteration.

Rik <rik5>
Group administrator
Wed 20 Jun 2018 08:46:08 AM UTC, comment #6: 

@Rik: my first thought was that info = -2 should be set when fjac is small AND fvec is still large. So I tried


[~, ~, flag] = fsolve (detB, pi)


(pi is the right solution), but I get flag = -3. I would conclude that:

1) fsolve should check whether the initial solution is already fine before the ## Outer loop

2) when it computes a new fjac, it means that the previous fvec was not fine. Therefore it is correct just to check whether the new fjac is not too small, as you did

3) I would use norm (fjac, 1). The 1-norm is easier to compute.

4) For unconverged cases, Matlab gives a message (not a warning).

Marco Caliari <caliari>
Group Member
Tue 19 Jun 2018 09:44:01 PM UTC, comment #5: 

@Marco: Can you review this patch for fsolve?


diff -r b7db401e1a99 scripts/optimization/fsolve.m
--- a/scripts/optimization/fsolve.m        Tue Jun 19 09:18:44 2018 -0700
+++ b/scripts/optimization/fsolve.m        Tue Jun 19 14:23:19 2018 -0700
@@ -304,10 +304,18 @@ function [x, fvec, info, output, fjac] =

       ## Get trust-region model (dogleg) minimizer.
       if (useqr)
+        if (abs (r) < macheps)
+          info = -2;
+          break;
+        endif
         qtf = q'*fvec;
         s = - __dogleg__ (r, qtf, dg, delta);
         w = qtf + r * s;
       else
+        if (abs (fjac) < macheps)
+          info = -2;
+          break;
+        endif
         s = - __dogleg__ (fjac, fvec, dg, delta);
         w = fvec + fjac * s;
       endif


If the Jacobian is very small this stops fsolve before calling _dogleg_ with a zero value for the first argument which would eventually lead to an ill-conditioned matrix.

One question, should the test be for any component of fjac or should we be taking the norm first?  In other words, this


if (any (abs (fjac(:)) < macheps))
  info = -2;
  break;
endif


OR


if (norm (fjac) < macheps)
  info = -2;
  break;
endif


I expect it is something like the latter.  The first thing that the __dogleg__routine does is calculate


x = fjac \ fvec;


If fjac is a vector or matrix having a zero would be okay.  Only an all zero array would be a problem.

Rik <rik5>
Group administrator
Wed 13 Jun 2018 07:14:27 PM UTC, comment #4: 

I think there are two separate bug reports here.  The first is that calling balance can lead to a segfault on Windows platforms.

Using the test code from comment #1 on a Linux machine I get


 ** On entry to DGEBAL parameter number  3 had an illegal value
A =

     0     1     0     0
     0     0     1     0
     0     0     0     1
     0     0   NaN     0

error: balance: unknown error in fortran subroutine


The segfault is probably back to the inability to overide the XERBLA handler on Windows (bug #39000).  The question is how far we want to keep patching things.  It seems like possibly every routine which ever interacts with the BLAS or LAPACK library would potentially need to look for NaN values before continuing.  At the very least, it would seem that this should be restricted to Windows platforms which would require putting in #ifdef tests.

The second bug report is about fsolve not returning a correct error code when it encounters a stationary point.  This should be fixed, and it shouldn't be that difficult.

Rik <rik5>
Group administrator
Mon 28 May 2018 08:55:06 AM UTC, comment #3: 

I think the problem is in the algorithm itself. Matlab returns with the exit flag -2 which means


Converged to a point that is not a root.


This happens when the Jacobian fjac has a small norm (0, in this report) but the solution fvec is far from 0. This can be detected before calling _dogleg_ at line 307.

Marco Caliari <caliari>
Group Member
Sun 27 May 2018 08:51:29 PM UTC, comment #2: 

The problem is in the call to DGEBAL in liboctave/numeric/aepbalance.cc

Avinoam Kalma <avinoam>
Group Member
Sun 27 May 2018 08:38:49 PM UTC, comment #1: 

The line that cause the crash is


A = [0     1     0     0;
     0     0     1     0;
     0     0     0     1;
     0     0   NaN     0]
balance(A);


The bug seems to be related to bug #39000

Avinoam Kalma <avinoam>
Group Member
Sun 27 May 2018 11:23:35 AM UTC, original submission:  

Trying to solve equation det (B(lam)) = 0 for lam, if starting point lam0 = 0, Octave crashes without error or warning.
If lam0 = 0.1, solution found correctly.

Here the script, with a plot of B vs. lam:


mu = 0;

A = @(lam) [0 1 0 0; 0 0 1 0; 0 0 0 1; mu^2 0 -lam^2 0];
C = [1 0 0 0; 0 0 1 0];
B = @(lam) [C*expm(A(lam)*0); C*expm(A(lam)*1)];
detB = @(lam) det(B(lam));

lam0 = 0;
lcrit = fsolve (detB, lam0);

lambda = 0:0.01:4*pi;
for j=1:length(lambda)
  d(j) = detB(lambda(j));
end

plot(lambda,d)
grid on


Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #44228:  detB_vs_lam.pdf added by None (188KiB - application/pdf)
file #44227:  B_vs_lam.ps added by None (185KiB - application/postscript)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by avinoam (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    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 group members can vote.

     

    Follow 5 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-06-25 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2018-05-27 avinoam StatusNone Confirmed
    2018-05-27 None Attached File- Added detB_vs_lam.pdf, #44228
    2018-05-27 None Attached File- Added B_vs_lam.ps, #44227

    Back to the top

    Powered by Savane 3.13-3230.
    Corresponding source code