the generalized eigenvalue problem over RDF/CDF

[mwageringel]

This ticket adds support for solving the generalized eigenvalue problem A v = λ B v for two matrices A,B over RDF or CDF.

Example:

sage: A = matrix.identity(RDF, 2)
sage: B = matrix(RDF, [[3, 5], [6, 10]])
sage: D, V = A.eigenmatrix_right(B); D   # tol 1e-14
[0.07692307692307694                 0.0]
[                0.0           +infinity]

sage: λ = D[0, 0]
sage: v = V[:, 0]
sage: (A * v  - B * v * λ).norm() < 1e-14
True

This is implemented using scipy.linalg.eig.

The changes include:

• an optional argument other is added to eigenmatrix_right/left in matrix2.pyx as well as related functions in matrix_double_dense.pyx
• an optional keyword homogeneous is added to obtain the generalized eigenvalues in terms of homogeneous coordinates
• improvements to the documentation

Component: linear algebra

Keywords: scipy

Author: Markus Wageringel

Branch/Commit: 254c6e9

Reviewer: Sébastien Labbé

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

[mwageringel]

Author: Markus Wageringel

[mwageringel]

Branch: u/gh-mwageringel/29243

[mwageringel]

New commits:

9a1161f

29243: implement generalized eigenvalues and eigenvectors over RDF/CDF

[mwageringel]

Commit: 9a1161f

[mwageringel]

comment:2

The patchbot marks

from six import string_types

as invalid in Cython code. Why is that, and what is a suitable replacement for checking whether a Python object is a string?

[videlec]

comment:3

Not exactly related to the purpose of this ticket... In the documentation, the following looks ugly

sage: spectrum[3]  # tol 1e-13
(-1.0000000000000406, [(1.0, -0.49999999999996264, 1.9999999999998617, 0.499999999999958)], 1)

Because of # abs tol 1e-13 flag, it can safely be replaced with the much more readable

sage: spectrum[3]  # abs tol 1e-13
(-1.0, [(1.0, -0.5, 2.0, 0.5)], 1)

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

fab67c9

29243: decrease accuracy in doctest output

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 9a1161f to fab67c9

[mwageringel]

comment:5

Replying to @videlec:

Because of # abs tol 1e-13 flag, it can safely be replaced with the much more readable

sage: spectrum[3]  # abs tol 1e-13
(-1.0, [(1.0, -0.5, 2.0, 0.5)], 1)

I was worried this might be a bit misleading because the actual result will never look like that. However, it makes the documentation much more readable indeed, so I have changed it as you suggested.

That reminds me though: it would be nice if Sage supported the IPython %precision directive in order to get rid of numerical noise like that from the printed IPython output.

[fchapoton]

comment:6

Replying to @mwageringel:

The patchbot marks

from six import string_types

as invalid in Cython code. Why is that, and what is a suitable replacement for checking whether a Python object is a string?

NEVER import six in a cython file. You can use "str".

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

a8dc499

29243: remove the use of six in pyx-file

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from fab67c9 to a8dc499

[mwageringel]

comment:8

Replying to @fchapoton:

NEVER import six in a cython file. You can use "str".

Thanks for the clarification. It is fixed now.

[mkoeppe]

comment:9

Batch modifying tickets that will likely not be ready for 9.1, based on a review of the ticket title, branch/review status, and last modification date.

[seblabbe]

comment:10

The code looks okay. The patchbot says green light.

Personnaly, I do not like so much this part of the changes:

-        evecs = self.eigenvectors_left()
+        if other is None:
+            evecs = self.eigenvectors_left()
+        else:
+            try:
+                evecs = self.eigenvectors_left(other=other)
+            except TypeError as e:
+                raise NotImplementedError('generalized eigenvector '
+                                          'decomposition is implemented '
+                                          'for RDF and CDF, but not for %s'
+                                          % self.base_ring()) from e

Why self.eigenvectors_left() can't handle self.eigenvectors_left(other) if other is None?

Why should a TypeError be raised if the base ring is not RDF or CDF?

I think the raise NotImplementedError should not be raised there, but deeper in the code. No?

[seblabbe]

Reviewer: Sébastien Labbé

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from a8dc499 to b2f0433

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

b2f0433

29243: improve error handling

[mwageringel]

comment:13

Replying to @seblabbe:

Why self.eigenvectors_left() can't handle self.eigenvectors_left(other) if other is None?

Why should a TypeError be raised if the base ring is not RDF or CDF?

I think the raise NotImplementedError should not be raised there, but deeper in the code. No?

Not every implementation of eigenvectors_left supported the argument other, in which case a TypeError would be raised. I do agree though that it is cleaner to add the additional argument to all implementations of eigenvectors_left and handle it there.

I have removed that section of the code and updated the branch accordingly. Thank you for the suggestion.

[seblabbe]

comment:14

Few more suggestions.

1. I would document the input other everywhere, something like:

-        - ``other`` -- not supported
+        - ``other`` -- a square matrix `B` (default: ``None``) in a generalized
+          eigenvalue problem; if ``None``, an ordinary eigenvalue problem is
+          solved (currently supported only if the base ring of ``self`` is ``RDF`` 
+          or ``CDF``)

2. In the following change:

-        if not self.is_square():
+        if not self.is_square() or other is not None and not other.is_square():
             msg = 'matrix must be square, not {0} x {1}'
