]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_element.py
eja: more orthonormalization fixes.
[sage.d.git] / mjo / eja / eja_element.py
index 120870b7fa08b1ee0fcb936b92c62a9470a27754..d3e9a33ceba6e1fc4e58b417f3b953e2e2d1c3d7 100644 (file)
@@ -2,10 +2,10 @@ from sage.matrix.constructor import matrix
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 
-from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
+from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 from mjo.eja.eja_utils import _mat2vec
 
-class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
+class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
     """
     An element of a Euclidean Jordan algebra.
     """
@@ -227,9 +227,9 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         Ditto for the quaternions::
 
-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)
             sage: J.one().inner_product(J.one())
-            3
+            2
 
         TESTS:
 
@@ -340,6 +340,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  TrivialEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
             ....:                                  random_eja)
 
         EXAMPLES::
@@ -387,6 +389,37 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: x,y = J.random_elements(2)
             sage: (x*y).det() == x.det()*y.det()
             True
+
+        The determinant in matrix algebras is just the usual determinant::
+
+            sage: set_random_seed()
+            sage: X = matrix.random(QQ,3)
+            sage: X = X + X.T
+            sage: J1 = RealSymmetricEJA(3)
+            sage: J2 = RealSymmetricEJA(3,field=QQ,orthonormalize=False)
+            sage: expected = X.det()
+            sage: actual1 = J1(X).det()
+            sage: actual2 = J2(X).det()
+            sage: actual1 == expected
+            True
+            sage: actual2 == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J1 = ComplexHermitianEJA(2)
+            sage: J2 = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: X = matrix.random(GaussianIntegers(), 2)
+            sage: X = X + X.H
+            sage: expected = AA(X.det())
+            sage: actual1 = J1(J1.real_embed(X)).det()
+            sage: actual2 = J2(J2.real_embed(X)).det()
+            sage: expected == actual1
+            True
+            sage: expected == actual2
+            True
+
         """
         P = self.parent()
         r = P.rank()
@@ -411,8 +444,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
-        We appeal to the quadratic representation as in Koecher's
-        Theorem 12 in Chapter III, Section 5.
+        In general we appeal to the quadratic representation as in
+        Koecher's Theorem 12 in Chapter III, Section 5. But if the
+        parent algebra's "characteristic polynomial of" coefficients
+        happen to be cached, then we use Proposition II.2.4 in Faraut
+        and Korányi which gives a formula for the inverse based on the
+        characteristic polynomial and the Cayley-Hamilton theorem for
+        Euclidean Jordan algebras::
 
         SETUP::
 
@@ -482,22 +520,19 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             ....:    x.operator().inverse()(J.one()) == x.inverse() )
             True
 
-        Proposition II.2.4 in Faraut and Korányi gives a formula for
-        the inverse based on the characteristic polynomial and the
-        Cayley-Hamilton theorem for Euclidean Jordan algebras::
+        Check that the fast (cached) and slow algorithms give the same
+        answer::
 
-            sage: set_random_seed()
-            sage: J = ComplexHermitianEJA(3)
-            sage: x = J.random_element()
-            sage: while not x.is_invertible():
-            ....:     x = J.random_element()
-            sage: r = J.rank()
-            sage: a = x.characteristic_polynomial().coefficients(sparse=False)
-            sage: expected  = (-1)^(r+1)/x.det()
-            sage: expected *= sum( a[i+1]*x^i for i in range(r) )
-            sage: x.inverse() == expected
+            sage: set_random_seed()                              # long time
+            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+            sage: x = J.random_element()                         # long time
+            sage: while not x.is_invertible():                   # long time
+            ....:     x = J.random_element()                     # long time
+            sage: slow = x.inverse()                             # long time
+            sage: _ = J._charpoly_coefficients()                 # long time
+            sage: fast = x.inverse()                             # long time
+            sage: slow == fast                                   # long time
             True
-
         """
         if not self.is_invertible():
             raise ValueError("element is not invertible")
@@ -523,7 +558,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         zero, but we need the characteristic polynomial for the
         determinant. The minimal polynomial is a lot easier to get,
         so we use Corollary 2 in Chapter V of Koecher to check
-        whether or not the paren't algebra's zero element is a root
+        whether or not the parent algebra's zero element is a root
         of this element's minimal polynomial.
 
         That is... unless the coefficients of our algebra's
