]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: improve projection/inclusion implementation for DirectSumEJA.
[sage.d.git] / mjo / eja / eja_algebra.py
index c40c8be4c371255f9bc174f99862ef0a7a204555..7857190a2e458f0d955a49de8df7e7a9052bfe4a 100644 (file)
@@ -3,6 +3,17 @@ Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
 are used in optimization, and have some additional nice methods beyond
 what can be supported in a general Jordan Algebra.
+
+
+SETUP::
+
+    sage: from mjo.eja.eja_algebra import random_eja
+
+EXAMPLES::
+
+    sage: random_eja()
+    Euclidean Jordan algebra of dimension...
+
 """
 
 from itertools import repeat
@@ -14,7 +25,6 @@ from sage.matrix.constructor import matrix
 from sage.matrix.matrix_space import MatrixSpace
 from sage.misc.cachefunc import cached_method
 from sage.misc.lazy_import import lazy_import
-from sage.misc.prandom import choice
 from sage.misc.table import table
 from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
@@ -23,6 +33,7 @@ from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
 from mjo.eja.eja_element import FiniteDimensionalEuclideanJordanAlgebraElement
 lazy_import('mjo.eja.eja_subalgebra',
             'FiniteDimensionalEuclideanJordanSubalgebra')
+from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 from mjo.eja.eja_utils import _mat2vec
 
 class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
@@ -216,25 +227,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         coords =  W.coordinate_vector(_mat2vec(elt))
         return self.from_vector(coords)
 
-    @staticmethod
-    def _max_random_instance_size():
-        """
-        Return an integer "size" that is an upper bound on the size of
-        this algebra when it is used in a random test
-        case. Unfortunately, the term "size" is quite vague -- when
-        dealing with `R^n` under either the Hadamard or Jordan spin
-        product, the "size" refers to the dimension `n`. When dealing
-        with a matrix algebra (real symmetric or complex/quaternion
-        Hermitian), it refers to the size of the matrix, which is
-        far less than the dimension of the underlying vector space.
-
-        We default to five in this class, which is safe in `R^n`. The
-        matrix algebra subclasses (or any class where the "size" is
-        interpreted to be far less than the dimension) should override
-        with a smaller number.
-        """
-        raise NotImplementedError
-
     def _repr_(self):
         """
         Return a string representation of ``self``.
@@ -842,16 +834,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return tuple( self.random_element(thorough)
                       for idx in range(count) )
 
-    @classmethod
-    def random_instance(cls, field=AA, **kwargs):
-        """
-        Return a random instance of this type of algebra.
-
-        Beware, this will crash for "most instances" because the
-        constructor below looks wrong.
-        """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
-        return cls(n, field, **kwargs)
 
     @cached_method
     def _charpoly_coefficients(self):
@@ -1013,32 +995,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
 
-
-def random_eja(field=AA):
-    """
-    Return a "random" finite-dimensional Euclidean Jordan Algebra.
-
-    SETUP::
-
-        sage: from mjo.eja.eja_algebra import random_eja
-
-    TESTS::
-
-        sage: random_eja()
-        Euclidean Jordan algebra of dimension...
-
-    """
-    classname = choice([TrivialEJA,
-                        HadamardEJA,
-                        JordanSpinEJA,
-                        RealSymmetricEJA,
-                        ComplexHermitianEJA,
-                        QuaternionHermitianEJA])
-    return classname.random_instance(field=field)
-
-
-
-
 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
     Algebras whose basis consists of vectors with rational
@@ -1098,12 +1054,72 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
         return tuple(map(lambda x: x.change_ring(self.base_ring()), a))
 
 
