]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: begin working on a new constructor.
[sage.d.git] / mjo / eja / eja_algebra.py
index 21718ee85e168d9ead63b061a3d623cb3ab6c8cc..44d83147bbeceb4acae52d573749545dec3c4b99 100644 (file)
@@ -24,15 +24,12 @@ from sage.combinat.free_module import CombinatorialFreeModule
 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.table import table
 from sage.modules.free_module import FreeModule, VectorSpace
 from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF,
                             PolynomialRing,
                             QuadraticField)
 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
 
@@ -57,7 +54,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J(1)
             Traceback (most recent call last):
             ...
-            ValueError: not a naturally-represented algebra element
+            ValueError: not an element of this algebra
 
         """
         return None
@@ -67,7 +64,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
                  mult_table,
                  prefix='e',
                  category=None,
-                 natural_basis=None,
+                 matrix_basis=None,
                  check_field=True,
                  check_axioms=True):
         """
@@ -118,7 +115,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             if not all( len(l) == n for l in mult_table ):
                 raise ValueError("multiplication table is not square")
 
-        self._natural_basis = natural_basis
+        self._matrix_basis = matrix_basis
 
         if category is None:
             category = MagmaticAlgebras(field).FiniteDimensional()
@@ -137,10 +134,15 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         # long run to have the multiplication table be in terms of
         # algebra elements. We do this after calling the superclass
         # constructor so that from_vector() knows what to do.
-        self._multiplication_table = [
-            list(map(lambda x: self.from_vector(x), ls))
-            for ls in mult_table
-        ]
+        self._multiplication_table = [ [ self.vector_space().zero()
+                                         for i in range(n) ]
+                                       for j in range(n) ]
+        # take advantage of symmetry
+        for i in range(n):
+            for j in range(i+1):
+                elt = self.from_vector(mult_table[i][j])
+                self._multiplication_table[i][j] = elt
+                self._multiplication_table[j][i] = elt
 
         if check_axioms:
             if not self._is_commutative():
@@ -152,7 +154,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     def _element_constructor_(self, elt):
         """
-        Construct an element of this algebra from its natural
+        Construct an element of this algebra from its vector or matrix
         representation.
 
         This gets called only after the parent element _call_ method
@@ -180,13 +182,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J(A)
             Traceback (most recent call last):
             ...
-            ArithmeticError: vector is not in free module
+            ValueError: not an element of this algebra
 
         TESTS:
 
         Ensure that we can convert any element of the two non-matrix
-        simple algebras (whose natural representations are their usual
-        vector representations) back and forth faithfully::
+        simple algebras (whose matrix representations are columns)
+        back and forth faithfully::
 
             sage: set_random_seed()
             sage: J = HadamardEJA.random_instance()
@@ -197,9 +199,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: x = J.random_element()
             sage: J(x.to_vector().column()) == x
             True
-
         """
-        msg = "not a naturally-represented algebra element"
+        msg = "not an element of this algebra"
         if elt == 0:
             # The superclass implementation of random_element()
             # needs to be able to coerce "0" into the algebra.
@@ -211,20 +212,23 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 
-        natural_basis = self.natural_basis()
-        basis_space = natural_basis[0].matrix_space()
-        if elt not in basis_space:
+        if elt not in self.matrix_space():
             raise ValueError(msg)
 
         # Thanks for nothing! Matrix spaces aren't vector spaces in
-        # Sage, so we have to figure out its natural-basis coordinates
+        # Sage, so we have to figure out its matrix-basis coordinates
         # ourselves. We use the basis space's ring instead of the
         # element's ring because the basis space might be an algebraic
         # closure whereas the base ring of the 3-by-3 identity matrix
         # could be QQ instead of QQbar.
