]> gitweb.michael.orlitzky.com - sage.d.git/commitdiff
eja: refactor some of the basis and inner-product stuff.
authorMichael Orlitzky <michael@orlitzky.com>
Wed, 21 Aug 2019 20:10:14 +0000 (16:10 -0400)
committerMichael Orlitzky <michael@orlitzky.com>
Wed, 21 Aug 2019 20:10:14 +0000 (16:10 -0400)
This is movement towards eventually cheating on the charpoly
coefficients, which we should be able to compute in the "nice" basis
and then scale to the normalized one. The coefficients are polynomials
in "the coordinates of x", and those coordinates change only by a
scalar multiple when we normalize the basis.

mjo/eja/eja_algebra.py

index a7e56187550111a928823d536dc3ff88347b7424..d0e7b074ce41a00f66d6f6dcd487653f3a8b1674 100644 (file)
@@ -416,9 +416,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(CombinatorialFreeModule):
             True
 
         """
-        if (not x in self) or (not y in self):
-            raise TypeError("arguments must live in this algebra")
-        return x.trace_inner_product(y)
+        X = x.natural_representation()
+        Y = y.natural_representation()
+        return self.__class__.natural_inner_product(X,Y)
 
 
     def is_trivial(self):
@@ -537,6 +537,20 @@ 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):
         """
@@ -748,7 +762,30 @@ class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field, mult_table, rank=n, **kwargs)
 
     def inner_product(self, x, y):
-        return _usual_ip(x,y)
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealCartesianProductEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: J = RealCartesianProductEJA(n)
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())
 
 
 def random_eja():
@@ -802,7 +839,7 @@ def random_eja():
 
 
 
-def _real_symmetric_basis(n, field, normalize):
+def _real_symmetric_basis(n, field):
     """
     Return a basis for the space of real symmetric n-by-n matrices.
 
@@ -814,7 +851,7 @@ def _real_symmetric_basis(n, field, normalize):
 
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
-        sage: B = _real_symmetric_basis(n, QQbar, False)
+        sage: B = _real_symmetric_basis(n, QQ)
         sage: all( M.is_symmetric() for M in  B)
         True
 
@@ -829,13 +866,11 @@ def _real_symmetric_basis(n, field, normalize):
                 Sij = Eij
             else:
                 Sij = Eij + Eij.transpose()
-            if normalize:
-                Sij = Sij / _real_symmetric_matrix_ip(Sij,Sij).sqrt()
             S.append(Sij)
     return tuple(S)
 
 
-def _complex_hermitian_basis(n, field, normalize):
+def _complex_hermitian_basis(n, field):
     """
     Returns a basis for the space of complex Hermitian n-by-n matrices.
 
@@ -854,7 +889,7 @@ def _complex_hermitian_basis(n, field, normalize):
         sage: set_random_seed()
         sage: n = ZZ.random_element(1,5)
         sage: field = QuadraticField(2, 'sqrt2')
-        sage: B = _complex_hermitian_basis(n, field, False)
+        sage: B = _complex_hermitian_basis(n, field)
         sage: all( M.is_symmetric() for M in  B)
         True
 
@@ -877,8 +912,7 @@ def _complex_hermitian_basis(n, field, normalize):
                 Sij = _embed_complex_matrix(Eij)
                 S.append(Sij)
             else:
-                # Beware, orthogonal but not normalized! The second one
-                # has a minus because it's conjugated.
+                # The second one has a minus because it's conjugated.
                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
                 S.append(Sij_real)
                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
@@ -886,11 +920,7 @@ def _complex_hermitian_basis(n, field, normalize):
 
     # Since we embedded these, we can drop back to the "field" that we
     # started with instead of the complex extension "F".
-    S = [ s.change_ring(field) for s in S ]
-    if normalize:
-        S = [ s / _complex_hermitian_matrix_ip(s,s).sqrt() for s in S ]
-
-    return tuple(S)
+    return tuple( s.change_ring(field) for s in S )
 
 
 
@@ -1206,10 +1236,6 @@ def _unembed_quaternion_matrix(M):
     return matrix(Q, n/4, elements)
 
 
-# The usual inner product on R^n.
-def _usual_ip(x,y):
-    return x.to_vector().inner_product(y.to_vector())
-
 # The inner product used for the real symmetric simple EJA.
 # We keep it as a separate function because e.g. the complex
 # algebra uses the same inner product, except divided by 2.