@@ -554,6 +589,18 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: (not J.is_trivial()) and J.zero().is_invertible()
             False
 
+        Test that the fast (cached) and slow algorithms give the same
+        answer::
+
+            sage: set_random_seed()                              # long time
+            sage: J = random_eja(field=QQ, orthonormalize=False) # long time
+            sage: x = J.random_element()                         # long time
+            sage: slow = x.is_invertible()                       # long time
+            sage: _ = J._charpoly_coefficients()                 # long time
+            sage: fast = x.is_invertible()                       # long time
+            sage: slow == fast                                   # long time
+            True
+
         """
         if self.is_zero():
             if self.parent().is_trivial():
@@ -783,10 +830,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
-        For now, we skip the messy minimal polynomial computation
-        and instead return the dimension of the vector space spanned
-        by the powers of this element. The latter is a bit more
-        straightforward to compute.
+        .........
 
         SETUP::
 
@@ -834,12 +878,59 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
-        if self.is_zero() and not self.parent().is_trivial():
+        n = self.parent().dimension()
+
+        if n == 0:
+            # The minimal polynomial is an empty product, i.e. the
+            # constant polynomial "1" having degree zero.
+            return 0
+        elif self.is_zero():
             # The minimal polynomial of zero in a nontrivial algebra
-            # is "t"; in a trivial algebra it's "1" by convention
-            # (it's an empty product).
+            # is "t", and is of degree one.
+            return 1
+        elif n == 1:
+            # If this is a nonzero element of a nontrivial algebra, it
+            # has degree at least one. It follows that, in an algebra
+            # of dimension one, the degree must be actually one.
             return 1
-        return self.subalgebra_generated_by().dimension()
+
+        # BEWARE: The subalgebra_generated_by() method uses the result
+        # of this method to construct a basis for the subalgebra. That
+        # means, in particular, that we cannot implement this method
+        # as ``self.subalgebra_generated_by().dimension()``.
+
+        # Algorithm: keep appending (vector representations of) powers
+        # self as rows to a matrix and echelonizing it. When its rank
+        # stops increasing, we've reached a redundancy.
+
+        # Given the special cases above, we can assume that "self" is
+        # nonzero, the algebra is nontrivial, and that its dimension
+        # is at least two.
+        M = matrix([(self.parent().one()).to_vector()])
+        old_rank = 1
+
+        # Specifying the row-reduction algorithm can e.g.  help over
+        # AA because it avoids the RecursionError that gets thrown
+        # when we have to look too hard for a root.
+        #
+        # Beware: QQ supports an entirely different set of "algorithm"
+        # keywords than do AA and RR.
+        algo = None
+        from sage.rings.all import QQ
+        if self.parent().base_ring() is not QQ:
+            algo = "scaled_partial_pivoting"
+
+        for d in range(1,n):
+            M = matrix(M.rows() + [(self**d).to_vector()])
+            M.echelonize(algo)
+            new_rank = M.rank()
+            if new_rank == old_rank:
+                return new_rank
+            else:
+                old_rank = new_rank
+
+        return n
+
 
 
     def left_matrix(self):
@@ -945,7 +1036,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: n_max = RealSymmetricEJA._max_random_instance_size()
             sage: n = ZZ.random_element(1, n_max)
             sage: J1 = RealSymmetricEJA(n)
-            sage: J2 = RealSymmetricEJA(n,normalize_basis=False)
+            sage: J2 = RealSymmetricEJA(n,orthonormalize=False)
             sage: X = random_matrix(AA,n)
             sage: X = X*X.transpose()
             sage: x1 = J1(X)
@@ -966,20 +1057,21 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
                 # in the "normal" case without us having to think about it.
                 return self.operator().minimal_polynomial()
 
-        A = self.subalgebra_generated_by(orthonormalize_basis=False)
+        A = self.subalgebra_generated_by(orthonormalize=False)
         return A(self).operator().minimal_polynomial()
 
 
 
-    def natural_representation(self):
+    def to_matrix(self):
         """
-        Return a more-natural representation of this element.
+        Return an (often more natural) representation of this element as a
+        matrix.
 
