]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: support conversion of naturally-represented elements into an EJA.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index 8d9b27e974fff75f769579853d7ac60716d26298..b25839815289ba6f60630a14b790dbfe29f8969b 100644 (file)
@@ -21,7 +21,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                               assume_associative=False,
                               category=None,
                               rank=None,
-                              natural_basis=None):
+                              natural_basis=None,
+                              inner_product=None):
         n = len(mult_table)
         mult_table = [b.base_extend(field) for b in mult_table]
         for b in mult_table:
@@ -45,7 +46,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                  names=names,
                                  category=cat,
                                  rank=rank,
-                                 natural_basis=natural_basis)
+                                 natural_basis=natural_basis,
+                                 inner_product=inner_product)
 
 
     def __init__(self, field,
@@ -54,7 +56,8 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                  assume_associative=False,
                  category=None,
                  rank=None,
-                 natural_basis=None):
+                 natural_basis=None,
+                 inner_product=None):
         """
         EXAMPLES:
 
@@ -70,6 +73,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         """
         self._rank = rank
         self._natural_basis = natural_basis
+        self._inner_product = inner_product
         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
         fda.__init__(field,
                      mult_table,
@@ -85,6 +89,20 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         return fmt.format(self.degree(), self.base_ring())
 
 
+    def inner_product(self, x, y):
+        """
+        The inner product associated with this Euclidean Jordan algebra.
+
+        Will default to the trace inner product if nothing else.
+        """
+        if (not x in self) or (not y in self):
+            raise TypeError("arguments must live in this algebra")
+        if self._inner_product is None:
+            return x.trace_inner_product(y)
+        else:
+            return self._inner_product(x,y)
+
+
     def natural_basis(self):
         """
         Return a more-natural representation of this algebra's basis.
@@ -142,6 +160,51 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         An element of a Euclidean Jordan algebra.
         """
 
+        def __init__(self, A, elt=None):
+            """
+            EXAMPLES:
+
+            The identity in `S^n` is converted to the identity in the EJA::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: I = identity_matrix(QQ,3)
+                sage: J(I) == J.one()
+                True
+
+            This skew-symmetric matrix can't be represented in the EJA::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: A = matrix(QQ,3, lambda i,j: i-j)
+                sage: J(A)
+                Traceback (most recent call last):
+                ...
+                ArithmeticError: vector is not in free module
+
+            """
+            # Goal: if we're given a matrix, and if it lives in our
+            # parent algebra's "natural ambient space," convert it
+            # into an algebra element.
+            #
+            # The catch is, we make a recursive call after converting
+            # the given matrix into a vector that lives in the algebra.
+            # This we need to try the parent class initializer first,
+            # to avoid recursing forever if we're given something that
+            # already fits into the algebra, but also happens to live
+            # in the parent's "natural ambient space" (this happens with
+            # vectors in R^n).
+            try:
+                FiniteDimensionalAlgebraElement.__init__(self, A, elt)
+            except ValueError:
+                natural_basis = A.natural_basis()
+                if elt in natural_basis[0].matrix_space():
+                    # Thanks for nothing! Matrix spaces aren't vector
+                    # spaces in Sage, so we have to figure out its
+                    # natural-basis coordinates ourselves.
+                    V = VectorSpace(elt.base_ring(), elt.nrows()**2)
+                    W = V.span( _mat2vec(s) for s in natural_basis )
+                    coords =  W.coordinates(_mat2vec(elt))
+                    FiniteDimensionalAlgebraElement.__init__(self, A, coords)
+
         def __pow__(self, n):
             """
             Return ``self`` raised to the power ``n``.
@@ -208,6 +271,62 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                 raise NotImplementedError('irregular element')
 
 
+        def inner_product(self, other):
+            """
+            Return the parent algebra's inner product of myself and ``other``.
+
+            EXAMPLES:
+
+            The inner product in the Jordan spin algebra is the usual
+            inner product on `R^n` (this example only works because the
+            basis for the Jordan algebra is the standard basis in `R^n`)::
+
+                sage: J = JordanSpinSimpleEJA(3)
+                sage: x = vector(QQ,[1,2,3])
+                sage: y = vector(QQ,[4,5,6])
+                sage: x.inner_product(y)
+                32
+                sage: J(x).inner_product(J(y))
+                32
+
+            The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
+            multiplication is the usual matrix multiplication in `S^n`,
+            so the inner product of the identity matrix with itself
+            should be the `n`::
+
+                sage: J = RealSymmetricSimpleEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
+            Likewise, the inner product on `C^n` is `<X,Y> =
+            Re(trace(X*Y))`, where we must necessarily take the real
+            part because the product of Hermitian matrices may not be
+            Hermitian::
+
+                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J.one().inner_product(J.one())
+                3
+
+            TESTS:
+
+            Ensure that we can always compute an inner product, and that
+            it gives us back a real number::
+
+                sage: set_random_seed()
+                sage: J = random_eja()
+                sage: x = J.random_element()
+                sage: y = J.random_element()
+                sage: x.inner_product(y) in RR
+                True
+
+            """
+            P = self.parent()
+            if not other in P:
+                raise TypeError("'other' must live in the same algebra")
+
+            return P.inner_product(self, other)
+
+
         def operator_commutes_with(self, other):
             """
             Return whether or not this element operator-commutes
@@ -238,7 +357,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
 
             """
             if not other in self.parent():
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             A = self.operator_matrix()
             B = other.operator_matrix()
@@ -328,7 +447,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             # TODO: we can do better once the call to is_invertible()
             # doesn't crash on irregular elements.
             #if not self.is_invertible():
-            #    raise ArgumentError('element is not invertible')
+            #    raise ValueError('element is not invertible')
 
             # We do this a little different than the usual recursive
             # call to a finite-dimensional algebra element, because we
@@ -684,7 +803,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             if other is None:
                 other=self
             elif not other in self.parent():
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             L = self.operator_matrix()
             M = other.operator_matrix()
@@ -847,7 +966,7 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             Return the trace inner product of myself and ``other``.
             """
             if not other in self.parent():
-                raise ArgumentError("'other' must live in the same algebra")
+                raise TypeError("'other' must live in the same algebra")
 
             return (self*other).trace()
 
@@ -885,7 +1004,10 @@ def eja_rn(dimension, field=QQ):
     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
            for i in xrange(dimension) ]
 