@@ -1218,15 +1244,6 @@ def _matrix_ip(X,Y):
     Y_mat = Y.natural_representation()
     return (X_mat*Y_mat).trace()
 
-def _real_symmetric_matrix_ip(X,Y):
-    return (X*Y).trace()
-
-def _complex_hermitian_matrix_ip(X,Y):
-    # This takes EMBEDDED matrices.
-    Xu = _unembed_complex_matrix(X)
-    Yu = _unembed_complex_matrix(Y)
-    # The trace need not be real; consider Xu = (i*I) and Yu = I.
-    return ((Xu*Yu).trace()).vector()[0] # real part, I guess
 
 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
@@ -1310,6 +1327,8 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     """
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
+        S = _real_symmetric_basis(n, field)
+
         if n > 1 and normalize_basis:
             # We'll need sqrt(2) to normalize the basis, and this
             # winds up in the multiplication table, so the whole
@@ -1319,8 +1338,12 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
             p = z**2 - 2
             if p.is_irreducible():
                 field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
+            S = [ s.change_ring(field) for s in S ]
+            self._basis_denormalizers = tuple(
+                self.__class__.natural_inner_product(s,s).sqrt()
+                for s in S )
+            S = tuple( s/c for (s,c) in zip(S,self._basis_denormalizers) )
 
-        S = _real_symmetric_basis(n, field, normalize_basis)
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(RealSymmetricEJA, self)
@@ -1330,10 +1353,6 @@ class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               natural_basis=S,
                               **kwargs)
 
-    def inner_product(self, x, y):
-        X = x.natural_representation()
-        Y = y.natural_representation()
-        return _real_symmetric_matrix_ip(X,Y)
 
 
 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
@@ -1408,6 +1427,8 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
 
     """
     def __init__(self, n, field=QQ, normalize_basis=True, **kwargs):
+        S = _complex_hermitian_basis(n, field)
+
         if n > 1 and normalize_basis:
             # We'll need sqrt(2) to normalize the basis, and this
             # winds up in the multiplication table, so the whole
@@ -1417,8 +1438,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
             p = z**2 - 2
             if p.is_irreducible():
                 field = NumberField(p, 'sqrt2', embedding=RLF(2).sqrt())
+            S = [ s.change_ring(field) for s in S ]
+            self._basis_denormalizers = tuple(
+                self.__class__.natural_inner_product(s,s).sqrt()
+                for s in S )
+            S = tuple( s/c for (s,c) in zip(S,self._basis_denormalizers) )
 
-        S = _complex_hermitian_basis(n, field, normalize_basis)
         Qs = _multiplication_table_from_matrix_basis(S)
 
         fdeja = super(ComplexHermitianEJA, self)
@@ -1429,11 +1454,12 @@ class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
                               **kwargs)
 
 
-    def inner_product(self, x, y):
-        X = x.natural_representation()
-        Y = y.natural_representation()
-        return _complex_hermitian_matrix_ip(X,Y)
-
+    @staticmethod
+    def natural_inner_product(X,Y):
+        Xu = _unembed_complex_matrix(X)
+        Yu = _unembed_complex_matrix(Y)
+        # The trace need not be real; consider Xu = (i*I) and Yu = I.
+        return ((Xu*Yu).trace()).vector()[0] # real part, I guess
 
 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
     """
@@ -1586,4 +1612,27 @@ class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
         return fdeja.__init__(field, mult_table, rank=min(n,2), **kwargs)
 
     def inner_product(self, x, y):
-        return _usual_ip(x,y)
+        """
+        Faster to reimplement than to use natural representations.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import JordanSpinEJA
+
+        TESTS:
+
+        Ensure that this is the usual inner product for the algebras
+        over `R^n`::
+
+            sage: set_random_seed()
+            sage: n = ZZ.random_element(1,5)
+            sage: J = JordanSpinEJA(n)
+            sage: x = J.random_element()
+            sage: y = J.random_element()
+            sage: X = x.natural_representation()
+            sage: Y = y.natural_representation()
+            sage: x.inner_product(y) == J.__class__.natural_inner_product(X,Y)
+            True
+
+        """
+        return x.to_vector().inner_product(y.to_vector())