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

 
 

bug #35779: Incorrect results of logm

Submitter:  None
Submitted:  Sat 10 Mar 2012 01:38:15 AM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  7 - High Item Group:  Incorrect Result
Status:  Fixed Assigned to:  jordigh
Originator Name:  Masaru Muranishi Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.6.1
Operating System:  * Microsoft Windows Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

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>
Group 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>
Group Member
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>
Group administrator
Tue 13 Mar 2012 05:13:48 PM UTC, comment #8: 

Some other tips for commit messages:

Please write


* file.cc (class::fcn): A change.
(class::other_fcn): Some other change.
* other_file.cc (fcn): Yet another change.


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>
Group 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>
Group Member
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>
Group 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>
Group Member
Tue 13 Mar 2012 12:12:22 PM UTC, comment #4: 

My first hg changeset!

(file #25354)

Marco Caliari <caliari>
Group 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

if (real_eig)
  s = real (s);
endif

with

if (real_eig && isreal(A))
  s = real (s);
endif

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

c(j+1) = 1

with

if (j+1 < n) then
  c(j+1) = 1
end if


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

Marco

Marco Caliari <caliari>
Group 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.


A = [1, 2; 3, 4];
expm (logm (A))
warning: logm: principal matrix logarithm is not defined for matrices with negative eigenvalues; computing non-principal logarithm
ans =

   1.0000 - 0.0000i   2.0000 + 0.0000i
   3.0000 + 0.0000i   4.0000 + 0.0000i


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.


A = [0 -1; 1 0];
eig (A)
ans =

   0 + 1i
   0 - 1i

[u, s] = schur (A)
u =

   1   0
   0   1

s =

   0  -1
   1   0

u * s * u'
ans =

   0  -1
   1   0

<- This is correct.  The matrix u is unitary and I can recover matrix A ->

[U, S] = rsf2csf (u, s)
U =

   0.70711 + 0.00000i   0.00000 + 1.00000i
   0.00000 + 1.00000i   0.70711 + 0.00000i

S =

   0.00000 + 1.00000i   0.00000 + 0.00000i
   0.00000 + 0.00000i   0.00000 - 1.00000i

<- This is only partially correct.  S is correct and the eigenvalues are +1i, -1i. ->

U * U'

ans =

   1.50000   0.00000
   0.00000   1.50000

<- This is incorrect.  U should be unitary and this result should be [1 0; 0 1] ->

U * S * U'
ans =

   0.00000 - 0.50000i   1.41421 - 0.00000i
  -1.41421 - 0.00000i   0.00000 + 0.50000i

<- And because U is incorrect it is not possible to recover the matrix A ->



Rik <rik5>
Group 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.

rsf2csf(eye(2),[0,-1;1,0])


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

Marco

Marco Caliari <caliari>
Group 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 Files:
   
   
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
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by jordigh (Updated the item)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 8 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2012-03-13 jordigh StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2012-03-13 jordigh Attached File- Added logm-2.diff, #25356
        StatusConfirmed Patch Submitted
        Assigned toNone jordigh
    2012-03-13 caliari Attached File- Added logm.cset, #25354
    2012-03-10 rik5 Priority5 - Normal 7 - High
        StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-4448.
    Corresponding source code