bugGNU Octave - Bugs: bug #35779, Incorrect results of logm

 
 

bug #35779: Incorrect results of logm

Submitted by:  None
Submitted on:  Sat 10 Mar 2012 01:38:15 AM UTC  
 
Category: NoneSeverity: 3 - Normal
Priority: 7 - HighItem Group: Incorrect Result
Status: FixedAssigned to: Jordi Gutiérrez Hermoso <jordigh>
Originator Name: Masaru MuranishiOriginator Email: -unavailable-
Open/Closed: ClosedRelease: 3.6.1
Operating System: Microsoft Windows

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 Jump to the original submission)

Sun 18 Mar 2012 04:57:14 PM UTC, comment #11:

Done. Grafted on to the stable branch at 14481.

Rik <rik5>
Project Administrator
Sat 17 Mar 2012 06:34:20 PM UTC, comment #10:

Rik,

Perhaps it's good to transplant this change to stable. I didn't think of it when I applied it. Do you want to do it?

If you're using hg 2.0 or later, you might like this command:

http://www.selenic.com/mercurial/hg.1.html#graft

Otherwise, I can do it.

Jordi Gutiérrez Hermoso <jordigh>
Project AdministratorIn charge of this item.
Tue 13 Mar 2012 10:19:40 PM UTC, comment #9:

A small question, should this logm patch have been committed on stable? The behavior is a regression from 3.2.4, which worked, to 3.4.0 and beyond where it fails. It also seems pretty serious in that users will be generating incorrect results without any indication that the results of their analysis are flawed.

Worst case, someone designs a bridge using Octave and the results say it will work. When the bridge is built it collapses because the underlying algorithm was fine but Octave couldn't calculate things correctly.

Not calculating things correctly is a real black eye for anyone. Remember Intel and the FDIV bug on the Pentium (http://en.wikipedia.org/wiki/Pentium_FDIV_bug)?

Rik <rik5>
Project Administrator
Tue 13 Mar 2012 05:13:48 PM UTC, comment #8:

Some other tips for commit messages:

Please write

I.e., the description of each change should start on a separate line.

For style points, please also write "New function" instead of "Added function" or "Return retval" instead of "Changed to return retval".

And please, never "Fixed bug" or similar. That doesn't add any specific information about what was changed.

Finally, the changelog info in the commit message should say what changed, not why. If you need to explain why, that info should probably go in comments in the code.

John W. Eaton <jwe>
Project Administrator
Tue 13 Mar 2012 01:32:40 PM UTC, comment #7:

Great, I've pushed it:

http://hg.savannah.gnu.org/hgweb/octave/rev/6c3441f3146b

Closing this report.

Jordi Gutiérrez Hermoso <jordigh>
Project AdministratorIn charge of this item.
Tue 13 Mar 2012 01:03:24 PM UTC, comment #6:

Dear Jordi,

your changeset is fine for me. I'll try to make a changeset for the spline bug #35739.

Marco

Marco Caliari <caliari>
Project Member
Tue 13 Mar 2012 12:51:48 PM UTC, comment #5:

Thanks Marco. This is much more convenient for me to analyse and apply.

Attached is a modified version of your patch. I made the following changes:

1) The commit message has a one-line description including the number of the bug you're squashing.

2) Individual files you touched get their changes separately enumerated. If there a particular C or C++ function you changed, the general format is

  • file.cc (class1::function1): Made change foo. (class2::function2): Made change bar.

I omitted the parenthetical function names here because your changes aren't inside any C or C++ function.

3) Don't use norm or max(abs(...)) or similar for the assert function The assert function already has a 3-argument calling form for testing both absolute and relative error, and this form produces a more informative diagnostic when it fails.

I hope it's clear that these are improvements. The format of the commit message is important. In source tarballs, we generate changelogs from the commit messages using the "hg log --style changelog" command, and so the style of the commit messages has to match accordingly.

If you agree that these changes are correct, I will push my version of the patch. When you pull from Savannah, it will appear to you that you have a duplicate patch in your local repository, because you will have both your patch and my version.
I will explain later how to handle this situation.

