]> gitweb.michael.orlitzky.com - sage.d.git/commitdiff
eja: refactor some matrix algebra stuff and break the tests.
authorMichael Orlitzky <michael@orlitzky.com>
Thu, 3 Dec 2020 19:24:44 +0000 (14:24 -0500)
committerMichael Orlitzky <michael@orlitzky.com>
Thu, 3 Dec 2020 19:24:44 +0000 (14:24 -0500)
mjo/eja/TODO
mjo/eja/eja_algebra.py

index 427a9539bb69aa4ab25e0687f0336a09377d57e0..cf870aedc40b1901ca81be9442ebf3d5b5f855cd 100644 (file)
@@ -15,10 +15,5 @@ RealSymmetricEJA(4):
 sage: F = J.base_ring()
 sage: a0 = (1/4)*X[4]**2*X[6]**2 - (1/2)*X[2]*X[5]*X[6]**2 - (1/2)*X[3]*X[4]*X[6]*X[7] + (F(2).sqrt()/2)*X[1]*X[5]*X[6]*X[7] + (1/4)*X[3]**2*X[7]**2 - (1/2)*X[0]*X[5]*X[7]**2 + (F(2).sqrt()/2)*X[2]*X[3]*X[6]*X[8] - (1/2)*X[1]*X[4]*X[6*X[8] - (1/2)*X[1]*X[3]*X[7]*X[8] + (F(2).sqrt()/2)*X[0]*X[4]*X[7]*X[8] + (1/4)*X[1]**2*X[8]**2 - (1/2)*X[0]*X[2]*X[8]**2 - (1/2)*X[2]*X[3]**2*X[9] + (F(2).sqrt()/2)*X[1]*X[3]*X[4]*X[9] - (1/2)*X[0]*X[4]**2*X[9] - (1/2)*X[1]**2*X[5]*X[9] + X[0]*X[2]*X[5]*X[9]
 
-5. Compute the scalar in the general natural_inner_product() for
-   matrices, so no overrides are necessary. Actually, this is
-   probably better implemented as a dimension_over_reals() method
-   that returns 1, 2, or 4.
-
-6. The main EJA element constructor is happy to convert between
+5. The main EJA element constructor is happy to convert between
    e.g. HadamardEJA(3) and JordanSpinEJA(3).
index e075ed21ab287fcefcfbf7c9674a39281b48b69c..5dee5d012d6c19d6fc5d79f804ad8aff2757c3d0 100644 (file)
@@ -1459,7 +1459,22 @@ class ConcreteEuclideanJordanAlgebra(RationalBasisEuclideanJordanAlgebra):
 
 class MatrixEuclideanJordanAlgebra:
     @staticmethod
-    def real_embed(M):
+    def dimension_over_reals():
+        r"""
+        The dimension of this matrix's base ring over the reals.
+
+        The reals are dimension one over themselves, obviously; that's
+        just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi`
+        have dimension two. Finally, the quaternions have dimension
+        four over the reals.
+
+        This is used to determine the size of the matrix returned from
+        :meth:`real_embed`, among other things.
+        """
+        raise NotImplementedError
+
+    @classmethod
+    def real_embed(cls,M):
         """
         Embed the matrix ``M`` into a space of real matrices.
 
@@ -1472,15 +1487,21 @@ class MatrixEuclideanJordanAlgebra:
           real_embed(M*N) = real_embed(M)*real_embed(N)
 
         """
-        raise NotImplementedError
+        if M.ncols() != M.nrows():
+            raise ValueError("the matrix 'M' must be square")
+        return M
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of :meth:`real_embed`.
         """
-        raise NotImplementedError
+        if M.ncols() != M.nrows():
+            raise ValueError("the matrix 'M' must be square")
+        if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero():
+            raise ValueError("the matrix 'M' must be a real embedding")
+        return M
 
     @staticmethod
     def jordan_product(X,Y):
@@ -1488,36 +1509,65 @@ class MatrixEuclideanJordanAlgebra:
 
     @classmethod
     def trace_inner_product(cls,X,Y):
+        r"""
+        Compute the trace inner-product of two real-embeddings.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
+            ....:                                  ComplexHermitianEJA,
+            ....:                                  QuaternionHermitianEJA)
+
+        EXAMPLES::
+
+        This gives the same answer as it would if we computed the trace
+        from the unembedded (original) matrices::
+
+            sage: set_random_seed()
+            sage: J = ComplexHermitianEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            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.trace_inner_product(Xe,Ye)
+            sage: actual == expected
+            True
+
+        ::
+
+            sage: set_random_seed()
+            sage: J = QuaternionHermitianEJA.random_instance()
+            sage: x,y = J.random_elements(2)
+            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.trace_inner_product(Xe,Ye)
+            sage: actual == expected
+            True
+
+        """
         Xu = cls.real_unembed(X)
         Yu = cls.real_unembed(Y)
         tr = (Xu*Yu).trace()
 
         try:
             # Works in QQ, AA, RDF, et cetera.
-            return tr.real()
+            return tr.real() / cls.dimension_over_reals()
         except AttributeError:
             # A quaternion doesn't have a real() method, but does
             # have coefficient_tuple() method that returns the
             # coefficients of 1, i, j, and k -- in that order.
-            return tr.coefficient_tuple()[0]
+            return tr.coefficient_tuple()[0] / cls.dimension_over_reals()
 
 
 class RealMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
-    def real_embed(M):
-        """
-        The identity function, for embedding real matrices into real
-        matrices.
-        """
-        return M
-
-    @staticmethod
-    def real_unembed(M):
-        """
-        The identity function, for unembedding real matrices from real
-        matrices.
-        """
-        return M
+    def dimension_over_reals():
+        return 1
 
 
 class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
@@ -1639,13 +1689,23 @@ class RealSymmetricEJA(ConcreteEuclideanJordanAlgebra,
                                                self.jordan_product,
                                                self.trace_inner_product,
                                                **kwargs)
+
+        # TODO: this could be factored out somehow, but is left here
+        # because the MatrixEuclideanJordanAlgebra is not presently
+        # a subclass of the FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        self.one.set_cache(self(matrix.identity(ZZ,n)))
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
+
 
 
 class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
-    def real_embed(M):
+    def dimension_over_reals():
+        return 2
+
+    @classmethod
+    def real_embed(cls,M):
         """
         Embed the n-by-n complex matrix ``M`` into the space of real
         matrices of size 2n-by-2n via the map the sends each entry `z = a +
@@ -1687,9 +1747,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(ComplexMatrixEuclideanJordanAlgebra,cls).real_embed(M)
         n = M.nrows()
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
 
         # We don't need any adjoined elements...
         field = M.base_ring().base_ring()
@@ -1703,8 +1762,8 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         return matrix.block(field, n, blocks)
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of _embed_complex_matrix().
 
@@ -1735,11 +1794,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(ComplexMatrixEuclideanJordanAlgebra,cls).real_unembed(M)
         n = ZZ(M.nrows())
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-        if not n.mod(2).is_zero():
-            raise ValueError("the matrix 'M' must be a complex embedding")
+        d = cls.dimension_over_reals()
 
         # If "M" was normalized, its base ring might have roots
         # adjoined and they can stick around after unembedding.
@@ -1757,9 +1814,9 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         # Go top-left to bottom-right (reading order), converting every
         # 2-by-2 block we see to a single complex element.
         elements = []
-        for k in range(n/2):
-            for j in range(n/2):
-                submat = M[2*k:2*k+2,2*j:2*j+2]
+        for k in range(n/d):
+            for j in range(n/d):
+                submat = M[d*k:d*k+d,d*j:d*j+d]
                 if submat[0,0] != submat[1,1]:
                     raise ValueError('bad on-diagonal submatrix')
                 if submat[0,1] != -submat[1,0]:
@@ -1767,38 +1824,7 @@ class ComplexMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
                 z = submat[0,0] + submat[0,1]*i
                 elements.append(z)
 
-        return matrix(F, n/2, elements)
-
-
-    @classmethod
-    def trace_inner_product(cls,X,Y):
-        """
-        Compute a matrix inner product in this algebra directly from
-        its real embedding.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import ComplexHermitianEJA
-
-        TESTS:
-
-        This gives the same answer as the slow, default method implemented
-        in :class:`MatrixEuclideanJordanAlgebra`::
-
-            sage: set_random_seed()
-            sage: J = ComplexHermitianEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            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.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/2
+        return matrix(F, n/d, elements)
 
 
 class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
@@ -1924,8 +1950,12 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
                                                   self.jordan_product,
                                                   self.trace_inner_product,
                                                   **kwargs)
+        # TODO: this could be factored out somehow, but is left here
+        # because the MatrixEuclideanJordanAlgebra is not presently
+        # a subclass of the FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        # TODO: pre-cache the identity!
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
 
     @staticmethod
     def _max_random_instance_size():
@@ -1941,7 +1971,11 @@ class ComplexHermitianEJA(ConcreteEuclideanJordanAlgebra,
 
 class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
     @staticmethod
-    def real_embed(M):
+    def dimension_over_reals():
+        return 4
+
+    @classmethod
+    def real_embed(cls,M):
         """
         Embed the n-by-n quaternion matrix ``M`` into the space of real
         matrices of size 4n-by-4n by first sending each quaternion entry `z
