]> gitweb.michael.orlitzky.com - sage.d.git/commitdiff
eja: begin refactoring to allow noncanonical inner products.
authorMichael Orlitzky <michael@orlitzky.com>
Tue, 24 Nov 2020 16:34:58 +0000 (11:34 -0500)
committerMichael Orlitzky <michael@orlitzky.com>
Tue, 24 Nov 2020 16:34:58 +0000 (11:34 -0500)
mjo/eja/eja_algebra.py
mjo/eja/eja_subalgebra.py

index 6252cc29a729472ef7f182cd3393ab3db98b6b58..c569098b3d3e5d6fbceef6352e22861dd125a7f8 100644 (file)
@@ -385,7 +385,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
 
         SETUP::
 
-            sage: from mjo.eja.eja_algebra import random_eja
+            sage: from mjo.eja.eja_algebra import (random_eja,
+            ....:                                  HadamardEJA,
+            ....:                                  BilinearFormEJA)
 
         EXAMPLES:
 
@@ -398,10 +400,34 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             sage: (x*y).inner_product(z) == y.inner_product(x*z)
             True
 
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: J = HadamardEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            sage: actual = x.inner_product(y)
+            sage: expected = x.to_vector().inner_product(y.to_vector())
+            sage: actual == expected
+            True
+
+        Ensure that this is one-half of the trace inner-product in a
+        BilinearFormEJA that isn't just the reals (when ``n`` isn't
+        one). This is in Faraut and Koranyi, and also my "On the
+        symmetry..." paper::
+
+            sage: set_random_seed()
+            sage: J = BilinearFormEJA.random_instance()
+            sage: n = J.dimension()
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
+            True
         """
-        X = x.natural_representation()
-        Y = y.natural_representation()
-        return self.natural_inner_product(X,Y)
+        B = self._inner_product_matrix
+        return (B*x.to_vector()).inner_product(y.to_vector())
 
 
     def is_trivial(self):
@@ -531,20 +557,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             return self._natural_basis[0].matrix_space()
 
 
-    @staticmethod
-    def natural_inner_product(X,Y):
-        """
-        Compute the inner product of two naturally-represented elements.
-
-        For example in the real symmetric matrix EJA, this will compute
-        the trace inner-product of two n-by-n symmetric matrices. The
-        default should work for the real cartesian product EJA, the
-        Jordan spin EJA, and the real symmetric matrices. The others
-        will have to be overridden.
-        """
-        return (X.conjugate_transpose()*Y).trace()
-
-
     @cached_method
     def one(self):
         """
@@ -1135,6 +1147,10 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         basis = tuple(basis)
 
         algebra_dim = len(basis)
+        degree = 0 # size of the matrices
+        if algebra_dim > 0:
+            degree = basis[0].nrows()
+
         if algebra_dim > 1 and normalize_basis:
             # We'll need sqrt(2) to normalize the basis, and this
             # winds up in the multiplication table, so the whole
@@ -1149,10 +1165,37 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
                 ~(self.natural_inner_product(s,s).sqrt()) for s in basis )
             basis = tuple(s*c for (s,c) in zip(basis,self._basis_normalizers))
 
-        Qs = self.multiplication_table_from_matrix_basis(basis)
+        # Now compute the multiplication and inner product tables.
+        # We have to do this *after* normalizing the basis, because
+        # scaling affects the answers.
+        V = VectorSpace(field, degree**2)
+        W = V.span_of_basis( _mat2vec(s) for s in basis )
+        mult_table = [[W.zero() for j in range(algebra_dim)]
+                                for i in range(algebra_dim)]
+        ip_table = [[W.zero() for j in range(algebra_dim)]
+                              for i in range(algebra_dim)]
+        for i in range(algebra_dim):
+            for j in range(algebra_dim):
+                mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
+                mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
+
+                try:
+                    # 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],
+                                                                basis[j])
+                except:
+                    pass
+
+        try:
+            # HACK PART DEUX
+            self._inner_product_matrix = matrix(field,ip_table)
+        except:
+            pass
 
         super(MatrixEuclideanJordanAlgebra, self).__init__(field,
-                                                           Qs,
+                                                           mult_table,
                                                            natural_basis=basis,
                                                            **kwargs)
 
@@ -1211,39 +1254,6 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         return tuple( a_i.subs(subs_dict) for a_i in a )
 
 
-    @staticmethod
-    def multiplication_table_from_matrix_basis(basis):
-        """
-        At least three of the five simple Euclidean Jordan algebras have the
-        symmetric multiplication (A,B) |-> (AB + BA)/2, where the
-        multiplication on the right is matrix multiplication. Given a basis
-        for the underlying matrix space, this function returns a
-        multiplication table (obtained by looping through the basis
-        elements) for an algebra of those matrices.
-        """
-        # In S^2, for example, we nominally have four coordinates even
-        # though the space is of dimension three only. The vector space V
-        # is supposed to hold the entire long vector, and the subspace W
-        # of V will be spanned by the vectors that arise from symmetric
-        # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
-        if len(basis) == 0:
-            return []
-
-        field = basis[0].base_ring()
-        dimension = basis[0].nrows()
-
-        V = VectorSpace(field, dimension**2)
-        W = V.span_of_basis( _mat2vec(s) for s in basis )
-        n = len(basis)
-        mult_table = [[W.zero() for j in range(n)] for i in range(n)]
-        for i in range(n):
-            for j in range(n):
-                mat_entry = (basis[i]*basis[j] + basis[j]*basis[i])/2
-                mult_table[i][j] = W.coordinate_vector(_mat2vec(mat_entry))
-
-        return mult_table
-
-
     @staticmethod
     def real_embed(M):
         """
@@ -1268,7 +1278,6 @@ class MatrixEuclideanJordanAlgebra(FiniteDimensionalEuclideanJordanAlgebra):
         """
         raise NotImplementedError
 
