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 FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: ConfirmedAssigned to: None
Originator Name: vathomasOriginator Email: -unavailable-
Open/Closed: OpenRelease: 3.6.4
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.


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

    Date Changed By Updated Field Previous Value => Replaced By
    Thu 29 Aug 2013 12:08:18 PM UTCcaliariStatusNone=>Confirmed
    Tue 27 Aug 2013 01:57:48 PM UTCNoneAttached File-=>Added expm_data, #28941
      Attached File-=>Added expm_test.m, #28942

    Back to the top

    Powered by Savane 3.1-cleanup