-        V = VectorSpace(basis_space.base_ring(), elt.nrows()*elt.ncols())
-        W = V.span_of_basis( _mat2vec(s) for s in natural_basis )
-        coords =  W.coordinate_vector(_mat2vec(elt))
+        V = VectorSpace(self.base_ring(), elt.nrows()*elt.ncols())
+        W = V.span_of_basis( _mat2vec(s) for s in self.matrix_basis() )
+
+        try:
+            coords =  W.coordinate_vector(_mat2vec(elt))
+        except ArithmeticError:  # vector is not in free module
+            raise ValueError(msg)
+
         return self.from_vector(coords)
 
     def _repr_(self):
@@ -513,18 +517,33 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         return table(M, header_row=True, header_column=True, frame=True)
 
 
-    def natural_basis(self):
+    def matrix_basis(self):
         """
-        Return a more-natural representation of this algebra's basis.
+        Return an (often more natural) representation of this algebras
+        basis as an ordered tuple of matrices.
+
+        Every finite-dimensional Euclidean Jordan Algebra is a, up to
+        Jordan isomorphism, a direct sum of five simple
+        algebras---four of which comprise Hermitian matrices. And the
+        last type of algebra can of course be thought of as `n`-by-`1`
+        column matrices (ambiguusly called column vectors) to avoid
+        special cases. As a result, matrices (and column vectors) are
+        a natural representation format for Euclidean Jordan algebra
+        elements.
 
-        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" basis
-        for our underlying vector space. (Typically, the natural basis
-        is used to construct the multiplication table in the first place.)
+        But, when we construct an algebra from a basis of matrices,
+        those matrix representations are lost in favor of coordinate
+        vectors *with respect to* that basis. We could eventually
+        convert back if we tried hard enough, but having the original
+        representations handy is valuable enough that we simply store
+        them and return them from this method.
 
-        Note that this will always return a matrix. The standard basis
-        in `R^n` will be returned as `n`-by-`1` column matrices.
+        Why implement this for non-matrix algebras? Avoiding special
+        cases for the :class:`BilinearFormEJA` pays with simplicity in
+        its own right. But mainly, we would like to be able to assume
+        that elements of a :class:`DirectSumEJA` can be displayed
+        nicely, without having to have special classes for direct sums
+        one of whose components was a matrix algebra.
 
         SETUP::
 
@@ -536,7 +555,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = RealSymmetricEJA(2)
             sage: J.basis()
             Finite family {0: e0, 1: e1, 2: e2}
-            sage: J.natural_basis()
+            sage: J.matrix_basis()
             (
             [1 0]  [                  0 0.7071067811865475?]  [0 0]
             [0 0], [0.7071067811865475?                   0], [0 1]
@@ -547,36 +566,38 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: J = JordanSpinEJA(2)
             sage: J.basis()
             Finite family {0: e0, 1: e1}
-            sage: J.natural_basis()
+            sage: J.matrix_basis()
             (
             [1]  [0]
             [0], [1]
             )
-
         """
-        if self._natural_basis is None:
-            M = self.natural_basis_space()
+        if self._matrix_basis is None:
+            M = self.matrix_space()
             return tuple( M(b.to_vector()) for b in self.basis() )
         else:
-            return self._natural_basis
+            return self._matrix_basis
 
 
-    def natural_basis_space(self):
+    def matrix_space(self):
         """
-        Return the matrix space in which this algebra's natural basis
-        elements live.
+        Return the matrix space in which this algebra's elements live, if
+        we think of them as matrices (including column vectors of the
+        appropriate size).
 
         Generally this will be an `n`-by-`1` column-vector space,
         except when the algebra is trivial. There it's `n`-by-`n`
-        (where `n` is zero), to ensure that two elements of the
-        natural basis space (empty matrices) can be multiplied.
+        (where `n` is zero), to ensure that two elements of the matrix
+        space (empty matrices) can be multiplied.
+
+        Matrix algebras override this with something more useful.
         """
         if self.is_trivial():
             return MatrixSpace(self.base_ring(), 0)