-class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
+class ConcreteEuclideanJordanAlgebra:
+    r"""
+    A class for the Euclidean Jordan algebras that we know by name.
+
+    These are the Jordan algebras whose basis, multiplication table,
+    rank, and so on are known a priori. More to the point, they are
+    the Euclidean Jordan algebras for which we are able to conjure up
+    a "random instance."
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import ConcreteEuclideanJordanAlgebra
+
+    TESTS:
+
+    Our natural basis is normalized with respect to the natural inner
+    product unless we specify otherwise::
+
+        sage: set_random_seed()
+        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: all( b.norm() == 1 for b in J.gens() )
+        True
+
+    Since our natural basis is normalized with respect to the natural
+    inner product, and since we know that this algebra is an EJA, any
+    left-multiplication operator's matrix will be symmetric because
+    natural->EJA basis representation is an isometry and within the EJA
+    the operator is self-adjoint by the Jordan axiom::
+
+        sage: set_random_seed()
+        sage: J = ConcreteEuclideanJordanAlgebra.random_instance()
+        sage: x = J.random_element()
+        sage: x.operator().matrix().is_symmetric()
+        True
+
+    """
+
     @staticmethod
     def _max_random_instance_size():
-        # Play it safe, since this will be squared and the underlying
-        # field can have dimension 4 (quaternions) too.
-        return 2
+        """
+        Return an integer "size" that is an upper bound on the size of
+        this algebra when it is used in a random test
+        case. Unfortunately, the term "size" is ambiguous -- when
+        dealing with `R^n` under either the Hadamard or Jordan spin
+        product, the "size" refers to the dimension `n`. When dealing
+        with a matrix algebra (real symmetric or complex/quaternion
+        Hermitian), it refers to the size of the matrix, which is far
+        less than the dimension of the underlying vector space.
+
+        This method must be implemented in each subclass.
+        """
+        raise NotImplementedError
+
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+
+        This method should be implemented in each subclass.
+        """
+        from sage.misc.prandom import choice
+        eja_class = choice(cls.__subclasses__())
+        return eja_class.random_instance(field)
+
+
+class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
     def __init__(self, field, basis, normalize_basis=True, **kwargs):
         """
@@ -1287,7 +1303,8 @@ class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return M
 
 
-class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
+class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
+                       ConcreteEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of real symmetric n-by-n
     matrices, the usual symmetric Jordan product, and the trace inner
@@ -1346,25 +1363,6 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
         sage: RealSymmetricEJA(3, prefix='q').gens()
         (q0, q1, q2, q3, q4, q5)
 
-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
-
-        sage: set_random_seed()
-        sage: J = RealSymmetricEJA.random_instance()
-        sage: all( b.norm() == 1 for b in J.gens() )
-        True
-
-    Since our natural basis is normalized with respect to the natural
-    inner product, and since we know that this algebra is an EJA, any
-    left-multiplication operator's matrix will be symmetric because
-    natural->EJA basis representation is an isometry and within the EJA
-    the operator is self-adjoint by the Jordan axiom::
-
-        sage: set_random_seed()
-        sage: x = RealSymmetricEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
-
     We can construct the (trivial) algebra of rank zero::
 
         sage: RealSymmetricEJA(0)
@@ -1407,6 +1405,13 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra):
     def _max_random_instance_size():
         return 4 # Dimension 10
 
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
         basis = self._denormalized_basis(n, field)
@@ -1450,8 +1455,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()
-            sage: n_max = ComplexMatrixEuclideanJordanAlgebra._max_random_instance_size()
-            sage: n = ZZ.random_element(n_max)
+            sage: n = ZZ.random_element(3)
             sage: F = QuadraticField(-1, 'I')
             sage: X = random_matrix(F, n)
             sage: Y = random_matrix(F, n)