+            m = self if not self.is_square() else other
+            raise ValueError(msg.format(m.nrows(), m.ncols()))

I would keep one if for self.is_square() and write another if below for other is not None and other.is_square(). Merging the two needs to run the code self.is_square() twice, which I find not clean.

3. I am worried about the changes not being backward compatible, mainly because of the changes in the ordering of the input arguments. eigenvalues methods are used a lot everywhere in the code of everyone. There are surely some places where people did not write the keyword argument keyword=value. It seems you made sure of being backward compatible. First, in that change of matrix2.pyx:

-    def eigenvectors_left(self,extend=True):
+    def eigenvectors_left(self, other=None, extend=True):

you check if other is a boolean type and you send a deprecation warning:

+        if other is not None:
+            if isinstance(other, bool):
+                # for backward compatibility
+                from sage.misc.superseded import deprecation
+                deprecation(29243,
+                            '"extend" should be used as keyword argument')
+                extend = other

But in matrix_double_dense.pyx,

-    def eigenvalues(self, algorithm='default', tol=None):
+    def eigenvalues(self, other=None, algorithm='default', tol=None,
+                    homogeneous=False):

you solve the problem like this:

+        if isinstance(other, str):
+            # for backward compatibilty, allow algorithm to be passed as first
+            # positional argument
+            algorithm = other
+            other = None

with no deprecation warning. Can you tell why?

[mwageringel]

comment:15

I will make the changes for 1 and 2.

Regarding the deprecation in 3, I am following what Travis has recently suggested in #29543 comment:2. In matrix_double_dense.pyx I have not used a deprecation warning just because the code was written before that. I will change it to issue a deprecation warning as well, as this additional argument processing is undesirable in general, and will make sure the function is written in a backward compatible way.

In the long run, Sage should make use of keyword-only arguments, a Python 3 feature, to avoid this kind of problem. See also #16607.

Related to that, I will also make the following change in matrix2.pyx in order to make sure that extend is only used as a keyword, from now on. The code is already written in a way that makes this backward compatible.

-    def eigenvectors_left(self, other=None, extend=True):
+    def eigenvectors_left(self, other=None, *, extend=True):

This does not work for the methods with more than one optional keyword argument, though – I hope that #16607 will find a more general solution for this eventually.

[mkoeppe]

comment:16

Replying to @mwageringel:

In the long run, Sage should make use of keyword-only arguments, a Python 3 feature, to avoid this kind of problem. See also #16607.

I agree. I have already started to make use of keyword-only arguments for new keywords in some tickets. Deprecating positional use of existing keywords arguments, on the other hand, is another story -- and the performance impact of trying to do that is what stopped #16607 if I remember correctly.

[mwageringel]

comment:17

It would be good to recommend this for new code, for example in the developer guide.

As for existing code, performance can be a problem, particularly if we wanted to change it across the whole library. In that case, a non-backward compatible release might be a better option.

However, I am thinking mostly of a decorator that helps in the few cases where we actually want to add a new positional argument, such as:

-    def eigenvalues(self, algorithm='default', tol=None):
+    def eigenvalues(self, other=None, *, algorithm='default', tol=None, homogeneous=False):

Manually writing this in a backward compatible way quickly gets out of hand if there are several optional keywords. Using a decorator here for one year of deprecation would not harm performance too much.

[mkoeppe]

comment:18

I agree that preparing a decorator for this type of API change with deprecation is the way to go if we are careful about using it in places that are not performance-critical. Should we try to revive #16607?

[seblabbe]

comment:19

Replying to @mkoeppe:

Should we try to revive #16607?

+1, I think there will be a need for that in the future.

[seblabbe]

comment:20

Replying to @mwageringel:

Related to that, I will also make the following change in matrix2.pyx in order to make sure that extend is only used as a keyword, from now on. The code is already written in a way that makes this backward compatible.

-    def eigenvectors_left(self, other=None, extend=True):
+    def eigenvectors_left(self, other=None, *, extend=True):

+1

Looks good. I did not know that feature.

I wait for your changes before looking at the branch again.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from b2f0433 to 081643c

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

081643c

29243: fix details

[mwageringel]

comment:22

Replying to @seblabbe:

Replying to @mkoeppe:

Should we try to revive #16607?

+1, I think there will be a need for that in the future.

Yes, I think so, too.

Replying to @seblabbe:

I wait for your changes before looking at the branch again.

I have updated the branch. The changes should now be completely backward compatible, raising deprecation warnings for positional use of keyword-only arguments.

[seblabbe]

comment:24

I have a question. In the following changes inside of the function def eigenvalues(self, other=None, algorithm='default', tol=None, *, homogeneous=False): of the file src/sage/matrix/matrix_double_dense.pyx, you make the following change:

+            deprecation(29243, '"extend" and "tol" should be used as '
+                               'keyword argument only')

Is it an error that the deprecation warning mentions "extend"? It is not an argument of that method...

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Branch pushed to git repo; I updated commit sha1. New commits:

254c6e9

29243: fix typo

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 081643c to 254c6e9

[mwageringel]

comment:27

Good catch. That was a mistake, indeed.

[seblabbe]

comment:28

Good to go.

[mwageringel]

comment:29

Thank you for the review.

[vbraun]

Changed branch from u/gh-mwageringel/29243 to 254c6e9