-        elif self._natural_basis is None or len(self._natural_basis) == 0:
+        elif self._matrix_basis is None or len(self._matrix_basis) == 0:
             return MatrixSpace(self.base_ring(), self.dimension(), 1)
         else:
-            return self._natural_basis[0].matrix_space()
+            return self._matrix_basis[0].matrix_space()
 
 
     @cached_method
@@ -703,22 +724,22 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             Vector space of degree 6 and dimension 2...
             sage: J1
             Euclidean Jordan algebra of dimension 3...
-            sage: J0.one().natural_representation()
+            sage: J0.one().to_matrix()
             [0 0 0]
             [0 0 0]
             [0 0 1]
             sage: orig_df = AA.options.display_format
             sage: AA.options.display_format = 'radical'
-            sage: J.from_vector(J5.basis()[0]).natural_representation()
+            sage: J.from_vector(J5.basis()[0]).to_matrix()
             [          0           0 1/2*sqrt(2)]
             [          0           0           0]
             [1/2*sqrt(2)           0           0]
-            sage: J.from_vector(J5.basis()[1]).natural_representation()
+            sage: J.from_vector(J5.basis()[1]).to_matrix()
             [          0           0           0]
             [          0           0 1/2*sqrt(2)]
             [          0 1/2*sqrt(2)           0]
             sage: AA.options.display_format = orig_df
-            sage: J1.one().natural_representation()
+            sage: J1.one().to_matrix()
             [1 0 0]
             [0 1 0]
             [0 0 0]
@@ -768,6 +789,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         if not c.is_idempotent():
             raise ValueError("element is not idempotent: %s" % c)
 
+        from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
+
         # Default these to what they should be if they turn out to be
         # trivial, because eigenspaces_left() won't return eigenvalues
         # corresponding to trivial spaces (e.g. it returns only the
@@ -971,38 +994,13 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
         Ensure that computing the rank actually works, since the ranks
         of all simple algebras are known and will be cached by default::
 
-            sage: J = HadamardEJA(4)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            4
-
-        ::
-
-            sage: J = JordanSpinEJA(4)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            2
-
-        ::
-
-            sage: J = RealSymmetricEJA(3)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            3
-
-        ::
-
-            sage: J = ComplexHermitianEJA(2)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            2
-
-        ::
+            sage: set_random_seed()    # long time
+            sage: J = random_eja()     # long time
+            sage: caches = J.rank()    # long time
+            sage: J.rank.clear_cache() # long time
+            sage: J.rank() == cached   # long time
+            True
 
-            sage: J = QuaternionHermitianEJA(2)
-            sage: J.rank.clear_cache()
-            sage: J.rank()
-            2
         """
         return len(self._charpoly_coefficients())
 
@@ -1027,6 +1025,75 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
     Element = FiniteDimensionalEuclideanJordanAlgebraElement
 
+class RationalBasisEuclideanJordanAlgebraNg(FiniteDimensionalEuclideanJordanAlgebra):
+    def __init__(self,
+                 field,
+                 basis,
+                 jordan_product,
+                 inner_product,
+                 orthonormalize=True,
+                 prefix='e',
+                 category=None,
+                 check_field=True,
+                 check_axioms=True):
+
+        n = len(basis)
+        vector_basis = basis
+
+        from sage.matrix.matrix import is_Matrix
+        basis_is_matrices = False
+
+        degree = 0
+        if n > 0:
+            if is_Matrix(basis[0]):
+                basis_is_matrices = True
+                vector_basis = tuple( map(_mat2vec,basis) )
+                degree = basis[0].nrows()**2
+            else:
+                degree = basis[0].degree()
+
+        V = VectorSpace(field, degree)
+
+        self._deorthonormalization_matrix = matrix.identity(field,n)
+        if orthonormalize:
+            A = matrix(field, vector_basis)
+            # uh oh, this is only the "usual" inner product
+            Q,R = A.gram_schmidt(orthonormal=True)
+            self._deorthonormalization_matrix = R.inverse().transpose()
+            vector_basis = Q.rows()
+            W = V.span_of_basis( vector_basis )
+            if basis_is_matrices:
+                from mjo.eja.eja_utils import _vec2mat
+                basis = tuple( map(_vec2mat,vector_basis) )
+
+        mult_table = [ [0 for i in range(n)] for j in range(n) ]
+        ip_table = [ [0 for i in range(n)] for j in range(n) ]
+
+        for i in range(n):
+            for j in range(i+1):
+                # do another mat2vec because the multiplication
+                # table is in terms of vectors
+                elt = _mat2vec(jordan_product(basis[i],basis[j]))
+                elt = W.coordinate_vector(elt)
+                mult_table[i][j] = elt
+                mult_table[j][i] = elt
+                ip = inner_product(basis[i],basis[j])
+                ip_table[i][j] = ip
+                ip_table[j][i] = ip
+
+        self._inner_product_matrix = matrix(field,ip_table)
+
+        if basis_is_matrices:
+            for m in basis:
+                m.set_immutable()
+
+        super().__init__(field,
+                         mult_table,
+                         prefix,
+                         category,
+                         basis,
+                         check_field,
+                         check_axioms)
 
 class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
     r"""
