bugGNU Octave - Bugs: bug #38577, Eig returns non-unitary...

 
 

bug #38577: Eig returns non-unitary transformation matrix

Submitter:  None
Submitted:  Sat 23 Mar 2013 04:44:18 PM UTC
   
 
Category:  None Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Incorrect Result
Status:  Invalid / Not an Octave Bug Assigned to:  None
Originator Name:  Peter Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * 3.6.2
Operating System:  * GNU/Linux Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 28 Mar 2013 08:49:33 AM UTC, comment #10: 

Hi,

indeed, I have LaPack 3.4.1 installed.

Using LD_LIBRARY_PATH to make octave using
ACML-5.3 I get correct results.

Amazing.

Thanks for the clarification,
Peter

Anonymous
Thu 28 Mar 2013 03:42:45 AM UTC, comment #9: 

I also see this bug with LAPACK 3.4.1 but not with LAPACK 3.4.2 from Debian experimental. Looks resolved.

Mike Miller <mtmiller>
Group Member
Wed 27 Mar 2013 01:52:18 PM UTC, comment #8: 

I managed to get access to my development machine, and I confirm that Octave is calling LAPACK's dsyev correctly. The problem is likely therefore to be in LAPACK:

http://www.netlib.org/lapack/bug_list.html

This seems like bug 0098, or at least related to it.

I am running LAPACK 3.4.1, where this bug is present. What LAPACK version are you using? I am tentatively closing this report. If you still believe that the bug is in Octave and not in LAPACK, we may reopen it.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Wed 27 Mar 2013 01:19:10 PM UTC, comment #7: 


> Hermitian matrices can never be spectrally defective.


But they can be very close to defective. A defective matrix is somewhat numerically invisible, since any perturbation destroys its defect. This is why the Jordan form is somewhat useless, since it's numerically unstable.

> Eig returning a zero eigenvector for a hermitian matrix is a sign of a bug.


Not necessarily, a defective matrix can be close to being Hermitian, e.g. [1 eps; 0 1].

The eigenvalue routine that Octave uses isn't implemented in Octave. It's just calling LAPACK's dgeev family of functions:

http://hg.savannah.gnu.org/hgweb/octave/file/a2f65b8f1955/liboctave/numeric/EIG.cc#l31

I'm not sure how your matrix is getting close to being defective, and it seems to depend on the algorithm since the Schur factorisation doesn't exhibit this problem. I have had reports that different Octave builds on Windows don't have this problem.

I don't have a development machine right now to test on, but if you're able to build Octave and step through it with a debugger, you should check if indeed the LAPACK dseev symmetric eigenvalue algorithm is being called for your matrix. If it isn't, then there's a possible bug in Octave. If it is, then the problem is in LAPACK.

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Wed 27 Mar 2013 09:08:59 AM UTC, comment #6: 

Hermitian matrices can never be spectrally defective. They are guaranteed to have an orthogonal basis of eigenvectors. Eig returning a zero eigenvector for a hermitian matrix is a sign of a bug.

For hermitian matrices eig and schur decompositions coincide. The schur decomposition
S = U' A U
implies that the triangular matrix S is hermitian if A is. Being both triangular and hermitian forces it to be diagonal.

Anonymous
Mon 25 Mar 2013 04:06:27 PM UTC, comment #5: 

I can confirm this. My first instinct is that the matrix is numerically defective, so that v not being of full rank would make sense:


a = sparse ([200,100], [100, 200], [-1,-1]);
[v,~] = eig (a);
rank (v)


That's of rank 199, not 200. It is akin to the following situation:


b = [0 1; 0 0];
[v,~] = eig (b);
rank (v)


However, this doesn't seem to explain it. The original 200x200 matrix is not apparently numerically defective:


rank (null (a)) ## Since 0 is the repeated eigenvalue...


Compare with


rank (null (b))


which is 1, showing the spectral defect.

At any rate, eig returning a zero eigenvector is a sign of a defective matrix. I would need to study deeper the Schur decomposition method to understand this... Curiously enough, as a workaround, you can use schur() instead of eig, since in this case it does produce the eigenvectors (the upper triangular S matrix ends up being not just upper triangular, but diagonal).

