]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/euclidean_jordan_algebra.py
eja: add matrix inner products, and factor out the two used in R^n.
[sage.d.git] / mjo / eja / euclidean_jordan_algebra.py
index e2d644a903d8eca09a585064f96e31278197e4d6..a73cfaf689c7206f7644e10021006e4bcf6f08f7 100644 (file)
@@ -20,7 +20,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                               names='e',
                               assume_associative=False,
                               category=None,
-                              rank=None):
+                              rank=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:
@@ -43,7 +45,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                                  assume_associative=assume_associative,
                                  names=names,
                                  category=cat,
-                                 rank=rank)
+                                 rank=rank,
+                                 natural_basis=natural_basis,
+                                 inner_product=inner_product)
 
 
     def __init__(self, field,
@@ -51,7 +55,9 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
                  names='e',
                  assume_associative=False,
                  category=None,
-                 rank=None):
+                 rank=None,
+                 natural_basis=None,
+                 inner_product=None):
         """
         EXAMPLES:
 
@@ -66,6 +72,8 @@ 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,
@@ -80,6 +88,63 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
         fmt = "Euclidean Jordan algebra of degree {} over {}"
         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.
+
+        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.)
+
+        Note that this will always return a matrix. The standard basis
+        in `R^n` will be returned as `n`-by-`1` column matrices.
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricSimpleEJA(2)
+            sage: J.basis()
+            Family (e0, e1, e2)
+            sage: J.natural_basis()
+            (
+            [1 0]  [0 1]  [0 0]
+            [0 0], [1 0], [0 1]
+            )
+
+        ::
+
+            sage: J = JordanSpinSimpleEJA(2)
+            sage: J.basis()
+            Family (e0, e1)
+            sage: J.natural_basis()
+            (
+            [1]  [0]
+            [0], [1]
+            )
+
+        """
+        if self._natural_basis is None:
+            return tuple( b.vector().column() for b in self.basis() )
+        else:
+            return self._natural_basis
+
+
     def rank(self):
         """
         Return the rank of this EJA.
@@ -161,6 +226,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
@@ -191,7 +312,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()
@@ -281,7 +402,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
@@ -414,6 +535,92 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return self.span_of_powers().dimension()
 
 
+        def minimal_polynomial(self):
+            """
+            EXAMPLES::
+
+                sage: set_random_seed()
+                sage: x = random_eja().random_element()
+                sage: x.degree() == x.minimal_polynomial().degree()
+                True
+
+            ::
+
+                sage: set_random_seed()
+                sage: x = random_eja().random_element()
+                sage: x.degree() == x.minimal_polynomial().degree()
+                True
+
+            The minimal polynomial and the characteristic polynomial coincide
+            and are known (see Alizadeh, Example 11.11) for all elements of
+            the spin factor algebra that aren't scalar multiples of the
+            identity::
+
+                sage: set_random_seed()
+                sage: n = ZZ.random_element(2,10)
+                sage: J = JordanSpinSimpleEJA(n)
+                sage: y = J.random_element()
+                sage: while y == y.coefficient(0)*J.one():
+                ....:     y = J.random_element()
+                sage: y0 = y.vector()[0]
+                sage: y_bar = y.vector()[1:]
+                sage: actual = y.minimal_polynomial()
+                sage: x = SR.symbol('x', domain='real')
+                sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
+                sage: bool(actual == expected)
+                True
+
+            """
+            # The element we're going to call "minimal_polynomial()" on.
+            # Either myself, interpreted as an element of a finite-
+            # dimensional algebra, or an element of an associative
+            # subalgebra.
+            elt = None
+
+            if self.parent().is_associative():
+                elt = FiniteDimensionalAlgebraElement(self.parent(), self)
+            else:
+                V = self.span_of_powers()
+                assoc_subalg = self.subalgebra_generated_by()
+                # Mis-design warning: the basis used for span_of_powers()
+                # and subalgebra_generated_by() must be the same, and in
+                # the same order!
+                elt = assoc_subalg(V.coordinates(self.vector()))
+
+            # Recursive call, but should work since elt lives in an
+            # associative algebra.
+            return elt.minimal_polynomial()
+
+
+        def natural_representation(self):
+            """
+            Return a more-natural representation of this element.
+
+            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" representation of this element as a Hermitian
+            matrix, if it has one. If not, you get the usual representation.
+
+            EXAMPLES::
+
+                sage: J = ComplexHermitianSimpleEJA(3)
+                sage: J.one()
+                e0 + e5 + e8
+                sage: J.one().natural_representation()
+                [1 0 0 0 0 0]
+                [0 1 0 0 0 0]
+                [0 0 1 0 0 0]
+                [0 0 0 1 0 0]
+                [0 0 0 0 1 0]
+                [0 0 0 0 0 1]
+
+            """
+            B = self.parent().natural_basis()
+            W = B[0].matrix_space()
+            return W.linear_combination(zip(self.vector(), B))
+
+
         def operator_matrix(self):
             """
             Return the matrix that represents left- (or right-)