@@ -1072,7 +1139,7 @@ class RationalBasisEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebr
             return superclass._charpoly_coefficients()
 
         mult_table = tuple(
-            map(lambda x: x.to_vector(), ls)
+            tuple(map(lambda x: x.to_vector(), ls))
             for ls in self._multiplication_table
         )
 
@@ -1102,26 +1169,25 @@ class ConcreteEuclideanJordanAlgebra:
 
     TESTS:
 
-    Our natural basis is normalized with respect to the natural inner
-    product unless we specify otherwise::
+    Our basis is normalized with respect to the algebra's 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
+    Since our basis is orthonormal with respect to the algebra's 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::
+    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()
+        sage: x.operator().is_self_adjoint()
         True
-
     """
 
     @staticmethod
@@ -1183,7 +1249,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                 field = field.extension(p, 'sqrt2', embedding=RLF(2).sqrt())
                 basis = tuple( s.change_ring(field) for s in basis )
             self._basis_normalizers = tuple(
-                ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
+                ~(self.matrix_inner_product(s,s).sqrt()) for s in basis )
             basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
 
         # Now compute the multiplication and inner product tables.
@@ -1204,7 +1270,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                     # HACK: ignore the error here if we don't need the
                     # inner product (as is the case when we construct
                     # a dummy QQ-algebra for fast charpoly coefficients.
-                    ip_table[i][j] = self.natural_inner_product(basis[i],
+                    ip_table[i][j] = self.matrix_inner_product(basis[i],
                                                                 basis[j])
                 except:
                     pass
@@ -1217,7 +1283,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
 
         super(MatrixEuclideanJordanAlgebra, self).__init__(field,
                                                            mult_table,
-                                                           natural_basis=basis,
+                                                           matrix_basis=basis,
                                                            **kwargs)
 
         if algebra_dim == 0:
@@ -1240,7 +1306,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
             # entries in a nice field already. Just compute the thing.
             return super(MatrixEuclideanJordanAlgebra, self)._charpoly_coefficients()
 
-        basis = ( (b/n) for (b,n) in zip(self.natural_basis(),
+        basis = ( (b/n) for (b,n) in zip(self.matrix_basis(),
                                          self._basis_normalizers) )
 
         # Do this over the rationals and convert back at the end.
@@ -1300,7 +1366,7 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         raise NotImplementedError
 
     @classmethod
-    def natural_inner_product(cls,X,Y):
+    def matrix_inner_product(cls,X,Y):
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
@@ -1379,9 +1445,9 @@ class RealSymmetricEJA(RealMatrixEuclideanJordanAlgebra,
         sage: set_random_seed()
         sage: J = RealSymmetricEJA.random_instance()
         sage: x,y = J.random_elements(2)
-        sage: actual = (x*y).natural_representation()
-        sage: X = x.natural_representation()
-        sage: Y = y.natural_representation()
+        sage: actual = (x*y).to_matrix()
+        sage: X = x.to_matrix()
+        sage: Y = y.to_matrix()
         sage: expected = (X*Y + Y*X)/2
         sage: actual == expected
         True
@@ -1580,9 +1646,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @classmethod
-    def natural_inner_product(cls,X,Y):
+    def matrix_inner_product(cls,X,Y):
         """