-    return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
+    return FiniteDimensionalEuclideanJordanAlgebra(field,
+                                                   Qs,
+                                                   rank=dimension,
+                                                   inner_product=_usual_ip)
 
 
 
@@ -980,6 +1102,12 @@ def _complex_hermitian_basis(n, field=QQ):
     return tuple(S)
 
 
+def _mat2vec(m):
+        return vector(m.base_ring(), m.list())
+
+def _vec2mat(v):
+        return matrix(v.base_ring(), sqrt(v.degree()), v.list())
+
 def _multiplication_table_from_matrix_basis(basis):
     """
     At least three of the five simple Euclidean Jordan algebras have the
@@ -1000,19 +1128,13 @@ def _multiplication_table_from_matrix_basis(basis):
     field = basis[0].base_ring()
     dimension = basis[0].nrows()
 
-    def mat2vec(m):
-        return vector(field, m.list())
-
-    def vec2mat(v):
-        return matrix(field, dimension, v.list())
-
     V = VectorSpace(field, dimension**2)
-    W = V.span( mat2vec(s) for s in basis )
+    W = V.span( _mat2vec(s) for s in basis )
 
     # Taking the span above reorders our basis (thanks, jerk!) so we
     # need to put our "matrix basis" in the same order as the
     # (reordered) vector basis.
-    S = tuple( vec2mat(b) for b in W.basis() )
+    S = tuple( _vec2mat(b) for b in W.basis() )
 
     Qs = []
     for s in S:
@@ -1025,7 +1147,7 @@ def _multiplication_table_from_matrix_basis(basis):
         # why we're computing rows here and not columns.
         Q_rows = []
         for t in S:
-            this_row = mat2vec((s*t + t*s)/2)
+            this_row = _mat2vec((s*t + t*s)/2)
             Q_rows.append(W.coordinates(this_row))
         Q = matrix(field, W.dimension(), Q_rows)
         Qs.append(Q)
@@ -1057,7 +1179,7 @@ def _embed_complex_matrix(M):
     """
     n = M.nrows()
     if M.ncols() != n:
