]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: add brute-force associativity test.
[sage.d.git] / mjo / eja / eja_algebra.py
index d38ba87d37d3e6bdb01e23c7a9631bf98ee2c96c..3b64ddd7374d8931f80fdefe4b0fb1734b682995 100644 (file)
@@ -1,9 +1,53 @@
 """
-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.
-
+Representations and constructions for Euclidean Jordan algebras.
+
+A Euclidean Jordan algebra is a Jordan algebra that has some
+additional properties:
+
+  1.   It is finite-dimensional.
+  2.   Its scalar field is the real numbers.
+  3a.  An inner product is defined on it, and...
+  3b.  That inner product is compatible with the Jordan product
+       in the sense that `<x*y,z> = <y,x*z>` for all elements
+       `x,y,z` in the algebra.
+
+Every Euclidean Jordan algebra is formally-real: for any two elements
+`x` and `y` in the algebra, `x^{2} + y^{2} = 0` implies that `x = y =
+0`. Conversely, every finite-dimensional formally-real Jordan algebra
+can be made into a Euclidean Jordan algebra with an appropriate choice
+of inner-product.
+
+Formally-real Jordan algebras were originally studied as a framework
+for quantum mechanics. Today, Euclidean Jordan algebras are crucial in
+symmetric cone optimization, since every symmetric cone arises as the
+cone of squares in some Euclidean Jordan algebra.
+
+It is known that every Euclidean Jordan algebra decomposes into an
+orthogonal direct sum (essentially, a Cartesian product) of simple
+algebras, and that moreover, up to Jordan-algebra isomorphism, there
+are only five families of simple algebras. We provide constructions
+for these simple algebras:
+
+  * :class:`BilinearFormEJA`
+  * :class:`RealSymmetricEJA`
+  * :class:`ComplexHermitianEJA`
+  * :class:`QuaternionHermitianEJA`
+
+Missing from this list is the algebra of three-by-three octononion
+Hermitian matrices, as there is (as of yet) no implementation of the
+octonions in SageMath. In addition to these, we provide two other
+example constructions,
+
+  * :class:`HadamardEJA`
+  * :class:`TrivialEJA`
+
+The Jordan spin algebra is a bilinear form algebra where the bilinear
+form is the identity. The Hadamard EJA is simply a Cartesian product
+of one-dimensional spin algebras. And last but not least, the trivial
+EJA is exactly what you think. Cartesian products of these are also
+supported using the usual ``cartesian_product()`` function; as a
+result, we support (up to isomorphism) all Euclidean Jordan algebras
+that don't involve octonions.
 
 SETUP::
 
@@ -13,7 +57,6 @@ EXAMPLES::
 
     sage: random_eja()
     Euclidean Jordan algebra of dimension...
-
 """
 
 from itertools import repeat
@@ -41,25 +84,34 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     INPUT:
 