@@ -1576,7 +1580,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
 
 
-class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
+class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
+                          ConcreteEuclideanJordanAlgebra):
     """
     The rank-n simple EJA consisting of complex Hermitian n-by-n
     matrices over the real numbers, the usual symmetric Jordan product,
@@ -1627,25 +1632,6 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
         sage: ComplexHermitianEJA(2, prefix='z').gens()
         (z0, z1, z2, z3)
 
-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
-
-        sage: set_random_seed()
-        sage: J = ComplexHermitianEJA.random_instance()
-        sage: all( b.norm() == 1 for b in J.gens() )
-        True
-
-    Since our natural basis is normalized with respect to the natural
-    inner product, and since we know that this algebra is an EJA, any
-    left-multiplication operator's matrix will be symmetric because
-    natural->EJA basis representation is an isometry and within the EJA
-    the operator is self-adjoint by the Jordan axiom::
-
-        sage: set_random_seed()
-        sage: x = ComplexHermitianEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
-
     We can construct the (trivial) algebra of rank zero::
 
         sage: ComplexHermitianEJA(0)
@@ -1715,6 +1701,17 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra):
                                                  **kwargs)
         self.rank.set_cache(n)
 
+    @staticmethod
+    def _max_random_instance_size():
+        return 3 # Dimension 9
+
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
@@ -1746,8 +1743,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         Embedding is a homomorphism (isomorphism, in fact)::
 
             sage: set_random_seed()
-            sage: n_max = QuaternionMatrixEuclideanJordanAlgebra._max_random_instance_size()
-            sage: n = ZZ.random_element(n_max)
+            sage: n = ZZ.random_element(2)
             sage: Q = QuaternionAlgebra(QQ,-1,-1)
             sage: X = random_matrix(Q, n)
             sage: Y = random_matrix(Q, n)
@@ -1879,8 +1875,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
 
 
-class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
-    """
+class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
+                             ConcreteEuclideanJordanAlgebra):
+    r"""
     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
     matrices, the usual symmetric Jordan product, and the
     real-part-of-trace inner product. It has dimension `2n^2 - n` over
@@ -1930,25 +1927,6 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
         sage: QuaternionHermitianEJA(2, prefix='a').gens()
         (a0, a1, a2, a3, a4, a5)
 
-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
-
-        sage: set_random_seed()
-        sage: J = QuaternionHermitianEJA.random_instance()
-        sage: all( b.norm() == 1 for b in J.gens() )
-        True
-
-    Since our natural basis is normalized with respect to the natural
-    inner product, and since we know that this algebra is an EJA, any
-    left-multiplication operator's matrix will be symmetric because
-    natural->EJA basis representation is an isometry and within the EJA
-    the operator is self-adjoint by the Jordan axiom::
-
-        sage: set_random_seed()
-        sage: x = QuaternionHermitianEJA.random_instance().random_element()
-        sage: x.operator().matrix().is_symmetric()
-        True
-
     We can construct the (trivial) algebra of rank zero::
 
         sage: QuaternionHermitianEJA(0)
@@ -2019,8 +1997,24 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra):
                                                     **kwargs)
         self.rank.set_cache(n)
 
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum rank of a random QuaternionHermitianEJA.
+        """
+        return 2 # Dimension 6
 
-class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)
+
+
+class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
+                  ConcreteEuclideanJordanAlgebra):
     """
     Return the Euclidean Jordan Algebra corresponding to the set
     `R^n` under the Hadamard product.
@@ -2078,8 +2072,20 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
 
     @staticmethod
     def _max_random_instance_size():
+        r"""
+        The maximum dimension of a random HadamardEJA.
+        """
         return 5
 
+    @classmethod
+    def random_instance(cls, field=AA, **kwargs):
+        """
+        Return a random instance of this type of algebra.
+        """
+        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        return cls(n, field, **kwargs)
+
+
     def inner_product(self, x, y):
         """
         Faster to reimplement than to use natural representations.
@@ -2105,7 +2111,8 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra):
         return x.to_vector().inner_product(y.to_vector())
 
 
-class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
+class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
+                      ConcreteEuclideanJordanAlgebra):
     r"""
     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
     with the half-trace inner product and jordan product ``x*y =
@@ -2222,6 +2229,9 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra):
 
     @staticmethod
     def _max_random_instance_size():
+        r"""
+        The maximum dimension of a random BilinearFormEJA.
+        """
         return 5
 
     @classmethod
@@ -2334,6 +2344,13 @@ class JordanSpinEJA(BilinearFormEJA):
         B = matrix.identity(field, n)
         super(JordanSpinEJA, self).__init__(B, field, **kwargs)
 
+    @staticmethod
+    def _max_random_instance_size():
+        r"""
+        The maximum dimension of a random JordanSpinEJA.
+        """
+        return 5
+
     @classmethod
     def random_instance(cls, field=AA, **kwargs):
         """
@@ -2345,7 +2362,8 @@ class JordanSpinEJA(BilinearFormEJA):
         return cls(n, field, **kwargs)
 
 
