#29243 closed enhancement (fixed)

the generalized eigenvalue problem over RDF/CDF

Reported by: gh-mwageringel Owned by:
Priority: major Milestone: sage-9.2
Component: linear algebra Keywords: scipy
Cc: Merged in:
Authors: Markus Wageringel Reviewers: Sébastien Labbé
Report Upstream: N/A Work issues:
Branch: 254c6e9 (Commits, GitHub, GitLab) Commit: 254c6e932fa9f5b8f96dde02dda63936877157ea
Dependencies: Stopgaps:

Status badges

Description

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

Change History (30)

comment:1 Changed 19 months ago by gh-mwageringel

  • Authors set to Markus Wageringel
  • Branch set to u/gh-mwageringel/29243
  • Commit set to 9a1161fdad1a85e0035b70ec5aec218c816deeee
  • Status changed from new to needs_review

New commits:

9a1161f29243: implement generalized eigenvalues and eigenvectors over RDF/CDF

comment:2 follow-up: Changed 19 months ago by gh-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?

Last edited 19 months ago by gh-mwageringel (previous) (diff)

comment:3 follow-up: Changed 19 months ago by vdelecroix

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)

comment:4 Changed 19 months ago by git

  • Commit changed from 9a1161fdad1a85e0035b70ec5aec218c816deeee to fab67c99e341b6cf27829d2b77e5e761c1fd6ab4

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

fab67c929243: decrease accuracy in doctest output

comment:5 in reply to: ↑ 3 Changed 19 months ago by gh-mwageringel

Replying to vdelecroix:

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.

comment:6 in reply to: ↑ 2 ; follow-up: Changed 19 months ago by chapoton

Replying to gh-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".

comment:7 Changed 19 months ago by git

  • Commit changed from fab67c99e341b6cf27829d2b77e5e761c1fd6ab4 to a8dc4997962682ff164da76306010b280b941fc2

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

a8dc49929243: remove the use of six in pyx-file

comment:8 in reply to: ↑ 6 Changed 19 months ago by gh-mwageringel

Replying to chapoton:

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

Thanks for the clarification. It is fixed now.

comment:9 Changed 18 months ago by mkoeppe

  • Milestone changed from sage-9.1 to sage-9.2

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.

comment:10 follow-up: Changed 13 months ago by slabbe

  • Status changed from needs_review to needs_info

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?

comment:11 Changed 13 months ago by slabbe

  • Reviewers set to Sébastien Labbé

comment:12 Changed 13 months ago by git

  • Commit changed from a8dc4997962682ff164da76306010b280b941fc2 to b2f04334c89040c5fb80b2be2be366dedc497d2e

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

b2f043329243: improve error handling

comment:13 in reply to: ↑ 10 Changed 13 months ago by gh-mwageringel

  • Status changed from needs_info to needs_review

Replying to slabbe:

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.

comment:14 Changed 13 months ago by slabbe

  • Status changed from needs_review to needs_work

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``)
  1. 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.

  1. 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?

comment:15 follow-ups: Changed 13 months ago by gh-mwageringel

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.

comment:16 in reply to: ↑ 15 Changed 13 months ago by mkoeppe

Replying to gh-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.

comment:17 Changed 13 months ago by gh-mwageringel

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.

comment:18 follow-up: Changed 13 months ago by mkoeppe

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?

comment:19 in reply to: ↑ 18 Changed 13 months ago by slabbe

Replying to mkoeppe:

Should we try to revive #16607?

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

comment:20 in reply to: ↑ 15 ; follow-up: Changed 13 months ago by slabbe

Replying to gh-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.

comment:21 Changed 13 months ago by git

  • Commit changed from b2f04334c89040c5fb80b2be2be366dedc497d2e to 081643ceb77e33e83a340dfefeb61ae7e64f6bcd

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

081643c29243: fix details

comment:22 in reply to: ↑ 20 Changed 13 months ago by gh-mwageringel

Replying to slabbe:

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 slabbe:

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.

comment:23 Changed 13 months ago by gh-mwageringel

  • Status changed from needs_work to needs_review

comment:24 Changed 13 months ago by slabbe

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...

comment:25 Changed 13 months ago by slabbe

  • Status changed from needs_review to needs_info

comment:26 Changed 13 months ago by git

  • Commit changed from 081643ceb77e33e83a340dfefeb61ae7e64f6bcd to 254c6e932fa9f5b8f96dde02dda63936877157ea

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

254c6e929243: fix typo

comment:27 Changed 13 months ago by gh-mwageringel

  • Status changed from needs_info to needs_review

Good catch. That was a mistake, indeed.

comment:28 Changed 13 months ago by slabbe

  • Status changed from needs_review to positive_review

Good to go.

comment:29 Changed 13 months ago by gh-mwageringel

Thank you for the review.

comment:30 Changed 13 months ago by vbraun

  • Branch changed from u/gh-mwageringel/29243 to 254c6e932fa9f5b8f96dde02dda63936877157ea
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.