-      - basis -- a tuple of basis elements in "matrix form," which
-        must be the same form as the arguments to ``jordan_product``
-        and ``inner_product``. In reality, "matrix form" can be either
-        vectors, matrices, or a Cartesian product (ordered tuple)
-        of vectors or matrices. All of these would ideally be vector
-        spaces in sage with no special-casing needed; but in reality
-        we turn vectors into column-matrices and Cartesian products
-        `(a,b)` into column matrices `(a,b)^{T}` after converting
-        `a` and `b` themselves.
-
-      - jordan_product -- function of two ``basis`` elements (in
-        matrix form) that returns their jordan product, also in matrix
-        form; this will be applied to ``basis`` to compute a
-        multiplication table for the algebra.
-
-      - inner_product -- function of two ``basis`` elements (in matrix
-        form) that returns their inner product. This will be applied
-        to ``basis`` to compute an inner-product table (basically a
-        matrix) for this algebra.
+      - ``basis`` -- a tuple; a tuple of basis elements in "matrix
+        form," which must be the same form as the arguments to
+        ``jordan_product`` and ``inner_product``. In reality, "matrix
+        form" can be either vectors, matrices, or a Cartesian product
+        (ordered tuple) of vectors or matrices. All of these would
+        ideally be vector spaces in sage with no special-casing
+        needed; but in reality we turn vectors into column-matrices
+        and Cartesian products `(a,b)` into column matrices
+        `(a,b)^{T}` after converting `a` and `b` themselves.
+
+      - ``jordan_product`` -- a function; afunction of two ``basis``
+        elements (in matrix form) that returns their jordan product,
+        also in matrix form; this will be applied to ``basis`` to
+        compute a multiplication table for the algebra.
+
+      - ``inner_product`` -- a function; a function of two ``basis``
+        elements (in matrix form) that returns their inner
+        product. This will be applied to ``basis`` to compute an
+        inner-product table (basically a matrix) for this algebra.
+
+      - ``field`` -- a subfield of the reals (default: ``AA``); the scalar
+        field for the algebra.
+
+      - ``orthonormalize`` -- boolean (default: ``True``); whether or
+        not to orthonormalize the basis. Doing so is expensive and
+        generally rules out using the rationals as your ``field``, but
+        is required for spectral decompositions.
+
     """
     Element = FiniteDimensionalEJAElement
 
@@ -126,7 +178,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
 
         category = MagmaticAlgebras(field).FiniteDimensional()
-        category = category.WithBasis().Unital()
+        category = category.WithBasis().Unital().Commutative()
+
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
@@ -370,6 +423,18 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         """
         return "Associative" in self.category().axioms()
 
+    def _is_commutative(self):
+        r"""
+        Whether or not this algebra's multiplication table is commutative.
+
+        This method should of course always return ``True``, unless
+        this algebra was constructed with ``check_axioms=False`` and
+        passed an invalid multiplication table.
+        """
+        return all( self.product_on_basis(i,j) == self.product_on_basis(i,j)
+                    for i in range(self.dimension())
+                    for j in range(self.dimension()) )
+
     def _is_jordanian(self):
         r"""
         Whether or not this algebra's multiplication table respects the
@@ -377,7 +442,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
         We only check one arrangement of `x` and `y`, so for a
         ``True`` result to be truly true, you should also check
-        :meth:`is_commutative`. This method should of course always
+        :meth:`_is_commutative`. This method should of course always
         return ``True``, unless this algebra was constructed with
         ``check_axioms=False`` and passed an invalid multiplication table.
         """
@@ -387,6 +452,81 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                     for i in range(self.dimension())
                     for j in range(self.dimension()) )
 
+    def _jordan_product_is_associative(self):
+        r"""
+        Return whether or not this algebra's Jordan product is
+        associative; that is, whether or not `x*(y*z) = (x*y)*z`
+        for all `x,y,x`.
+
+        This method should agree with :meth:`is_associative` unless
+        you lied about the value of the ``associative`` parameter
+        when you constructed the algebra.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(4, orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = ComplexHermitianEJA(2,field=QQ,orthonormalize=False)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize=False)
+            sage: A._jordan_product_is_associative()
+            True
+
+        ::
+
+            sage: J = QuaternionHermitianEJA(2)
+            sage: J._jordan_product_is_associative()
+            False
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: A._jordan_product_is_associative()
+            True
+
+        """
+        R = self.base_ring()
+
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # I don't know of any examples that make this magnitude
+            # necessary because I don't know how to make an
+            # associative algebra when the element subalgebra
+            # construction is unreliable (as it is over RDF; we can't
+            # find the degree of an element because we can't compute
+            # the rank of a matrix). But even multiplication of floats
+            # is non-associative, so *some* epsilon is needed... let's
+            # just take the one from _inner_product_is_associative?
+            epsilon = 1e-15
+
+        for i in range(self.dimension()):
+            for j in range(self.dimension()):
+                for k in range(self.dimension()):
+                    x = self.gens()[i]
+                    y = self.gens()[j]
+                    z = self.gens()[k]
+                    diff = (x*y)*z - x*(y*z)
+
+                    if diff.norm() > epsilon:
+                        return False
+
+        return True
+
     def _inner_product_is_associative(self):
         r"""
         Return whether or not this algebra's inner product `B` is
@@ -396,11 +536,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         this algebra was constructed with ``check_axioms=False`` and
         passed an invalid Jordan or inner-product.
         """
