RDF/CDF eigenvalues: symmetric matrices, multiplicities

[rbeezer]

As the summary says, this patch improves eigenvalues of matrices with double-precision floating-point entries.

1. Uses SciPy's eigh routine for symmetric and Hermitian matrices.
2. Allows grouping of "equal" eigenvalues, based on a tolerance parameter.

Apply:

1. attachment: trac_11608-eigenvalues-symmetric-multiplicity-rebase.patch
2. attachment: trac_11608-eigenvalues-symmetric-multiplicity-update-v2.patch

CC: @jasongrout

Component: linear algebra

Keywords: sd40.5

Author: Rob Beezer

Reviewer: Volker Braun

Merged: sage-5.1.beta4

Issue created by migration from https://trac.sagemath.org/ticket/11608

[rbeezer]

Author: Rob Beezer

[loefflerd]

Work Issues: rebase

[loefflerd]

comment:3

Patch does not apply to 5.0.beta7 (see patchbot logs)

[rbeezer]

Attachment: trac_11608-eigenvalues-symmetric-multiplicity-rebase.patch.gz

Rebased on 5.0-beta10

[rbeezer]

Changed work issues from rebase to none

[dimpase]

comment:5

it's kind of weird that one coerces a complex matrix into real symmetric, but then treats it as a Hermitian one...
Surely with real symmetric matrices there might be better options, no?

[vbraun]

comment:6

Despite the name, scipy.linalg.eigh will call lapack syevr / heevr depending on whether the field is real / complex.

[rbeezer]

comment:7

The algorithm='symmetric' keyword is meant to allow the routine accomodate both base rings (RDF/CDF) while also allowing a user to specify that a matrix is known to be symmetric. In this case, the conversion from CDF to RDF is made by the routine and then the SciPy behavior that Volker notes will take place.

Any suggestions on making the docstring clearer? Maybe saying "then applies the algorithm for Hermitian matrices" should be stated differently?

[vbraun]

comment:8

On Fedora 16 x86_64, I get

sage -t  devel/sage-main/sage/matrix/matrix_double_dense.pyx
**********************************************************************
File "/home/vbraun/opt/sage-5.0.rc0/devel/sage-main/sage/matrix/matrix_double_dense.pyx", line 1570:
    sage: A.eigenvalues(algorithm='symmetric', tol=1.0e-20)
Expected:
    [(-2.0, 1), (-2.0, 2), (-2.0, 1), (1.0, 1), (1.0, 1), (1.0, 1), (1.0, 1), (1.0, 1), (3.0, 1)]
Got:
    [(-2.0, 1), (-2.0, 2), (-2.0, 1), (1.0, 1), (1.0, 1), (1.0, 2), (1.0, 1), (3.0, 1)]

ATLAS is from Fedora, CPU is a Sandy-Bridge i7 quad-core.

Minor nitpicks while we are at it:

• Can you return a tuple (=immutable) instead of a list?
• the symmetric docstring, how about: "converts self to a real matrix and applies the algorithm for Hermitian matrices". Note: technically you convert (explicit), not coerce (implicit).
• in the hermitian docstring, how about: "uses the :meth:~scipy.linalg.eigh method from SciPy, which applies only to real symmetric or complex Hermitian matrices. Since...". This would make it clear that the syevr implementation is used if possible.
• the ".. warning::" block is not indented, so it is not typeset correctly.

[vbraun]

Reviewer: Volker Braun

[rbeezer]

comment:9

Replying to @vbraun:

On Fedora 16 x86_64, I get...

Thanks, I'll increase the tolerance and that should work better across platforms. You'd think I'd get the hang of this numerical stuff...

Minor nitpicks while we are at it:

• Can you return a tuple (=immutable) instead of a list?

I'm not opposed, but other eigenvalue routines return mutable objects. I'd prefer to do one grand change across all matrix types, on a single-purpose ticket of its own.

sage: A = matrix(QQ, [[1]])
sage: ev=A.eigenvalues()   
sage: type(ev)
<class 'sage.structure.sequence.Sequence_generic'>
sage: ev[0]=5
sage: ev
[5]

• the symmetric docstring, how about: "converts self to a real matrix and applies the algorithm for Hermitian matrices". Note: technically you convert (explicit), not coerce (implicit).
• in the hermitian docstring, how about: "uses the :meth:~scipy.linalg.eigh method from SciPy, which applies only to real symmetric or complex Hermitian matrices. Since...". This would make it clear that the syevr implementation is used if possible.
• the ".. warning::" block is not indented, so it is not typeset correctly.

Thanks for the nits. I'll make some changes right now.

Rob

[rbeezer]

comment:10

"Update" patch:

1. Removes the 10^-20 tolerance doctest. Not really needed and will only cause trouble. (It was suppose to illustrate how a very small tolerance would split up "equal" eigenvalues.)
2. Adjusts documentation according to suggestions above.
3. Includes new pointers to routines to check symmetric and Hermitian properties.
4. No code changes.

Rob

[vbraun]

comment:11

Looks great!

[rbeezer]

comment:12

Replying to @vbraun:

Looks great!

Thanks, Volker!

This was built on various 5.0-beta's, so I just double-checked that it applies, tests and builds on 5.0-rc0.

Rob

[jdemeyer]

comment:14

The warning in the docstring is misformatted, see http://sagemath.org/doc/developer/conventions.html#documentation-strings

[rbeezer]

comment:15

Attachment: trac_11608-eigenvalues-symmetric-multiplicity-update-v2.patch.gz

I thought this was OK and I just triple-checked, and I think it is OK. One change: made "warning" all capital letters, this is the only change to the v2 version of the update patch.

The HTML documentation renders the warning blocks as they should. I'm going to move this back to positive review.

Jeroen - please let me know if I am really missing something here. Thanks.

Rob

[jdemeyer]

Changed keywords from none to sd40.5

[jdemeyer]

comment:16

Yes, it looks okay now. I don't remember what precisely was wrong before.

[jdemeyer]

Merged: sage-5.1.beta4