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


bug #39879: expm returns inaccurate results

Submitted by:  None
Submitted on:  Tue Aug 27 13:57:48 2013  
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: ConfirmedAssigned to: None
Originator Name: vathomasOriginator Email: -unavailable-
Open/Closed: OpenRelease: dev
Operating System: GNU/Linux

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)

Tue Nov 22 21:57:28 2016, 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 Dec 15 18:07:15 2015, 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 Jun 5 11:59:24 2014, 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 Jun 5 10:31:28 2014, 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 Jun 5 10:00:03 2014, 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 Jun 5 08:40:02 2014, 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 Jun 5 07:47:46 2014, 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 Jun 5 07:34:43 2014, comment #2:

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


Thu Aug 29 12:08:18 2013, 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 Aug 27 13:57:48 2013, 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 File(s):

Attached Files
file #28941:  expm_data added by None (3kB - 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
  • -unavailable- added by hardy (Posted a comment)
  • -unavailable- added by mtmiller (Updated the item)
  • -unavailable- added by mudit_sharma (Posted a comment)
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by None (Submitted the item)

    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 4 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Wed Nov 16 22:56:06 2016mtmillerRelease3.6.4=>dev
    Thu Aug 29 12:08:18 2013caliariStatusNone=>Confirmed
    Tue Aug 27 13:57:48 2013NoneAttached File-=>Added expm_data, #28941
      Attached File-=>Added expm_test.m, #28942

    Back to the top

    Powered by Savane 3.1-cleanup