-class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra):
+class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
+                 ConcreteEuclideanJordanAlgebra):
     """
     The trivial Euclidean Jordan algebra consisting of only a zero element.
 
@@ -2401,7 +2419,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     SETUP::
 
-        sage: from mjo.eja.eja_algebra import (HadamardEJA,
+        sage: from mjo.eja.eja_algebra import (random_eja,
+        ....:                                  HadamardEJA,
         ....:                                  RealSymmetricEJA,
         ....:                                  DirectSumEJA)
 
@@ -2415,8 +2434,25 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         sage: J.rank()
         5
 
+    TESTS:
+
+    The external direct sum construction is only valid when the two factors
+    have the same base ring; an error is raised otherwise::
+
+        sage: set_random_seed()
+        sage: J1 = random_eja(AA)
+        sage: J2 = random_eja(QQ)
+        sage: J = DirectSumEJA(J1,J2)
+        Traceback (most recent call last):
+        ...
+        ValueError: algebras must share the same base field
+
     """
-    def __init__(self, J1, J2, field=AA, **kwargs):
+    def __init__(self, J1, J2, **kwargs):
+        if J1.base_ring() != J2.base_ring():
+            raise ValueError("algebras must share the same base field")
+        field = J1.base_ring()
+
         self._factors = (J1, J2)
         n1 = J1.dimension()
         n2 = J2.dimension()
@@ -2489,8 +2525,11 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         """
         (J1,J2) = self.factors()
         n = J1.dimension()
-        pi_left  = lambda x: J1.from_vector(x.to_vector()[:n])
-        pi_right = lambda x: J2.from_vector(x.to_vector()[n:])
+        V_basis = self.vector_space().basis()
+        P1 = matrix(self.base_ring(), V_basis[:n])
+        P2 = matrix(self.base_ring(), V_basis[n:])
+        pi_left = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J1,P1)
+        pi_right = FiniteDimensionalEuclideanJordanAlgebraOperator(self,J2,P2)
         return (pi_left, pi_right)
 
     def inclusions(self):
@@ -2499,7 +2538,8 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  JordanSpinEJA,
             ....:                                  RealSymmetricEJA,
             ....:                                  DirectSumEJA)
 
@@ -2524,14 +2564,35 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
             sage: J.one().to_vector()
             (1, 0, 0, 1, 0, 1)
 
+        TESTS:
+
+        Composing a projection with the corresponding inclusion should
+        produce the identity map, and mismatching them should produce
+        the zero map::
+
+            sage: set_random_seed()
+            sage: J1 = random_eja()
+            sage: J2 = random_eja()
+            sage: J = DirectSumEJA(J1,J2)
+            sage: (iota_left, iota_right) = J.inclusions()
+            sage: (pi_left, pi_right) = J.projections()
+            sage: pi_left*iota_left == J1.one().operator()
+            True
+            sage: pi_right*iota_right == J2.one().operator()
+            True
+            sage: (pi_left*iota_right).is_zero()
+            True
+            sage: (pi_right*iota_left).is_zero()
+            True
+
         """
         (J1,J2) = self.factors()
         n = J1.dimension()
         V_basis = self.vector_space().basis()
         I1 = matrix.column(self.base_ring(), V_basis[:n])
         I2 = matrix.column(self.base_ring(), V_basis[n:])
-        iota_left = lambda x: self.from_vector(I1*x.to_vector())
-        iota_right = lambda x: self.from_vector(I2*+x.to_vector())
+        iota_left = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,self,I1)
+        iota_right = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,self,I2)
         return (iota_left, iota_right)
 
     def inner_product(self, x, y):
@@ -2550,7 +2611,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
         EXAMPLE::
 
-            sage: J1 = HadamardEJA(3)
+            sage: J1 = HadamardEJA(3,QQ)
             sage: J2 = QuaternionHermitianEJA(2,QQ,normalize_basis=False)
             sage: J = DirectSumEJA(J1,J2)
             sage: x1 = J1.one()
@@ -2572,3 +2633,7 @@ class DirectSumEJA(FiniteDimensionalEuclideanJordanAlgebra):
         y2 = pi_right(y)
 
         return (x1.inner_product(y1) + x2.inner_product(y2))
+
+
+
+random_eja = ConcreteEuclideanJordanAlgebra.random_instance