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.

 

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

Thu 05 Jun 2014 11:59:24 AM UTC, comment #7:

OK,

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

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

Thomas

Anonymous
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

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

Hi,
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
3.9589e-015

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.

Best

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

Hi,

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?

Thanks,

Marco

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.

Best

Anonymous
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

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

Hi,

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

Anonymous

 

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

Attach File(s):
   
   
Comment:
   

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