[Octave-bug-tracker] [bug #35779] Incorrect results of logm

[Rik]

Update of bug #35779 (project octave):

                Priority:              5 - Normal => 7 - High               
                  Status:                    None => Confirmed              

    _______________________________________________________

Follow-up 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 ->




    _______________________________________________________

Reply to this item at:

  <http://savannah.gnu.org/bugs/?35779>

_______________________________________________
  Message sent via/by Savannah
  http://savannah.gnu.org/