+        R = self.base_ring()
 
-        # Used to check whether or not something is zero in an inexact
-        # ring. This number is sufficient to allow the construction of
-        # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
-        epsilon = 1e-16
+        # Used to check whether or not something is zero.
+        epsilon = R.zero()
+        if not R.is_exact():
+            # This choice is sufficient to allow the construction of
+            # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True.
+            epsilon = 1e-15
 
         for i in range(self.dimension()):
             for j in range(self.dimension()):
@@ -410,12 +553,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                     z = self.gens()[k]
                     diff = (x*y).inner_product(z) - x.inner_product(y*z)
 
-                    if self.base_ring().is_exact():
-                        if diff != 0:
-                            return False
-                    else:
-                        if diff.abs() > epsilon:
-                            return False
+                    if diff.abs() > epsilon:
+                        return False
 
         return True
 
@@ -1724,9 +1863,9 @@ class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: RealSymmetricEJA(2, field=RDF)
+        sage: RealSymmetricEJA(2, field=RDF, check_axioms=True)
         Euclidean Jordan algebra of dimension 3 over Real Double Field
-        sage: RealSymmetricEJA(2, field=RR)
+        sage: RealSymmetricEJA(2, field=RR, check_axioms=True)
         Euclidean Jordan algebra of dimension 3 over Real Field with
         53 bits of precision
 
@@ -1992,9 +2131,9 @@ class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: ComplexHermitianEJA(2, field=RDF)
+        sage: ComplexHermitianEJA(2, field=RDF, check_axioms=True)
         Euclidean Jordan algebra of dimension 4 over Real Double Field
-        sage: ComplexHermitianEJA(2, field=RR)
+        sage: ComplexHermitianEJA(2, field=RR, check_axioms=True)
         Euclidean Jordan algebra of dimension 4 over Real Field with
         53 bits of precision
 
@@ -2289,9 +2428,9 @@ class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA):
 
     In theory, our "field" can be any subfield of the reals::
 
-        sage: QuaternionHermitianEJA(2, field=RDF)
+        sage: QuaternionHermitianEJA(2, field=RDF, check_axioms=True)
         Euclidean Jordan algebra of dimension 6 over Real Double Field
-        sage: QuaternionHermitianEJA(2, field=RR)
+        sage: QuaternionHermitianEJA(2, field=RR, check_axioms=True)
         Euclidean Jordan algebra of dimension 6 over Real Field with
         53 bits of precision
 
@@ -3210,4 +3349,42 @@ class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct,
 
 FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA
 
+class RationalBasisCartesianProductEJA(CartesianProductEJA,
+                                       RationalBasisEJA):
+    r"""
+    A separate class for products of algebras for which we know a
+    rational basis.
+
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+        ....:                                  RealSymmetricEJA)
+
+    EXAMPLES:
+
+    This gives us fast characteristic polynomial computations in
+    product algebras, too::
+
+
+        sage: J1 = JordanSpinEJA(2)
+        sage: J2 = RealSymmetricEJA(3)
+        sage: J = cartesian_product([J1,J2])
+        sage: J.characteristic_polynomial_of().degree()
+        5
+        sage: J.rank()
+        5
+
+    """
+    def __init__(self, algebras, **kwargs):
+        CartesianProductEJA.__init__(self, algebras, **kwargs)
+
+        self._rational_algebra = None
+        if self.vector_space().base_field() is not QQ:
+            self._rational_algebra = cartesian_product([
+                r._rational_algebra for r in algebras
+            ])
+
+
+RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
+
 random_eja = ConcreteEJA.random_instance