README: rewrite it, it was rather out-of-date

[sage.d.git] / mjo / eja / eja_element.py

diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py
index e1ea609efce334a9800ce162ba74d1342662da41..f4d5995ce9302d2fdc518dac2b0ca20f0fadf6d1 100644 (file)
--- a/mjo/eja/eja_element.py
+++ b/mjo/eja/eja_element.py
@@ -3,11 +3,11 @@ from sage.misc.cachefunc import cached_method
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 
 from sage.modules.free_module import VectorSpace
 from sage.modules.with_basis.indexed_element import IndexedFreeModuleElement
 

-from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
+from mjo.eja.eja_operator import EJAOperator

 from mjo.eja.eja_utils import _scale
 
 
 from mjo.eja.eja_utils import _scale
 
 

-class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
+class EJAElement(IndexedFreeModuleElement):

     """
     An element of a Euclidean Jordan algebra.
     """
     """
     An element of a Euclidean Jordan algebra.
     """


@@ -332,7 +332,8 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (AlbertEJA,
+            ....:                                  JordanSpinEJA,

             ....:                                  TrivialEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,
             ....:                                  TrivialEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  ComplexHermitianEJA,


@@ -395,6 +396,22 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: actual2 == expected
             True
 
             sage: actual2 == expected
             True
 

+        There's a formula for the determinant of the Albert algebra
+        (Yokota, Section 2.1)::
+
+            sage: def albert_det(x):
+            ....:     X = x.to_matrix()
+            ....:     res  = X[0,0]*X[1,1]*X[2,2]
+            ....:     res += 2*(X[1,2]*X[2,0]*X[0,1]).real()
+            ....:     res -= X[0,0]*X[1,2]*X[2,1]
+            ....:     res -= X[1,1]*X[2,0]*X[0,2]
+            ....:     res -= X[2,2]*X[0,1]*X[1,0]
+            ....:     return res.leading_coefficient()
+            sage: J = AlbertEJA(field=QQ, orthonormalize=False)
+            sage: xs = J.random_elements(10)
+            sage: all( albert_det(x) == x.det() for x in xs )
+            True
+

         """
         P = self.parent()
         r = P.rank()
         """
         P = self.parent()
         r = P.rank()


@@ -538,7 +555,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         If computing my determinant will be fast, we do so and compare
         with zero (Proposition II.2.4 in Faraut and
 
         If computing my determinant will be fast, we do so and compare
         with zero (Proposition II.2.4 in Faraut and

-        Koranyi). Otherwise, Proposition II.3.2 in Faraut and Koranyi
+        Korányi). Otherwise, Proposition II.3.2 in Faraut and Korányi

         reduces the problem to the invertibility of my quadratic
         representation.
 
         reduces the problem to the invertibility of my quadratic
         representation.
 


@@ -795,7 +812,23 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
 
         ALGORITHM:
 
 
         ALGORITHM:
 

-        .........
+        First we handle the special cases where the algebra is
+        trivial, this element is zero, or the dimension of the algebra
+        is one and this element is not zero. With those out of the
+        way, we may assume that ``self`` is nonzero, the algebra is
+        nontrivial, and that the dimension of the algebra is at least
+        two.
+
+        Beginning with the algebra's unit element (power zero), we add
+        successive (basis representations of) powers of this element
+        to a matrix, row-reducing at each step. After row-reducing, we
+        check the rank of the matrix. If adding a row and row-reducing
+        does not increase the rank of the matrix at any point, the row
+        we've just added lives in the span of the previous ones; thus
+        the corresponding power of ``self`` lives in the span of its
+        lesser powers. When that happens, the degree of the minimal
+        polynomial is the rank of the matrix; if it never happens, the
+        degree must be the dimension of the entire space.

 
         SETUP::
 
 
         SETUP::
 


@@ -838,7 +871,6 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: x = random_eja().random_element()
             sage: x.degree() == x.minimal_polynomial().degree()
             True
             sage: x = random_eja().random_element()
             sage: x.degree() == x.minimal_polynomial().degree()
             True

-

         """
         n = self.parent().dimension()
 
         """
         n = self.parent().dimension()
 


@@ -1147,7 +1179,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         P = self.parent()
         left_mult_by_self = lambda y: self*y
         L = P.module_morphism(function=left_mult_by_self, codomain=P)
         P = self.parent()
         left_mult_by_self = lambda y: self*y
         L = P.module_morphism(function=left_mult_by_self, codomain=P)

-        return FiniteDimensionalEJAOperator(P, P, L.matrix() )
+        return EJAOperator(P, P, L.matrix() )

 
 
     def quadratic_representation(self, other=None):
 
 
     def quadratic_representation(self, other=None):


@@ -1440,7 +1472,10 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         if self.is_nilpotent():
             raise ValueError("this only works with non-nilpotent elements!")
 
         if self.is_nilpotent():
             raise ValueError("this only works with non-nilpotent elements!")
 

-        J = self.subalgebra_generated_by()
+        # The subalgebra is transient (we return an element of the
+        # superalgebra, i.e. this algebra) so why bother
+        # orthonormalizing?
+        J = self.subalgebra_generated_by(orthonormalize=False)

         u = J(self)
 
         # The image of the matrix of left-u^m-multiplication
         u = J(self)
 
         # The image of the matrix of left-u^m-multiplication


@@ -1461,14 +1496,12 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         # subspace... or do we? Can't we just solve, knowing that
         # A(c) = u^(s+1) should have a solution in the big space,
         # too?
         # subspace... or do we? Can't we just solve, knowing that
         # A(c) = u^(s+1) should have a solution in the big space,
         # too?

-        #
-        # Beware, solve_right() means that we're using COLUMN vectors.
-        # Our FiniteDimensionalAlgebraElement superclass uses rows.

         u_next = u**(s+1)
         A = u_next.operator().matrix()
         c = J.from_vector(A.solve_right(u_next.to_vector()))
 
         u_next = u**(s+1)
         A = u_next.operator().matrix()
         c = J.from_vector(A.solve_right(u_next.to_vector()))
 

-        # Now c is the idempotent we want, but it still lives in the subalgebra.
+        # Now c is the idempotent we want, but it still lives in
+        # the subalgebra.

         return c.superalgebra_element()
 
 
         return c.superalgebra_element()
 
 


@@ -1542,7 +1575,75 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         # we want the negative of THAT for the trace.
         return -p(*self.to_vector())
 
         # we want the negative of THAT for the trace.
         return -p(*self.to_vector())
 

-    def operator_inner_product(self, other):
+
+    def trace_inner_product(self, other):
+        """
+        Return the trace inner product of myself and ``other``.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import random_eja
+
+        TESTS:
+
+        The trace inner product is commutative, bilinear, and associative::
+
+            sage: J = random_eja()
+            sage: x,y,z = J.random_elements(3)
+            sage: # commutative
+            sage: x.trace_inner_product(y) == y.trace_inner_product(x)
+            True
+            sage: # bilinear
+            sage: a = J.base_ring().random_element()
+            sage: actual = (a*(x+z)).trace_inner_product(y)
+            sage: expected = ( a*x.trace_inner_product(y) +
+            ....:              a*z.trace_inner_product(y) )
+            sage: actual == expected
+            True
+            sage: actual = x.trace_inner_product(a*(y+z))
+            sage: expected = ( a*x.trace_inner_product(y) +
+            ....:              a*x.trace_inner_product(z) )
+            sage: actual == expected
+            True
+            sage: # associative
+            sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
+            True
+
+        """
+        if not other in self.parent():
+            raise TypeError("'other' must live in the same algebra")
+
+        return (self*other).trace()
+
+
+    def trace_norm(self):
+        """
+        The norm of this element with respect to :meth:`trace_inner_product`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  HadamardEJA)
+
+        EXAMPLES::
+
+            sage: J = HadamardEJA(2)
+            sage: x = sum(J.gens())
+            sage: x.trace_norm()
+            1.414213562373095?
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: x = sum(J.gens())
+            sage: x.trace_norm()
+            2.828427124746190?
+
+        """
+        return self.trace_inner_product(self).sqrt()
+
+
+    def operator_trace_inner_product(self, other):

         r"""
         Return the operator inner product of myself and ``other``.
 
         r"""
         Return the operator inner product of myself and ``other``.
 


@@ -1554,10 +1655,10 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         Euclidean Jordan algebra, this is another associative inner
         product under which the cone of squares is symmetric.
 
         Euclidean Jordan algebra, this is another associative inner
         product under which the cone of squares is symmetric.
 

-        This *probably* works even if the basis hasn't been
-        orthonormalized because the eigenvalues of the corresponding
-        matrix don't change when the basis does (they're preserved by
-        any similarity transformation).
+        This works even if the basis hasn't been orthonormalized
+        because the eigenvalues of the corresponding matrix don't
+        change when the basis does (they're preserved by any
+        similarity transformation).

 
         SETUP::
 
 
         SETUP::
 


