bugGNU Octave - Bugs: bug #57196, eigs should detect...

 
 

bug #57196: eigs should detect symmetric/hermitian inputs and return real values only

Submitter:  Rik <rik5>
Submitted:  Sat 09 Nov 2019 03:26:45 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * dev Operating System:  * Any
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 04 Dec 2019 08:34:06 AM UTC, comment #8: 

Writing zhaupd for ARPACK and rewrite the ARPACK interface in octave (or even switch to Krylov-Schur) would be an excellent GSOC.

Marco Caliari <caliari>
Group Member
Thu 14 Nov 2019 10:07:00 PM UTC, comment #7: 

Made a small mistake in the first changeset which I corrected here (https://hg.savannah.gnu.org/hgweb/octave/rev/e557774a2604).  The symmetry test for complex matrices must be ishermitian(), not issymmetric().

Rik <rik5>
Group administrator
Wed 13 Nov 2019 06:53:05 PM UTC, comment #6: 

@jwe: In the end, it's not much of an inconvenience.  We would still need to check the input for symmetry in order to determine whether to call ZHAUPD or ZNAUPD: we pay the cost of that calculation already.

@David: Thank you for the review.

I checked in the change here https://hg.savannah.gnu.org/hgweb/octave/rev/6604b41ca63f.

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Wed 13 Nov 2019 05:11:48 PM UTC, comment #5: 

Since they acknowledge that ignoring non-zero imaginary parts is appropriate in this case, it's unfortunate that they didn't just supply a ZHAUPD subroutine that calls ZNAUPD and then zeros out the imaginary parts.

John W. Eaton <jwe>
Group administrator
Wed 13 Nov 2019 04:10:08 PM UTC, comment #4: 


The case near 502 is a real matrix that is promoted to complex because you are asked for the eigenvalues close to SIGMA where SIGMA has an imaginary part. The choice of SIGMA, even if complex shouldn't impact the symmmetric nature of the matrix used, so yes the patch should apply at that point as well

David

David Bateman <dbateman>
Group Member
Tue 12 Nov 2019 04:20:27 PM UTC, comment #3: 

I agree, I think it makes the most sense to take care of this in the C++ code. 

Attached is a patch which appears to do the job.  I'm slightly uncertain if I should have applied the fix to the second case beginning at line 502.  If you could review that it would be helpful.

Marking as "Patch Submitted".


(file #47848)

Rik <rik5>
Group administrator
Tue 12 Nov 2019 10:51:48 AM UTC, comment #2: 

Page 51 of the ARPACK user guide discuss why this is the case

https://pdfs.semanticscholar.org/1ac9/6899e1c355d78b55ad6b1cc038493de4a906.pdf

Basically the hermitian and non hermitian complex versions of ARPACK are the same cost in terms of computation and so they only implemented one of them. The manual just suggests ignoring the non zero imaginary part for hermitian problems. So using "real" on the eigenvalues is in fact the "right" solution, but should probably be implemented in _eigs_ rather than eigs.

David

David Bateman <dbateman>
Group Member
Tue 12 Nov 2019 10:41:30 AM UTC, comment #1: 

Rik,

The basic problem is that ARPACK doesn't supply a ZSAUPD function (or is that ZHAUPD for hermitian) that is the complex version of the DSAUPD function for real symmetric problems. So ARPACK doesn't profit from the complex symmetry in the solution of this type of problem.

The real solution would be to implement the ZSAUPD function,  but that seems too ambitious. Otherwise in _eigs_ it should be possible to detect hermitian matrices and at least move the conversion of the eigenvalues to real closer to the source of the problem

David

David Bateman <dbateman>
Group Member
Sat 09 Nov 2019 03:26:45 PM UTC, original submission:  

A Hermitian matrix has only real eigenvalues.  However, due to various well-known issues with numerical computation such as round-off error and limited precision, Octave will return a complex result.

Sample code:


octave:50> x = ones (10,10);
octave:51> z = complex (x, x);
octave:52> A = [sparse(x), z; z', sparse(x)];
octave:53> eigs (A)
ans =

   2.4142e+01 - 7.8571e-16i
  -4.1421e+00 - 1.2042e-15i
   3.6792e-16 + 2.2611e-16i
  -5.6318e-31 + 6.3019e-31i
  -2.4149e-32 + 9.7647e-32i
   1.9359e-32 - 9.3915e-32i

octave:54> imag (ans)
ans =

  -7.8571e-16
  -1.2042e-15
   2.2611e-16
   6.3019e-31
   9.7647e-32
  -9.3915e-32


Matlab, however, returns real values.  Apparently there is a check for symmetry/hermitian properties of the input and this decides whether real values or complex values are returned.  As a test with the above matrix A, calling eigs returns all real values, but a small change that disrupts the symmettry causes complex values to be returned.


A(1,2) = eps;
eigs (A)


We could implement something similar in Octave, and it would also resolve bug #57185.




Rik <rik5>
Group administrator

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #47848:  bug57196.cset added by rik5 (2KiB - application/octet-stream)

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by rik5 (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 6 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-11-13 rik5 StatusConfirmed Fixed
        Open/ClosedOpen Closed
    2019-11-13 dbateman StatusPatch Submitted Confirmed
    2019-11-12 rik5 Attached File- Added bug57196.cset, #47848
        StatusConfirmed Patch Submitted
    2019-11-09 rik5 Dependencies- bugs #57185 is dependent

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code