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 0983ea0607a438464b3e9ab6245a344ba59a5aef..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()
 
 


@@ -1561,7 +1594,7 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
             True
             sage: # bilinear
             sage: x.trace_inner_product(y) == y.trace_inner_product(x)
             True
             sage: # bilinear

-            sage: a = J.base_ring().random_element();
+            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 = (a*(x+z)).trace_inner_product(y)
             sage: expected = ( a*x.trace_inner_product(y) +
             ....:              a*z.trace_inner_product(y) )


@@ -1610,19 +1643,164 @@ class FiniteDimensionalEJAElement(IndexedFreeModuleElement):
         return self.trace_inner_product(self).sqrt()
 
 
         return self.trace_inner_product(self).sqrt()
 
 

-class CartesianProductEJAElement(FiniteDimensionalEJAElement):
-    def det(self):
+    def operator_trace_inner_product(self, other):

         r"""
         r"""

-        Compute the determinant of this product-element using the
-        determianants of its factors.
+        Return the operator inner product of myself and ``other``.
+
+        The "operator inner product," whose name is not standard, is
+        defined be the usual linear-algebraic trace of the
+        ``(x*y).operator()``.
+
+        Proposition III.1.5 in Faraut and Korányi shows that on any
+        Euclidean Jordan algebra, this is another associative inner
+        product under which the cone of squares is symmetric.
+
+        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::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  random_eja)
+
+        EXAMPLES:
+
+        Proposition III.4.2 of Faraut and Korányi shows that on a
+        simple algebra of rank `r` and dimension `n`, this inner
+        product is `n/r` times the canonical
+        :meth:`trace_inner_product`::
+
+            sage: J = JordanSpinEJA(4, field=QQ)
+            sage: x,y = J.random_elements(2)
+            sage: n = J.dimension()
+            sage: r = J.rank()
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = (n/r)*x.trace_inner_product(y)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: J = RealSymmetricEJA(3)
+            sage: x,y = J.random_elements(2)
+            sage: n = J.dimension()
+            sage: r = J.rank()
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = (n/r)*x.trace_inner_product(y)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False)
+            sage: x,y = J.random_elements(2)
+            sage: n = J.dimension()
+            sage: r = J.rank()
+            sage: actual = x.operator_trace_inner_product(y)
+            sage: expected = (n/r)*x.trace_inner_product(y)
+            sage: actual == expected
+            True
+
+        TESTS:
+
+        The operator inner product is commutative, bilinear, and
+        associative::
+
+            sage: J = random_eja()
+            sage: x,y,z = J.random_elements(3)
+            sage: # commutative
+            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()
+            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 = 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 = (x*y).operator_trace_inner_product(z)
+            sage: expected = y.operator_trace_inner_product(x*z)
+            sage: actual == expected
+            True
+
+        Despite the fact that the implementation uses a matrix representation,
+        the answer is independent of the basis used::
+
+            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

 
 

-        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() )
+        if not other in self.parent():
+            raise TypeError("'other' must live in the same algebra")
+
+        return (self*other).operator().matrix().trace()
+
+
+    def operator_trace_norm(self):
+        """
+        The norm of this element with respect to
+        :meth:`operator_trace_inner_product`.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            ....:                                  HadamardEJA)
+
+        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: x.operator_trace_norm()
+            1.414213562373095?
+
+        ::
+
+            sage: J = JordanSpinEJA(4)
+            sage: x = sum(J.gens())
+            sage: x.operator_trace_norm()
+            4
+
+        """
+        return self.operator_trace_inner_product(self).sqrt()
+
+
+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()


@@ -1630,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() )