-
     @classmethod
     def natural_inner_product(cls,X,Y):
         Xu = cls.real_unembed(X)
@@ -2059,6 +2068,16 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
         mult_table = [ [ V.gen(i)*(i == j) for j in range(n) ]
                        for i in range(n) ]
 
+        # Inner products are real numbers and not algebra
+        # elements, so once we turn the algebra element
+        # into a vector in inner_product(), we never go
+        # back. As a result -- contrary to what we do with
+        # self._multiplication_table -- we store the inner
+        # product table as a plain old matrix and not as
+        # an algebra operator.
+        ip_table = matrix.identity(field,n)
+        self._inner_product_matrix = ip_table
+
         super(HadamardEJA, self).__init__(field,
                                           mult_table,
                                           check_axioms=False,
@@ -2086,31 +2105,6 @@ class HadamardEJA(RationalBasisEuclideanJordanAlgebra,
         return cls(n, field, **kwargs)
 
 
-    def inner_product(self, x, y):
-        """
-        Faster to reimplement than to use natural representations.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import HadamardEJA
-
-        TESTS:
-
-        Ensure that this is the usual inner product for the algebras
-        over `R^n`::
-
-            sage: set_random_seed()
-            sage: J = HadamardEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            sage: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
-            True
-
-        """
-        return x.to_vector().inner_product(y.to_vector())
-
-
 class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
                       ConcreteEuclideanJordanAlgebra):
     r"""
@@ -2192,7 +2186,6 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
         True
     """
     def __init__(self, B, field=AA, **kwargs):
-        self._B = B
         n = B.nrows()
 
         if not B.is_positive_definite():
@@ -2213,13 +2206,24 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
                 z = V([z0] + zbar.list())
                 mult_table[i][j] = z
 
-        # The rank of this algebra is two, unless we're in a
-        # one-dimensional ambient space (because the rank is bounded
-        # by the ambient dimension).
+        # Inner products are real numbers and not algebra
+        # elements, so once we turn the algebra element
+        # into a vector in inner_product(), we never go
+        # back. As a result -- contrary to what we do with
+        # self._multiplication_table -- we store the inner
+        # product table as a plain old matrix and not as
+        # an algebra operator.
+        ip_table = B
+        self._inner_product_matrix = ip_table
+
         super(BilinearFormEJA, self).__init__(field,
                                               mult_table,
                                               check_axioms=False,
                                               **kwargs)
+
+        # The rank of this algebra is two, unless we're in a
+        # one-dimensional ambient space (because the rank is bounded
+        # by the ambient dimension).
         self.rank.set_cache(min(n,2))
 
         if n == 0:
@@ -2258,35 +2262,6 @@ class BilinearFormEJA(RationalBasisEuclideanJordanAlgebra,
 
         return cls(B, field, **kwargs)
 
-    def inner_product(self, x, y):
-        r"""
-        Half of the trace inner product.
-
-        This is defined so that the special case of the Jordan spin
-        algebra gets the usual inner product.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import BilinearFormEJA
-
-        TESTS:
-
-        Ensure that this is one-half of the trace inner-product when
-        the algebra isn't just the reals (when ``n`` isn't one). This
-        is in Faraut and Koranyi, and also my "On the symmetry..."
-        paper::
-
-            sage: set_random_seed()
-            sage: J = BilinearFormEJA.random_instance()
-            sage: n = J.dimension()
-            sage: x = J.random_element()
-            sage: y = J.random_element()
-            sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2)
-            True
-
-        """
-        return (self._B*x.to_vector()).inner_product(y.to_vector())
-
 
 class JordanSpinEJA(BilinearFormEJA):
     """
@@ -2332,9 +2307,9 @@ class JordanSpinEJA(BilinearFormEJA):
             sage: set_random_seed()
             sage: J = JordanSpinEJA.random_instance()
             sage: x,y = J.random_elements(2)
-            sage: X = x.natural_representation()
-            sage: Y = y.natural_representation()
-            sage: x.inner_product(y) == J.natural_inner_product(X,Y)
+            sage: actual = x.inner_product(y)
+            sage: expected = x.to_vector().inner_product(y.to_vector())
+            sage: actual == expected
             True
 
     """
@@ -2394,6 +2369,7 @@ class TrivialEJA(FiniteDimensionalEuclideanJordanAlgebra,
     """
     def __init__(self, field=AA, **kwargs):
         mult_table = []
+        self._inner_product_matrix = matrix(field,0)
         super(TrivialEJA, self).__init__(field,
                                          mult_table,
                                          check_axioms=False,
index 110b049573bbc8d8aeaaab768e9cebceadbf3c7e..6a9d10f65b7164627394c5ddb32bd682c323c5b2 100644 (file)
@@ -177,6 +177,10 @@ class FiniteDimensionalEuclideanJordanSubalgebra(FiniteDimensionalEuclideanJorda
 
         n = len(basis)
         mult_table = [[W.zero() for i in range(n)] for j in range(n)]
+        ip_table = [ [ self._superalgebra.inner_product(basis[i],basis[j])
+                       for i in range(n) ]
+                     for j in range(n) ]
+
         for i in range(n):
             for j in range(n):
                 product = basis[i]*basis[j]
@@ -187,6 +191,7 @@ class FiniteDimensionalEuclideanJordanSubalgebra(FiniteDimensionalEuclideanJorda
                 product_vector = V.from_vector(product.to_vector())
                 mult_table[i][j] = W.coordinate_vector(product_vector)
 
+        self._inner_product_matrix = matrix(field, ip_table)
         natural_basis = tuple( b.natural_representation() for b in basis )