I have heard reports that this problem doesn't exist in some Windows builds. Which version of LAPACK and which BLAS are you using?

Jordi Gutiérrez Hermoso <jordigh>
Group Member
Mon 25 Mar 2013 09:24:46 AM UTC, comment #4: 

Hi,

< http://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors > :
An eigenvector of a square matrix A is a non-zero _vecto_  v that, when multiplied by A, yields the original vector multiplied by a single number .

<http://mathworld.wolfram.com/Eigenvector.html > :
Eigenvectors may not be equal to the zero vector.

Peter

Anonymous
Mon 25 Mar 2013 07:26:44 AM UTC, comment #3: 

Hi,


> The zero vector for x solves this equation and is orthogonal
> to the other eigenvectors which are all zeros except for a single
> '1' on the diagonal. Octave appears to place this as the very
> last eigenvector >every time. You can check this with
> [...]
> So is this issue merely definitional?


Actually it is not merely definitional.
A * 0 = 0 is true for any matrix. The main point is that the
eigen vectors of a hermitian matrix form a complete basis.
With one zero vector this is not true, one vector is missing.

Actually, I didn't construct this example on purpose, but found it
while debugging a  script, where I perform base transformations into the basis of eigen vectors of a hermitian matrix. At least
in quantum mechanics this is essential.

Peter

P.S. I'm sorry, but I still have to read the bug-tracker documentation to understand the mark-up.


Anonymous
Sun 24 Mar 2013 06:54:52 PM UTC, comment #2: 

Apparently the markup language ate some of my last post.  Here it is again:

The eigenvalues in E are correct ([1, -1, lots of zeros]). This can be checked with sum (E(:)) which is always 0.

For the eigenvectors with zero eigenvalues one needs to solve


A*x = 0*x


The zero vector for x solves this equation and is orthogonal to the other eigenvectors which are all zeros except for a single '1' on the diagonal.  Octave appears to place this as the very last eigenvector every time. You can check this with



diag (U*U')(end-1).


So is this issue merely definitional?  Octave believes the zero vector is a valid eigenvector because it meets the criteria above, but ordinary usage does not include the zero vector as an eigenvector?

Rik <rik5>
Group administrator
Sun 24 Mar 2013 03:38:41 PM UTC, comment #1: 

The eigenvalues in E are correct ([1, -1, lots of zeros]).  This can be checked with sum (E(:)) which is always 0.

For the eigenvectors with zero eigenvalues one needs to solve Ax = 0*x.  The zero vector for x solves this equation and is orthogonal to the other eigenvectors which are all zeros except for a single '1' on the diagonal.  Octave appears to place this as the very last eigenvector every time.  You can check this with diag (U*U')(end-1).

So in this issue merely definitional?  Octave believes the zero vector is a valid eigenvector because it meets the criteria above, but ordinary usage does not include the zero vector as an eigenvector?


Rik <rik5>
Group administrator
Sat 23 Mar 2013 04:44:18 PM UTC, original submission:  

Dear All,

the following code


M=200;
A=zeros(M,M);

Index1= ( M - mod(M,2) ) / 2;

A(M,M-Index1) = -1.0;
A(M-Index1,M) = -1.0;

[ U, E ] = eig( A );

z = max ( abs( diag( U * U' ) - ones (M,1)))


produces a '1' as result, instead of  a numerical zero.
Tested on octave 3.2.4 and 3.6.2.
It looks like one eigen vector is just a zero vector.
On 3.6.2 it work for M<130, and fails for M >=130.
At least I checked until M=1000.

Peter

Anonymous

 

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

Attach Files:
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mtmiller (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  •  

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

    Date Changed by Updated Field Previous Value => Replaced by
    2013-03-27 jordigh StatusConfirmed Invalid / Not an Octave Bug
        Open/ClosedOpen Closed
    2013-03-25 jordigh StatusNone Confirmed

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code