(file #25356)

Jordi Gutiérrez Hermoso <jordigh>
Project AdministratorIn charge of this item.
Tue 13 Mar 2012 12:12:22 PM UTC, comment #4:

My first hg changeset!

(file #25354)

Marco Caliari <caliari>
Project Member
Mon 12 Mar 2012 01:10:15 PM UTC, comment #3:

Hi all.

For the case [0,-i;i,0] the fix is easy: just replace
with
in logm.m.
For the case [0,-1;1,0] the problem is in zrsf2csf.f (and in crsf2csf.f): around line 48 there is a off-by-one error. You should replace
with

Now I'm trying to build from hg sources in order to make an hg patch.

Marco

Marco Caliari <caliari>
Project Member
Sat 10 Mar 2012 05:41:42 PM UTC, comment #2:

I can confirm this bug on the most recent tip. I'm increasing the priority to high.

Thanks Marco for looking into this as the details are a bit beyond me. I have done some narrowing of the problem which I will share.

The bug was introduced between releases 3.2.4 and 3.4.0 when Octave switched from a logm() based on eig() to one based on schur() and a Pade approximation.

The bug is not, at first, obvious because it depends on the input matrix. The internal test cases used to validate the implementation in Octave did not happen to be of the correct type to trigger the bug.

For example, the original reporter can verify that logm works correctly on a great number of matrices.

This is the original matrix within round-off error.

The problem occurs when the matrix has complex eigenvalues. The initial Schur decomposition produces a quasi-upper triangular matrix which must be converted to a true upper triangular matrix by converting the real Schur form to the complex Schur form. The code that does this is written in Fortran and resides in libcruft/lapack-xtra/zrsf2csf.f.

The original reporter's matrix was A = [0 -1; -1 0]. Running the following code I can see that there is a problem with rsf2csf.

Rik <rik5>
Project Administrator
Sat 10 Mar 2012 11:05:29 AM UTC, comment #1:

Hi.

The first problem is easy to fix: it is due to the wrong assumption that if the diagonal of the Schur factorization is real, then the result has to be real. Removal of spurious imaginary parts has to be done in the correct way.

The second is actually a problem of rsf2csf.

is simply wrong because the result is not unitary. I will investigate in the next few days.

Marco

Marco Caliari <caliari>
Project Member
Sat 10 Mar 2012 01:38:15 AM UTC, original submission:

I'd like to report following.

In my understanding, logm is inverse function of expm, namely, following should be satisfied for any square matrix M.

M=logm(expm(M))

I tried following calculation using octave3.6.1.

A=[0 -i;i 0]
B=expm(A)
C=logm(B)

We expect A=C.
However, I got following results.

A=
0 + 0i -0 - 1i
0 + 1i 0 + 0i

B=
1.54308 + 0.00000i 0.00000 - 1.17520i
0.00000 + 1.17520i 1.54308 + 0.00000i

C=
-7.8496e-017 0.00000e+000
0.0000e+000 1.09991e-015

Obviously, C is different from A.

Moreover, I tried following.

A=[0 -1;1 0]
B=expm(A)
C=logm(B)

I got following.

A=
0 -1
1 0

B=
0.54030 -0.84147
0.84147 0.54030

C=
-5.4911e-016 1.4142e+000
-1.4142e+000 -1.8647e-016

Obviously, C is different from A.

I think that the matrix [0 -i;i 0] or [0 -1;1 0] are relatively simple matrix.
If octave can not give correct answer to this kind of relatively simple matrix, I concern this fact may spoile trust of octave.

For your information, I'd like to add results of same kind of calculation by scilab.

A=
0 - 1i
i 0

B= (=expm(A))
1.5430806 -1.1752012i
1.1752012i 1.5430806

C= (=logm(B))
0 -i
i -3.331E-16

A=
0 -1
1 0

B= (=expm(A))
0.5403023 -0.8414710
0.8414710 0.5403023

C= (=logm(B))
5.370D-14-2.220D-16i -1.+9.192D-17
1.-1.579D-16i 2.516D-16+2.22D-16i

I think scilab give correct answer to same calculation.

Anonymous

 

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

Attach File(s):
   
   
Comment:
   

Attached Files
file #25356:  logm-2.diff added by jordigh (3KiB - text/x-diff)
file #25354:  logm.cset added by caliari (4KiB - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -unavailable- added by jwe (Posted a comment)
  • -unavailable- added by jordigh (Updated the item)
  • -unavailable- added by rik5 (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 8 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Tue 13 Mar 2012 01:32:40 PM UTCjordighStatusPatch Submitted=>Fixed
      Open/ClosedOpen=>Closed
    Tue 13 Mar 2012 12:51:48 PM UTCjordighAttached File-=>Added logm-2.diff, #25356
      StatusConfirmed=>Patch Submitted
      Assigned toNone=>jordigh
    Tue 13 Mar 2012 12:12:22 PM UTCcaliariAttached File-=>Added logm.cset, #25354
    Sat 10 Mar 2012 05:41:42 PM UTCrik5Priority5 - Normal=>7 - High
      StatusNone=>Confirmed

    Back to the top


    Powered by Savane 3.1-cleanup1