@@ -1577,7 +1678,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: x,y = J.random_elements(2)
             sage: n = J.dimension()
             sage: r = J.rank()
             sage: x,y = J.random_elements(2)
             sage: n = J.dimension()
             sage: r = J.rank()

-            sage: actual = x.operator_inner_product(y)
+            sage: actual = x.operator_trace_inner_product(y)

             sage: expected = (n/r)*x.trace_inner_product(y)
             sage: actual == expected
             True
             sage: expected = (n/r)*x.trace_inner_product(y)
             sage: actual == expected
             True


@@ -1588,7 +1689,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: x,y = J.random_elements(2)
             sage: n = J.dimension()
             sage: r = J.rank()
             sage: x,y = J.random_elements(2)
             sage: n = J.dimension()
             sage: r = J.rank()

-            sage: actual = x.operator_inner_product(y)
+            sage: actual = x.operator_trace_inner_product(y)

             sage: expected = (n/r)*x.trace_inner_product(y)
             sage: actual == expected
             True
             sage: expected = (n/r)*x.trace_inner_product(y)
             sage: actual == expected
             True


@@ -1599,7 +1700,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: x,y = J.random_elements(2)
             sage: n = J.dimension()
             sage: r = J.rank()
             sage: x,y = J.random_elements(2)
             sage: n = J.dimension()
             sage: r = J.rank()

