Opened 2 years ago

Closed 17 months ago

# the generalized eigenvalue problem over RDF/CDF

Reported by: Owned by: gh-mwageringel major sage-9.2 linear algebra scipy Markus Wageringel Sébastien Labbé N/A 254c6e9 254c6e932fa9f5b8f96dde02dda63936877157ea

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

### comment:1 Changed 2 years ago by gh-mwageringel

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

New commits:

 ​9a1161f `29243: implement generalized eigenvalues and eigenvectors over RDF/CDF`

### comment:2 follow-up: ↓ 6 Changed 2 years 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 2 years ago by gh-mwageringel (previous) (diff)

### comment:3 follow-up: ↓ 5 Changed 2 years 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 2 years ago by git

• 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 gh-mwageringel

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: ↓ 8 Changed 2 years ago by chapoton

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 2 years ago by git

• 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 pyx-file`

### comment:8 in reply to: ↑ 6 Changed 2 years ago by gh-mwageringel

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 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: ↓ 13 Changed 18 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 18 months ago by slabbe

• Reviewers set to Sébastien Labbé

### comment:12 Changed 18 months ago by git

• 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 gh-mwageringel

• Status changed from needs_info to needs_review

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 18 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: ↓ 16 ↓ 20 Changed 18 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 18 months ago by mkoeppe

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 18 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: ↓ 19 Changed 18 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 18 months ago by slabbe

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: ↓ 22 Changed 18 months ago by slabbe

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 17 months ago by git

• 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 gh-mwageringel

Should we try to revive #16607?

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

Yes, I think so, too.

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 17 months ago by gh-mwageringel

• Status changed from needs_work to needs_review

### comment:24 Changed 17 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 17 months ago by slabbe

• Status changed from needs_review to needs_info

### comment:26 Changed 17 months ago by git

• 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 gh-mwageringel

• Status changed from needs_info to needs_review

Good catch. That was a mistake, indeed.

### comment:28 Changed 17 months ago by slabbe

• Status changed from needs_review to positive_review

Good to go.

### comment:29 Changed 17 months ago by gh-mwageringel

Thank you for the review.

### comment:30 Changed 17 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.