-        Every finite-dimensional Euclidean Jordan Algebra is a
-        direct sum of five simple algebras, four of which comprise
-        Hermitian matrices. This method returns the original
-        "natural" representation of this element as a Hermitian
-        matrix, if it has one. If not, you get the usual representation.
+        Every finite-dimensional Euclidean Jordan Algebra is a direct
+        sum of five simple algebras, four of which comprise Hermitian
+        matrices. This method returns a "natural" matrix
+        representation of this element as either a Hermitian matrix or
+        column vector.
 
         SETUP::
 
@@ -991,7 +1083,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: J = ComplexHermitianEJA(3)
             sage: J.one()
             e0 + e3 + e8
-            sage: J.one().natural_representation()
+            sage: J.one().to_matrix()
             [1 0 0 0 0 0]
             [0 1 0 0 0 0]
             [0 0 1 0 0 0]
@@ -1001,31 +1093,27 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         ::
 
-            sage: J = QuaternionHermitianEJA(3)
+            sage: J = QuaternionHermitianEJA(2)
             sage: J.one()
-            e0 + e5 + e14
-            sage: J.one().natural_representation()
-            [1 0 0 0 0 0 0 0 0 0 0 0]
-            [0 1 0 0 0 0 0 0 0 0 0 0]
-            [0 0 1 0 0 0 0 0 0 0 0 0]
-            [0 0 0 1 0 0 0 0 0 0 0 0]
-            [0 0 0 0 1 0 0 0 0 0 0 0]
-            [0 0 0 0 0 1 0 0 0 0 0 0]
-            [0 0 0 0 0 0 1 0 0 0 0 0]
-            [0 0 0 0 0 0 0 1 0 0 0 0]
-            [0 0 0 0 0 0 0 0 1 0 0 0]
-            [0 0 0 0 0 0 0 0 0 1 0 0]
-            [0 0 0 0 0 0 0 0 0 0 1 0]
-            [0 0 0 0 0 0 0 0 0 0 0 1]
+            e0 + e5
+            sage: J.one().to_matrix()
+            [1 0 0 0 0 0 0 0]
+            [0 1 0 0 0 0 0 0]
+            [0 0 1 0 0 0 0 0]
+            [0 0 0 1 0 0 0 0]
+            [0 0 0 0 1 0 0 0]
+            [0 0 0 0 0 1 0 0]
+            [0 0 0 0 0 0 1 0]
+            [0 0 0 0 0 0 0 1]
 
         """
-        B = self.parent().natural_basis()
-        W = self.parent().natural_basis_space()
+        B = self.parent().matrix_basis()
+        W = self.parent().matrix_space()
 
         # This is just a manual "from_vector()", but of course
         # matrix spaces aren't vector spaces in sage, so they
         # don't have a from_vector() method.
-        return W.linear_combination(zip(B,self.to_vector()))
+        return W.linear_combination( zip(B, self.to_vector()) )
 
 
     def norm(self):
@@ -1078,10 +1166,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         P = self.parent()
         left_mult_by_self = lambda y: self*y
         L = P.module_morphism(function=left_mult_by_self, codomain=P)
-        return FiniteDimensionalEuclideanJordanAlgebraOperator(
-                 P,
-                 P,
-                 L.matrix() )
+        return FiniteDimensionalEJAOperator(P, P, L.matrix() )
 
 
     def quadratic_representation(self, other=None):
@@ -1273,13 +1358,13 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             [(0, f2), (1, f0)]
 
         """
-        A = self.subalgebra_generated_by(orthonormalize_basis=True)
+        A = self.subalgebra_generated_by(orthonormalize=True)
         result = []
         for (evalue, proj) in A(self).operator().spectral_decomposition():
             result.append( (evalue, proj(A.one()).superalgebra_element()) )
         return result
 
-    def subalgebra_generated_by(self, orthonormalize_basis=False):
+    def subalgebra_generated_by(self, **kwargs):
         """
         Return the associative subalgebra of the parent EJA generated
         by this element.
@@ -1326,8 +1411,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         """
-        from mjo.eja.eja_element_subalgebra import FiniteDimensionalEuclideanJordanElementSubalgebra
-        return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)
+        from mjo.eja.eja_element_subalgebra import FiniteDimensionalEJAElementSubalgebra
+        return FiniteDimensionalEJAElementSubalgebra(self, **kwargs)
 
 
     def subalgebra_idempotent(self):