Opened 2 years ago
Closed 17 months ago
#29243 closed enhancement (fixed)
the generalized eigenvalue problem over RDF/CDF
Reported by:  ghmwageringel  Owned by:  

Priority:  major  Milestone:  sage9.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: 
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 1e14 [0.07692307692307694 0.0] [ 0.0 +infinity] sage: λ = D[0, 0] sage: v = V[:, 0] sage: (A * v  B * v * λ).norm() < 1e14 True
This is implemented using scipy.linalg.eig.
The changes include:
 an optional argument
other
is added toeigenmatrix_right/left
inmatrix2.pyx
as well as related functions inmatrix_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 2 years ago by
 Branch set to u/ghmwageringel/29243
 Commit set to 9a1161fdad1a85e0035b70ec5aec218c816deeee
 Status changed from new to needs_review
comment:2 followup: ↓ 6 Changed 2 years ago by
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?
comment:3 followup: ↓ 5 Changed 2 years ago by
Not exactly related to the purpose of this ticket... In the documentation, the following looks ugly
sage: spectrum[3] # tol 1e13 (1.0000000000000406, [(1.0, 0.49999999999996264, 1.9999999999998617, 0.499999999999958)], 1)
Because of # abs tol 1e13
flag, it can safely be replaced with the much more readable
sage: spectrum[3] # abs tol 1e13 (1.0, [(1.0, 0.5, 2.0, 0.5)], 1)
comment:4 Changed 2 years ago by
 Commit changed from 9a1161fdad1a85e0035b70ec5aec218c816deeee to fab67c99e341b6cf27829d2b77e5e761c1fd6ab4
Branch pushed to git repo; I updated commit sha1. New commits:
fab67c9  29243: decrease accuracy in doctest output

comment:5 in reply to: ↑ 3 Changed 2 years ago by
Replying to vdelecroix:
Because of
# abs tol 1e13
flag, it can safely be replaced with the much more readablesage: spectrum[3] # abs tol 1e13 (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 ; followup: ↓ 8 Changed 2 years ago by
Replying to ghmwageringel:
The patchbot marks
from six import string_typesas 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 2 years ago by
 Commit changed from fab67c99e341b6cf27829d2b77e5e761c1fd6ab4 to a8dc4997962682ff164da76306010b280b941fc2
Branch pushed to git repo; I updated commit sha1. New commits:
a8dc499  29243: remove the use of six in pyxfile

comment:8 in reply to: ↑ 6 Changed 2 years ago by
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 22 months ago by
 Milestone changed from sage9.1 to sage9.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 followup: ↓ 13 Changed 18 months ago by
 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 18 months ago by
 Reviewers set to Sébastien Labbé
comment:12 Changed 18 months ago by
 Commit changed from a8dc4997962682ff164da76306010b280b941fc2 to b2f04334c89040c5fb80b2be2be366dedc497d2e
Branch pushed to git repo; I updated commit sha1. New commits:
b2f0433  29243: improve error handling

comment:13 in reply to: ↑ 10 Changed 18 months ago by
 Status changed from needs_info to needs_review
Replying to slabbe:
Why
self.eigenvectors_left()
can't handleself.eigenvectors_left(other)
ifother
isNone
?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 18 months ago by
 Status changed from needs_review to needs_work
Few more suggestions.
 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``)
 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.
 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 argumentkeyword=value
. It seems you made sure of being backward compatible. First, in that change ofmatrix2.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 followups: ↓ 16 ↓ 20 Changed 18 months ago by
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 keywordonly 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 18 months ago by
Replying to ghmwageringel:
In the long run, Sage should make use of keywordonly arguments, a Python 3 feature, to avoid this kind of problem. See also #16607.
I agree. I have already started to make use of keywordonly 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 18 months ago by
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 nonbackward 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 followup: ↓ 19 Changed 18 months ago by
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 performancecritical. Should we try to revive #16607?
comment:19 in reply to: ↑ 18 Changed 18 months ago by
comment:20 in reply to: ↑ 15 ; followup: ↓ 22 Changed 18 months ago by
Replying to ghmwageringel:
Related to that, I will also make the following change in
matrix2.pyx
in order to make sure thatextend
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 17 months ago by
 Commit changed from b2f04334c89040c5fb80b2be2be366dedc497d2e to 081643ceb77e33e83a340dfefeb61ae7e64f6bcd
Branch pushed to git repo; I updated commit sha1. New commits:
081643c  29243: fix details

comment:22 in reply to: ↑ 20 Changed 17 months ago by
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 keywordonly arguments.
comment:23 Changed 17 months ago by
 Status changed from needs_work to needs_review
comment:24 Changed 17 months ago by
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 17 months ago by
 Status changed from needs_review to needs_info
comment:26 Changed 17 months ago by
 Commit changed from 081643ceb77e33e83a340dfefeb61ae7e64f6bcd to 254c6e932fa9f5b8f96dde02dda63936877157ea
Branch pushed to git repo; I updated commit sha1. New commits:
254c6e9  29243: fix typo

comment:27 Changed 17 months ago by
 Status changed from needs_info to needs_review
Good catch. That was a mistake, indeed.
comment:29 Changed 17 months ago by
Thank you for the review.
comment:30 Changed 17 months ago by
 Branch changed from u/ghmwageringel/29243 to 254c6e932fa9f5b8f96dde02dda63936877157ea
 Resolution set to fixed
 Status changed from positive_review to closed
New commits:
29243: implement generalized eigenvalues and eigenvectors over RDF/CDF