Trac #29243: the generalized eigenvalue problem over RDF/CDF

This ticket adds support for solving the [https://en.wikipedia.org/wiki/
Eigendecomposition_of_a_matrix#Generalized_eigenvalue_problem
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 [https://docs.scipy.org/doc/scipy/reference/ge
nerated/scipy.linalg.eig.html 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

URL: https://trac.sagemath.org/29243
Reported by: gh-mwageringel
Ticket author(s): Markus Wageringel
Reviewer(s): Sébastien Labbé

Author: Release Manager

src/sage/matrix/matrix2.pyx

@@ -6390,10 +6390,21 @@ cdef class Matrix(Matrix1):
         self.cache('eigenvalues', eigenvalues)
         return eigenvalues
 
-    def eigenvectors_left(self,extend=True):
+    def eigenvectors_left(self, other=None, *, extend=True):
         r"""
         Compute the left eigenvectors of a matrix.
 
+        INPUT:
+
+        - ``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``)
+
+        - ``extend`` -- boolean (default: ``True``)
+
+        OUTPUT:
+
         For each distinct eigenvalue, returns a list of the form (e,V,n)
         where e is the eigenvalue, V is a list of eigenvectors forming a
         basis for the corresponding left eigenspace, and n is the algebraic
@@ -6402,8 +6413,9 @@ cdef class Matrix(Matrix1):
         If the option extend is set to False, then only the eigenvalues that
         live in the base ring are considered.
 
-        EXAMPLES: We compute the left eigenvectors of a `3\times 3`
-        rational matrix.
+        EXAMPLES:
+
+        We compute the left eigenvectors of a `3\times 3` rational matrix.
 
         ::
 
@@ -6436,7 +6448,37 @@ cdef class Matrix(Matrix1):
             (0, 0, 1)
             ], 1)]
 
+        TESTS::
+
+            sage: A = matrix(QQ, [[1, 2], [3, 4]])
+            sage: B = matrix(QQ, [[1, 1], [0, 1]])
+            sage: A.eigenvectors_left(B)
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: generalized eigenvector decomposition is
+            implemented for RDF and CDF, but not for Rational Field
+
+        Check the deprecation::
+
+            sage: matrix(QQ, [[1, 2], [3, 4]]).eigenvectors_left(False)
+            doctest:...: DeprecationWarning: "extend" should be used as keyword argument
+            See https://trac.sagemath.org/29243 for details.
+            []
         """
+        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
+                other = None
+            else:
+                raise NotImplementedError('generalized eigenvector '
+                                          'decomposition is implemented '
+                                          'for RDF and CDF, but not for %s'
+                                          % self.base_ring())
+
         x = self.fetch('eigenvectors_left')
         if not x is None:
             return x
@@ -6476,10 +6518,21 @@ cdef class Matrix(Matrix1):
 
     left_eigenvectors = eigenvectors_left
 
-    def eigenvectors_right(self, extend=True):
+    def eigenvectors_right(self, other=None, *, extend=True):
         r"""
         Compute the right eigenvectors of a matrix.
 
+        INPUT:
+
+        - ``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``)
+
+        - ``extend`` -- boolean (default: ``True``)
+
+        OUTPUT:
+
         For each distinct eigenvalue, returns a list of the form (e,V,n)
         where e is the eigenvalue, V is a list of eigenvectors forming a
         basis for the corresponding right eigenspace, and n is the
@@ -6488,8 +6541,9 @@ cdef class Matrix(Matrix1):
         closure of the base field where this is implemented; otherwise
         it will restrict to eigenvalues in the base field.
 
-        EXAMPLES: We compute the right eigenvectors of a
-        `3\times 3` rational matrix.
+        EXAMPLES:
+
+        We compute the right eigenvectors of a `3\times 3` rational matrix.
 
         ::
 
@@ -6511,17 +6565,57 @@ cdef class Matrix(Matrix1):
             sage: delta = eval*evec - A*evec
             sage: abs(abs(delta)) < 1e-10
             True
+
+        TESTS::
+
+            sage: A = matrix(QQ, [[1, 2], [3, 4]])
+            sage: B = matrix(QQ, [[1, 1], [0, 1]])
+            sage: A.eigenvectors_right(B)
+            Traceback (most recent call last):
+            ...
+            NotImplementedError: generalized eigenvector decomposition is
+            implemented for RDF and CDF, but not for Rational Field
         """
-        return self.transpose().eigenvectors_left(extend=extend)
+        return self.transpose().eigenvectors_left(other=other, extend=extend)
 
     right_eigenvectors = eigenvectors_right
 
-    def eigenmatrix_left(self):
+    def eigenmatrix_left(self, other=None):
         r"""
-        Return matrices D and P, where D is a diagonal matrix of
-        eigenvalues and P is the corresponding matrix where the rows are
-        corresponding eigenvectors (or zero vectors) so that P\*self =
-        D\*P.
+        Return matrices `D` and `P`, where `D` is a diagonal matrix of
+        eigenvalues and the rows of `P` are corresponding eigenvectors
+        (or zero vectors).
+
+        INPUT:
+
+        - ``other`` -- a square matrix `B` (default: ``None``) in a generalized
+          eigenvalue problem; if ``None``, an ordinary eigenvalue problem is
+          solved
+
+        OUTPUT:
+
+        If ``self`` is a square matrix `A`, then the output is a diagonal
+        matrix `D` and a matrix `P` such that
+
+        .. MATH::
+
+            P A = D P,
+
+        where the rows of `P` are eigenvectors of `A` and the diagonal entries
+        of `D` are the corresponding eigenvalues.
+
+        If a matrix `B` is passed as optional argument, the output is a
+        solution to the generalized eigenvalue problem such that
+
+        .. MATH::
+
+            P A = D P B.
+
+        The ordinary eigenvalue problem is equivalent to the generalized one if
+        `B` is the identity matrix.
+
+        The generalized eigenvector decomposition is currently only implemented
+        for matrices over ``RDF`` and ``CDF``.
 
         EXAMPLES::
 
@@ -6541,14 +6635,14 @@ cdef class Matrix(Matrix1):
             sage: P*A == D*P
             True
 
-        Because P is invertible, A is diagonalizable.
+        Because `P` is invertible, `A` is diagonalizable.
 
         ::
 
             sage: A == (~P)*D*P
             True
 
-        The matrix P may contain zero rows corresponding to eigenvalues for
+        The matrix `P` may contain zero rows corresponding to eigenvalues for
         which the algebraic multiplicity is greater than the geometric
         multiplicity. In these cases, the matrix is not diagonalizable.
 
@@ -6558,7 +6652,6 @@ cdef class Matrix(Matrix1):
             [2 1 0]
             [0 2 1]
             [0 0 2]
-            sage: A = jordan_block(2,3)
             sage: D, P = A.eigenmatrix_left()
             sage: D
             [2 0 0]
@@ -6571,6 +6664,42 @@ cdef class Matrix(Matrix1):
             sage: P*A == D*P
             True
 
+        A generalized eigenvector decomposition::
+
+            sage: A = matrix(RDF, [[1, -2], [3, 4]])
+            sage: B = matrix(RDF, [[0, 7], [2, -3]])
+            sage: D, P = A.eigenmatrix_left(B)
+            sage: (P * A - D * P * B).norm() < 1e-14
+            True
+
+        The matrix `B` in a generalized eigenvalue problem may be singular::
+
+            sage: A = matrix.identity(CDF, 2)
+            sage: B = matrix(CDF, [[2, 1+I], [4, 2+2*I]])
+            sage: D, P = A.eigenmatrix_left(B)
+            sage: D.diagonal()  # tol 1e-14
+            [0.2 - 0.1*I, +infinity]
+
+        In this case, we can still verify the eigenvector equation for the
+        first eigenvalue and first eigenvector::
+
+            sage: l = D[0, 0]
+            sage: v = P[0, :]
+            sage: (v * A - l * v * B).norm() < 1e-14
+            True
+
+        The second eigenvector is contained in the left kernel of `B`::
+
+            sage: (P[1, :] * B).norm() < 1e-14
+            True
+
+        .. SEEALSO::
+
+            :meth:`eigenvalues`,
+            :meth:`eigenvectors_left`,
+            :meth:`.Matrix_double_dense.eigenvectors_left`,
+            :meth:`eigenmatrix_right`.
+
         TESTS:
 
         For matrices with floating point entries, some platforms will
@@ -6630,7 +6759,7 @@ cdef class Matrix(Matrix1):
         """
         from sage.misc.flatten import flatten
         from sage.matrix.constructor import diagonal_matrix, matrix
-        evecs = self.eigenvectors_left()
+        evecs = self.eigenvectors_left(other=other)
         D = diagonal_matrix(flatten([[e[0]]*e[2] for e in evecs]))
         rows = []
         for e in evecs:
@@ -6647,12 +6776,42 @@ cdef class Matrix(Matrix1):
 
     left_eigenmatrix = eigenmatrix_left
 
-    def eigenmatrix_right(self):
+    def eigenmatrix_right(self, other=None):
         r"""
-        Return matrices D and P, where D is a diagonal matrix of
-        eigenvalues and P is the corresponding matrix where the columns are
-        corresponding eigenvectors (or zero vectors) so that self\*P =
-        P\*D.
+        Return matrices `D` and `P`, where `D` is a diagonal matrix of
+        eigenvalues and the columns of `P` are corresponding eigenvectors
+        (or zero vectors).
+
+        INPUT:
+
+        - ``other`` -- a square matrix `B` (default: ``None``) in a generalized
+          eigenvalue problem; if ``None``, an ordinary eigenvalue problem is
+          solved
+
+        OUTPUT:
+
+        If ``self`` is a square matrix `A`, then the output is a diagonal
+        matrix `D` and a matrix `P` such that
+
+        .. MATH::
+
+            A P = P D,
+
+        where the columns of `P` are eigenvectors of `A` and the diagonal
+        entries of `D` are the corresponding eigenvalues.
+
+        If a matrix `B` is passed as optional argument, the output is a
+        solution to the generalized eigenvalue problem such that
+
+        .. MATH::
+
+            A P = B P D.
+
+        The ordinary eigenvalue problem is equivalent to the generalized one if
+        `B` is the identity matrix.
+
+        The generalized eigenvector decomposition is currently only implemented
+        for matrices over ``RDF`` and ``CDF``.
 
         EXAMPLES::
 
@@ -6672,14 +6831,14 @@ cdef class Matrix(Matrix1):
             sage: A*P == P*D
             True
 
-        Because P is invertible, A is diagonalizable.
+        Because `P` is invertible, `A` is diagonalizable.
 
         ::
 
             sage: A == P*D*(~P)
             True
 
-        The matrix P may contain zero columns corresponding to eigenvalues
+        The matrix `P` may contain zero columns corresponding to eigenvalues
         for which the algebraic multiplicity is greater than the geometric
         multiplicity. In these cases, the matrix is not diagonalizable.
 
@@ -6689,7 +6848,6 @@ cdef class Matrix(Matrix1):
             [2 1 0]
             [0 2 1]
             [0 0 2]
-            sage: A = jordan_block(2,3)
             sage: D, P = A.eigenmatrix_right()
             sage: D
             [2 0 0]
@@ -6702,6 +6860,42 @@ cdef class Matrix(Matrix1):
             sage: A*P == P*D
             True
 
+        A generalized eigenvector decomposition::
+
+            sage: A = matrix(RDF, [[1, -2], [3, 4]])
+            sage: B = matrix(RDF, [[0, 7], [2, -3]])
+            sage: D, P = A.eigenmatrix_right(B)
+            sage: (A * P - B * P * D).norm() < 1e-14
+            True
+
+        The matrix `B` in a generalized eigenvalue problem may be singular::
+
+            sage: A = matrix.identity(RDF, 2)
+            sage: B = matrix(RDF, [[3, 5], [6, 10]])
+            sage: D, P = A.eigenmatrix_right(B); D   # tol 1e-14
+            [0.07692307692307694                 0.0]
+            [                0.0           +infinity]
+
+        In this case, we can still verify the eigenvector equation for the
+        first eigenvalue and first eigenvector::
+
+            sage: l = D[0, 0]
+            sage: v = P[:, 0]
+            sage: (A * v  - B * v * l).norm() < 1e-14
+            True
+
+        The second eigenvector is contained in the right kernel of `B`::
+
+            sage: (B * P[:, 1]).norm() < 1e-14
+            True
+
+        .. SEEALSO::
+
+            :meth:`eigenvalues`,
+            :meth:`eigenvectors_right`,
+            :meth:`.Matrix_double_dense.eigenvectors_right`,
+            :meth:`eigenmatrix_left`.
+
         TESTS:
 
         For matrices with floating point entries, some platforms will
@@ -6744,7 +6938,8 @@ cdef class Matrix(Matrix1):
             True
 
         """
-        D,P = self.transpose().eigenmatrix_left()
+        D,P = self.transpose().eigenmatrix_left(None if other is None
+                                                else other.transpose())
         return D,P.transpose()
 
     right_eigenmatrix = eigenmatrix_right