-            sage: actual = x.operator_inner_product(y)
+            sage: actual = x.operator_trace_inner_product(y)

             sage: expected = (n/r)*x.trace_inner_product(y)
             sage: actual == expected
             True
             sage: expected = (n/r)*x.trace_inner_product(y)
             sage: actual == expected
             True


@@ -1612,113 +1713,94 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: # commutative
             sage: J = random_eja()
             sage: x,y,z = J.random_elements(3)
             sage: # commutative

-            sage: x.operator_inner_product(y) == y.operator_inner_product(x)
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = y.operator_trace_inner_product(x)
+            sage: actual == expected

             True
             sage: # bilinear
             sage: a = J.base_ring().random_element()
             True
             sage: # bilinear
             sage: a = J.base_ring().random_element()

-            sage: actual = (a*(x+z)).operator_inner_product(y)
-            sage: expected = ( a*x.operator_inner_product(y) +
-            ....:              a*z.operator_inner_product(y) )
+            sage: actual = (a*(x+z)).operator_trace_inner_product(y)
+            sage: expected = ( a*x.operator_trace_inner_product(y) +
+            ....:              a*z.operator_trace_inner_product(y) )

             sage: actual == expected
             True
             sage: actual == expected
             True

-            sage: actual = x.operator_inner_product(a*(y+z))
-            sage: expected = ( a*x.operator_inner_product(y) +
-            ....:              a*x.operator_inner_product(z) )
+            sage: actual = x.operator_trace_inner_product(a*(y+z))
+            sage: expected = ( a*x.operator_trace_inner_product(y) +
+            ....:              a*x.operator_trace_inner_product(z) )

             sage: actual == expected
             True
             sage: # associative
             sage: actual == expected
             True
             sage: # associative

-            sage: actual = (x*y).operator_inner_product(z)
-            sage: expected = y.operator_inner_product(x*z)
+            sage: actual = (x*y).operator_trace_inner_product(z)
+            sage: expected = y.operator_trace_inner_product(x*z)

             sage: actual == expected
             True
 
             sage: actual == expected
             True
 

-        """
-        if not other in self.parent():
-            raise TypeError("'other' must live in the same algebra")
-
-        return (self*other).operator().matrix().trace()
-
+        Despite the fact that the implementation uses a matrix representation,
+        the answer is independent of the basis used::

 
 

-    def trace_inner_product(self, other):
-        """
-        Return the trace inner product of myself and ``other``.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import random_eja
-
-        TESTS:
-
-        The trace inner product is commutative, bilinear, and associative::
-
-            sage: J = random_eja()
-            sage: x,y,z = J.random_elements(3)
-            sage: # commutative
-            sage: x.trace_inner_product(y) == y.trace_inner_product(x)
-            True
-            sage: # bilinear
-            sage: a = J.base_ring().random_element()
-            sage: actual = (a*(x+z)).trace_inner_product(y)
-            sage: expected = ( a*x.trace_inner_product(y) +
-            ....:              a*z.trace_inner_product(y) )
-            sage: actual == expected
-            True
-            sage: actual = x.trace_inner_product(a*(y+z))
-            sage: expected = ( a*x.trace_inner_product(y) +
-            ....:              a*x.trace_inner_product(z) )
+            sage: J = RealSymmetricEJA(3, field=QQ, orthonormalize=False)
+            sage: V = RealSymmetricEJA(3)
+            sage: x,y = J.random_elements(2)
+            sage: w = V(x.to_matrix())
+            sage: z = V(y.to_matrix())
+            sage: expected = x.operator_trace_inner_product(y)
+            sage: actual = w.operator_trace_inner_product(z)

             sage: actual == expected
             True
             sage: actual == expected
             True