@@ -1980,10 +2014,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(QuaternionMatrixEuclideanJordanAlgebra,cls).real_embed(M)
         quaternions = M.base_ring()
         n = M.nrows()
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
 
         F = QuadraticField(-1, 'I')
         i = F.gen()
@@ -2006,8 +2039,8 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
 
 
 
-    @staticmethod
-    def real_unembed(M):
+    @classmethod
+    def real_unembed(cls,M):
         """
         The inverse of _embed_quaternion_matrix().
 
@@ -2037,11 +2070,9 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
             True
 
         """
+        super(QuaternionMatrixEuclideanJordanAlgebra,cls).real_unembed(M)
         n = ZZ(M.nrows())
-        if M.ncols() != n:
-            raise ValueError("the matrix 'M' must be square")
-        if not n.mod(4).is_zero():
-            raise ValueError("the matrix 'M' must be a quaternion embedding")
+        d = cls.dimension_over_reals()
 
         # Use the base ring of the matrix to ensure that its entries can be
         # multiplied by elements of the quaternion algebra.
@@ -2053,10 +2084,10 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
         # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
         # quaternion block.
         elements = []
-        for l in range(n/4):
-            for m in range(n/4):
+        for l in range(n/d):
+            for m in range(n/d):
                 submat = ComplexMatrixEuclideanJordanAlgebra.real_unembed(
-                    M[4*l:4*l+4,4*m:4*m+4] )
+                    M[d*l:d*l+d,d*m:d*m+d] )
                 if submat[0,0] != submat[1,1].conjugate():
                     raise ValueError('bad on-diagonal submatrix')
                 if submat[0,1] != -submat[1,0].conjugate():