@@ -481,64 +688,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return fda_elt.matrix().transpose()
 
 
-
-        def minimal_polynomial(self):
-            """
-            EXAMPLES::
-
-                sage: set_random_seed()
-                sage: x = random_eja().random_element()
-                sage: x.degree() == x.minimal_polynomial().degree()
-                True
-
-            ::
-
-                sage: set_random_seed()
-                sage: x = random_eja().random_element()
-                sage: x.degree() == x.minimal_polynomial().degree()
-                True
-
-            The minimal polynomial and the characteristic polynomial coincide
-            and are known (see Alizadeh, Example 11.11) for all elements of
-            the spin factor algebra that aren't scalar multiples of the
-            identity::
-
-                sage: set_random_seed()
-                sage: n = ZZ.random_element(2,10)
-                sage: J = JordanSpinSimpleEJA(n)
-                sage: y = J.random_element()
-                sage: while y == y.coefficient(0)*J.one():
-                ....:     y = J.random_element()
-                sage: y0 = y.vector()[0]
-                sage: y_bar = y.vector()[1:]
-                sage: actual = y.minimal_polynomial()
-                sage: x = SR.symbol('x', domain='real')
-                sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
-                sage: bool(actual == expected)
-                True
-
-            """
-            # The element we're going to call "minimal_polynomial()" on.
-            # Either myself, interpreted as an element of a finite-
-            # dimensional algebra, or an element of an associative
-            # subalgebra.
-            elt = None
-
-            if self.parent().is_associative():
-                elt = FiniteDimensionalAlgebraElement(self.parent(), self)
-            else:
-                V = self.span_of_powers()
-                assoc_subalg = self.subalgebra_generated_by()
-                # Mis-design warning: the basis used for span_of_powers()
-                # and subalgebra_generated_by() must be the same, and in
-                # the same order!
-                elt = assoc_subalg(V.coordinates(self.vector()))
-
-            # Recursive call, but should work since elt lives in an
-            # associative algebra.
-            return elt.minimal_polynomial()
-
-
         def quadratic_representation(self, other=None):
             """
             Return the quadratic representation of this element.
@@ -609,7 +758,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()
@@ -772,7 +921,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()
 
@@ -810,7 +959,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)
 
 
 
@@ -865,7 +1017,7 @@ def _real_symmetric_basis(n, field=QQ):
                 # Beware, orthogonal but not normalized!
                 Sij = Eij + Eij.transpose()
             S.append(Sij)
-    return S
+    return tuple(S)
 
 
 def _complex_hermitian_basis(n, field=QQ):
@@ -902,7 +1054,7 @@ def _complex_hermitian_basis(n, field=QQ):
                 S.append(Sij_real)
                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
                 S.append(Sij_imag)
-    return S
+    return tuple(S)
 
 
 def _multiplication_table_from_matrix_basis(basis):
@@ -912,7 +1064,10 @@ def _multiplication_table_from_matrix_basis(basis):
     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.
+    elements) for an algebra of those matrices. A reordered copy
+    of the basis is also returned to work around the fact that
+    the ``span()`` in this function will change the order of the basis
+    from what we think it is, to... something else.
     """
     # In S^2, for example, we nominally have four coordinates even
     # though the space is of dimension three only. The vector space V
@@ -934,7 +1089,7 @@ def _multiplication_table_from_matrix_basis(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 = [ vec2mat(b) for b in W.basis() ]
+    S = tuple( vec2mat(b) for b in W.basis() )
 
     Qs = []
     for s in S:
@@ -952,7 +1107,7 @@ def _multiplication_table_from_matrix_basis(basis):
         Q = matrix(field, W.dimension(), Q_rows)
         Qs.append(Q)
 
-    return Qs
+    return (Qs, S)
 
 
 def _embed_complex_matrix(M):
@@ -979,7 +1134,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():
@@ -1007,9 +1162,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()
@@ -1021,14 +1176,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):
     """
@@ -1059,9 +1226,13 @@ def RealSymmetricSimpleEJA(n, field=QQ):
 
     """
     S = _real_symmetric_basis(n, field=field)
-    Qs = _multiplication_table_from_matrix_basis(S)
+    (Qs, T) = _multiplication_table_from_matrix_basis(S)
 
-    return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=n)
+    return FiniteDimensionalEuclideanJordanAlgebra(field,
+                                                   Qs,
+                                                   rank=n,
+                                                   natural_basis=T,
+                                                   inner_product=_matrix_ip)
 
 
 def ComplexHermitianSimpleEJA(n, field=QQ):
@@ -1083,8 +1254,22 @@ def ComplexHermitianSimpleEJA(n, field=QQ):
 
     """
     S = _complex_hermitian_basis(n)
-    Qs = _multiplication_table_from_matrix_basis(S)
-    return FiniteDimensionalEuclideanJordanAlgebra(field, Qs, rank=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,
+                                                   inner_product=ip)
 
 
 def QuaternionHermitianSimpleEJA(n):
@@ -1133,10 +1318,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 = []
@@ -1155,4 +1340,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)