-            sage: # associative
-            sage: (x*y).trace_inner_product(z) == y.trace_inner_product(x*z)
-            True

 
         """
         if not other in self.parent():
             raise TypeError("'other' must live in the same algebra")
 
 
         """
         if not other in self.parent():
             raise TypeError("'other' must live in the same algebra")
 

-        return (self*other).trace()
+        return (self*other).operator().matrix().trace()

 
 
 
 

-    def trace_norm(self):
+    def operator_trace_norm(self):

         """
         """

-        The norm of this element with respect to :meth:`trace_inner_product`.
+        The norm of this element with respect to
+        :meth:`operator_trace_inner_product`.

 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  HadamardEJA)
 
 
         SETUP::
 
             sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
             ....:                                  HadamardEJA)
 

-        EXAMPLES::
+        EXAMPLES:
+
+        On a simple algebra, this will differ from :meth:`trace_norm`
+        by the scalar factor ``(n/r).sqrt()``, where `n` is the
+        dimension of the algebra and `r` its rank. This follows from
+        the corresponding result (Proposition III.4.2 of Faraut and
+        Korányi) for the trace inner product::

 
             sage: J = HadamardEJA(2)
             sage: x = sum(J.gens())
 
             sage: J = HadamardEJA(2)
             sage: x = sum(J.gens())

-            sage: x.trace_norm()
+            sage: x.operator_trace_norm()

             1.414213562373095?
 
         ::
 
             sage: J = JordanSpinEJA(4)
             sage: x = sum(J.gens())
             1.414213562373095?
 
         ::
 
             sage: J = JordanSpinEJA(4)
             sage: x = sum(J.gens())

-            sage: x.trace_norm()
-            2.828427124746190?
+            sage: x.operator_trace_norm()
+            4

 
         """
 
         """

-        return self.trace_inner_product(self).sqrt()
+        return self.operator_trace_inner_product(self).sqrt()

 
 
 
 

-class CartesianProductEJAElement(FiniteDimensionalEJAElement):
-    def det(self):
-        r"""
-        Compute the determinant of this product-element using the
-        determianants of its factors.
-
-        This result Follows from the spectral decomposition of (say)
-        the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
-        0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
-        """
-        from sage.misc.misc_c import prod
-        return prod( f.det() for f in self.cartesian_factors() )
+class CartesianProductParentEJAElement(EJAElement):
+    r"""
+    An intermediate class for elements that have a Cartesian
+    product as their parent algebra.

 
 

+    This is needed because the ``to_matrix`` method (which gives you a
+    representation from the superalgebra) needs to do special stuff
+    for Cartesian products. Specifically, an EJA subalgebra of a
+    Cartesian product EJA will not itself be a Cartesian product (it
+    has its own basis) -- but we want ``to_matrix()`` to be able to
+    give us a Cartesian product representation.
+    """

     def to_matrix(self):
         # An override is necessary to call our custom _scale().
         B = self.parent().matrix_basis()
     def to_matrix(self):
         # An override is necessary to call our custom _scale().
         B = self.parent().matrix_basis()


@@ -1726,7 +1808,20 @@ class CartesianProductEJAElement(FiniteDimensionalEJAElement):
 
         # Aaaaand linear combinations don't work in Cartesian
         # product spaces, even though they provide a method with
 
         # Aaaaand linear combinations don't work in Cartesian
         # product spaces, even though they provide a method with

-        # that name. This is hidden behind an "if" because the
+        # that name. This is hidden in a subclass because the

         # _scale() function is slow.
         pairs = zip(B, self.to_vector())
         return W.sum( _scale(b, alpha) for (b,alpha) in pairs )
         # _scale() function is slow.
         pairs = zip(B, self.to_vector())
         return W.sum( _scale(b, alpha) for (b,alpha) in pairs )

+
+class CartesianProductEJAElement(CartesianProductParentEJAElement):
+    def det(self):
+        r"""
+        Compute the determinant of this product-element using the
+        determianants of its factors.
+
+        This result Follows from the spectral decomposition of (say)
+        the pair `(x,y)` in terms of the Jordan frame `\left\{ (c_1,
+        0),(c_2, 0),...,(0,d_1),(0,d_2),... \right\}.
+        """
+        from sage.misc.misc_c import prod
+        return prod( f.det() for f in self.cartesian_factors() )