@@ -2067,38 +2098,7 @@ class QuaternionMatrixEuclideanJordanAlgebra(MatrixEuclideanJordanAlgebra):
                 z += submat[0,1].imag()*k
                 elements.append(z)
 
-        return matrix(Q, n/4, elements)
-
-
-    @classmethod
-    def trace_inner_product(cls,X,Y):
-        """
-        Compute a matrix inner product in this algebra directly from
-        its real embedding.
-
-        SETUP::
-
-            sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA
-
-        TESTS:
-
-        This gives the same answer as the slow, default method implemented
-        in :class:`MatrixEuclideanJordanAlgebra`::
-
-            sage: set_random_seed()
-            sage: J = QuaternionHermitianEJA.random_instance()
-            sage: x,y = J.random_elements(2)
-            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.trace_inner_product(Xe,Ye)
-            sage: actual == expected
-            True
-
-        """
-        return RealMatrixEuclideanJordanAlgebra.trace_inner_product(X,Y)/4
+        return matrix(Q, n/d, elements)
 
 
 class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
@@ -2225,8 +2225,13 @@ class QuaternionHermitianEJA(ConcreteEuclideanJordanAlgebra,
                                                      self.jordan_product,
                                                      self.trace_inner_product,
                                                      **kwargs)
+        # TODO: this could be factored out somehow, but is left here
+        # because the MatrixEuclideanJordanAlgebra is not presently
+        # a subclass of the FDEJA class that defines rank() and one().
         self.rank.set_cache(n)
-        # TODO: cache one()!
+        idV = matrix.identity(ZZ, self.dimension_over_reals()*n)
+        self.one.set_cache(self(idV))
+
 
     @staticmethod
     def _max_random_instance_size():