-        Compute a natural inner product in this algebra directly from
+        Compute a matrix inner product in this algebra directly from
         its real embedding.
 
         SETUP::
@@ -1597,17 +1663,17 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: J = ComplexHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
-            sage: Xe = x.natural_representation()
-            sage: Ye = y.natural_representation()
+            sage: Xe = x.to_matrix()
+            sage: Ye = y.to_matrix()
             sage: X = ComplexHermitianEJA.real_unembed(Xe)
             sage: Y = ComplexHermitianEJA.real_unembed(Ye)
             sage: expected = (X*Y).trace().real()
-            sage: actual = ComplexHermitianEJA.natural_inner_product(Xe,Ye)
+            sage: actual = ComplexHermitianEJA.matrix_inner_product(Xe,Ye)
             sage: actual == expected
             True
 
         """
-        return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/2
+        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/2
 
 
 class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
@@ -1648,9 +1714,9 @@ class ComplexHermitianEJA(ComplexMatrixEuclideanJordanAlgebra,
         sage: set_random_seed()
         sage: J = ComplexHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
-        sage: actual = (x*y).natural_representation()
-        sage: X = x.natural_representation()
-        sage: Y = y.natural_representation()
+        sage: actual = (x*y).to_matrix()
+        sage: X = x.to_matrix()
+        sage: Y = y.to_matrix()
         sage: expected = (X*Y + Y*X)/2
         sage: actual == expected
         True
@@ -1875,9 +1941,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
     @classmethod
-    def natural_inner_product(cls,X,Y):
+    def matrix_inner_product(cls,X,Y):
         """
-        Compute a natural inner product in this algebra directly from
+        Compute a matrix inner product in this algebra directly from
         its real embedding.
 
         SETUP::
@@ -1892,17 +1958,17 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             sage: set_random_seed()
             sage: J = QuaternionHermitianEJA.random_instance()
             sage: x,y = J.random_elements(2)
-            sage: Xe = x.natural_representation()
-            sage: Ye = y.natural_representation()
+            sage: Xe = x.to_matrix()
+            sage: Ye = y.to_matrix()
             sage: X = QuaternionHermitianEJA.real_unembed(Xe)
             sage: Y = QuaternionHermitianEJA.real_unembed(Ye)
             sage: expected = (X*Y).trace().coefficient_tuple()[0]
-            sage: actual = QuaternionHermitianEJA.natural_inner_product(Xe,Ye)
+            sage: actual = QuaternionHermitianEJA.matrix_inner_product(Xe,Ye)
             sage: actual == expected
             True
 
         """
-        return RealMatrixEuclideanJordanAlgebra.natural_inner_product(X,Y)/4
+        return RealMatrixEuclideanJordanAlgebra.matrix_inner_product(X,Y)/4
 
 
 class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
@@ -1943,9 +2009,9 @@ class QuaternionHermitianEJA(QuaternionMatrixEuclideanJordanAlgebra,
         sage: set_random_seed()
         sage: J = QuaternionHermitianEJA.random_instance()
         sage: x,y = J.random_elements(2)
-        sage: actual = (x*y).natural_representation()
-        sage: X = x.natural_representation()
-        sage: Y = y.natural_representation()
+        sage: actual = (x*y).to_matrix()
+        sage: X = x.to_matrix()
+        sage: Y = y.to_matrix()
         sage: expected = (X*Y + Y*X)/2
         sage: actual == expected
         True