bugGNU Octave - Bugs: bug #39879, expm returns inaccurate results


bug #39879: expm returns inaccurate results

Submitted by:  None
Submitted on:  Tue 27 Aug 2013 01:57:48 PM UTC  
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  vathomas Originator Email:  -email is unavailable-
Open/Closed:  Closed Release:  dev
Operating System:  GNU/Linux

Add a New Comment(Rich Markup)

You are not logged in

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


( Jump to the original submission)

Wed 28 Mar 2018 04:08:14 AM UTC, comment #10:

I agree that this bug seems fixed now with current versions of LAPACK or OpenBLAS. I'm closing it as fixed. If there is still something that needs to be done in Octave, please comment and/or reopen this bug, or file a new bug.

Mike Miller <mtmiller>
Project Administrator
Tue 22 Nov 2016 09:57:28 PM UTC, comment #9:

This issue seems to be fixed in Octave 4.2.0 (unless you use a lapack version <3.5.0).

Running the sample code from comment #0 in Octave (on my not so new Ubuntu 14.04 with lapack 3.5.0) gives me an error of only 3.9540e-15. As said in comment #7 this is probably due to my (high enough) lapack version. Older lapack versions probably still have this issue. But 3.5.0 is now (as time has passed by) not very new any more, current version is 3.6.1 and mxe-octave is also built with 3.5.0 today.

I have not checked comment #8.

Is it too early to close this issue? This was always a lapack bug, no Octave bug, when I understand the discussion right. And this bug has been fixed in lapack more than 2.5 years ago, now.

Hartmut <hardy>
Tue 15 Dec 2015 06:07:15 PM UTC, comment #8:

Here is the case where expm from GNU Octave fails to
give the correct results:
A = [0 1e-8 0;-(2e10+4e8/6) -3 2e10;200/3 0 -200/3];

It becomes clear that when we vary the order of magnitude
of the input matrix elements the accuracy of
the results vary.

Mudit Sharma <mudit_sharma>
Thu 05 Jun 2014 11:59:24 AM UTC, comment #7:


I tested LAPACK from 3.4.1 to 3.5.0 (fortran version) with tmpa matrix and the problem is fixed in 3.5.0. So, I presume you have 3.5.0.
About Octave:

*) balance in LAPACK < 3.5.0 has a problem, therefore any Octave version linked to those libraries has the same problem
*) expm works even without balancing. Shall we remove the balancing/unbalancing steps? Shall we make a comment in expm.m?


Marco Caliari <caliari>
Project Member
Thu 05 Jun 2014 10:31:28 AM UTC, comment #6:

I do not see much in octave_config_info (except from -llapack flag that does not help I think).

Running your code gives the following answer:

ans = 0.25000
ans = 0.12500


Thu 05 Jun 2014 10:00:03 AM UTC, comment #5:

I don't know how to check lapack version (maybe octave_config_info?).
Anyway, can you post the result of

I get


Marco Caliari <caliari>
Project Member
Thu 05 Jun 2014 08:40:02 AM UTC, comment #4:

I am not sure how to tell which LAPACK I am using. But I run the test with the MXE for windows (3.8.1-2 installer from http://mxeoctave.osuv.de/) and the norm has the value

The norm had more or less the same value when I tested on my debian machine. There, I use the current "testing" package 3.8.1 which I think is using ATLAS.

If you tell me how I can provide more info, I would be happy to.


Thu 05 Jun 2014 07:47:46 AM UTC, comment #3:


the bug is not solved in 3.8.0. I still see the difference 0.0953. Can you please tell which LAPACK are you using?



Marco Caliari <caliari>
Project Member
Thu 05 Jun 2014 07:34:43 AM UTC, comment #2:

This bug is solved in version 3.8. Could you please close it.


Thu 29 Aug 2013 12:08:18 PM UTC, comment #1:

The problem is in the balancing of the matrix.
For instance, the first row of the "balanced" matrix has 1-norm 5e-2 and the first column 9e14 (for the original matrix they are 3.5e-1 and 1e-1).
I think balance uses xGEBAL and that the problem is there (I'm using LAPACK 3.4.1 ang getting similar results with Matlab, where it is written "Balancing can destroy the properties of certain matrices; use it with some care. If a matrix contains small elements that are due to roundoff error [our case], balancing might scale them up to make them as significant as the other elements of the original matrix.").
I see here http://www.netlib.org/lapack/Errata/index2.html that xGEBAL still has some problems in LAPACK 3.4.2.
If you remove the balancing stage, expm is fine. Possible solutions:

1) remove balancing from expm (anyone has a counterexample?)
2) try to fix balance: if, after balancing, the matrix is more unbalanced than the original (for whatever reason), give back the original matrix and a warning.


Marco Caliari <caliari>
Project Member
Tue 27 Aug 2013 01:57:48 PM UTC, original submission:


I have ecountered some problem using expm function. The results I obtain have low accuracy.

To reproduce the the problem, I have attached two files:
-expm_data is a data file containing the matrices
tmpa : test matrix
expm_tmpa: result of expm(tmpa) in my PC
expa : expected result (checked both with matlab and an external fortran library)
-expm_test.m is a small script file loading the data and calculating the norm between the expected and the actual result.

In my case the norm has value 0.0953.

best regards



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

Attach Files:

Attached Files
file #28941:  expm_data added by None (3KiB - application/octet-stream)
file #28942:  expm_test.m added by None (282B - text/x-objcsrc)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by hardy (Posted a comment)
  • -email is unavailable- added by mtmiller (Updated the item)
  • -email is unavailable- added by mudit_sharma (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)

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

    Only project members can vote.




    Follow 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2018-03-28 mtmiller StatusConfirmed => Fixed
        Open/ClosedOpen => Closed
    2016-11-16 mtmiller Release3.6.4 => dev
    2013-08-29 caliari StatusNone => Confirmed
    2013-08-27 None Attached File- => Added expm_data, #28941
        Attached File- => Added expm_test.m, #28942

    Back to the top

    Powered by Savane 3.3