-        raise ArgumentError("the matrix 'M' must be square")
+        raise ValueError("the matrix 'M' must be square")
     field = M.base_ring()
     blocks = []
     for z in M.list():
@@ -1085,9 +1207,9 @@ def _unembed_complex_matrix(M):
     """
     n = ZZ(M.nrows())
     if M.ncols() != n:
-        raise ArgumentError("the matrix 'M' must be square")
+        raise ValueError("the matrix 'M' must be square")
     if not n.mod(2).is_zero():
-        raise ArgumentError("the matrix 'M' must be a complex embedding")
+        raise ValueError("the matrix 'M' must be a complex embedding")
 
     F = QuadraticField(-1, 'i')
     i = F.gen()
@@ -1099,14 +1221,26 @@ def _unembed_complex_matrix(M):
         for j in xrange(n/2):
             submat = M[2*k:2*k+2,2*j:2*j+2]
             if submat[0,0] != submat[1,1]:
-                raise ArgumentError('bad real submatrix')
+                raise ValueError('bad real submatrix')
             if submat[0,1] != -submat[1,0]:
-                raise ArgumentError('bad imag submatrix')
+                raise ValueError('bad imag submatrix')
             z = submat[0,0] + submat[1,0]*i
             elements.append(z)
 
     return matrix(F, n/2, elements)
 
+# The usual inner product on R^n.
+def _usual_ip(x,y):
+    return x.vector().inner_product(y.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.
+def _matrix_ip(X,Y):
+    X_mat = X.natural_representation()
+    Y_mat = Y.natural_representation()
+    return (X_mat*Y_mat).trace()
+
 
 def RealSymmetricSimpleEJA(n, field=QQ):
     """
@@ -1142,7 +1276,8 @@ def RealSymmetricSimpleEJA(n, field=QQ):
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,
-                                                   natural_basis=T)
+                                                   natural_basis=T,
+                                                   inner_product=_matrix_ip)
 
 
 def ComplexHermitianSimpleEJA(n, field=QQ):
@@ -1165,10 +1300,21 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
     """
     S = _complex_hermitian_basis(n)
     (Qs, T) = _multiplication_table_from_matrix_basis(S)
+
+    # Since a+bi on the diagonal is represented as
+    #
+    #   a + bi  = [  a  b  ]
+    #             [ -b  a  ],
+    #
+    # we'll double-count the "a" entries if we take the trace of
+    # the embedding.
+    ip = lambda X,Y: _matrix_ip(X,Y)/2
+
     return FiniteDimensionalEuclideanJordanAlgebra(field,
                                                    Qs,
                                                    rank=n,
-                                                   natural_basis=T)
+                                                   natural_basis=T,
+                                                   inner_product=ip)
 
 
 def QuaternionHermitianSimpleEJA(n):
@@ -1217,10 +1363,10 @@ def JordanSpinSimpleEJA(n, field=QQ):
 
     In one dimension, this is the reals under multiplication::
 
-      sage: J1 = JordanSpinSimpleEJA(1)
-      sage: J2 = eja_rn(1)
-      sage: J1 == J2
-      True
+        sage: J1 = JordanSpinSimpleEJA(1)
+        sage: J2 = eja_rn(1)
+        sage: J1 == J2
+        True
 
     """
     Qs = []
@@ -1239,4 +1385,7 @@ def JordanSpinSimpleEJA(n, field=QQ):
     # The rank of the spin factor algebra is two, UNLESS we're in a
     # one-dimensional ambient space (the rank is bounded by the
     # ambient dimension).
-    return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=min(n,2))
+    return FiniteDimensionalEuclideanJordanAlgebra(field,
+                                                   Qs,
+                                                   rank=min(n,2),
+                                                   inner_product=_usual_ip)