]> 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 c0b7787a535b7754c446fcb0046c300b21695579..b25839815289ba6f60630a14b790dbfe29f8969b 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.
@@ -95,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``.
@@ -161,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
@@ -191,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()
@@ -281,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
@@ -414,7 +580,93 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             return self.span_of_powers().dimension()
 
 
-        def matrix(self):
+        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-)
             multiplication by this element in the parent algebra.
@@ -480,71 +732,6 @@ class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
             return fda_elt.matrix().transpose()
 
-        #
-        # The plan is to eventually phase out "matrix()", which sounds
-        # too much like "matrix_representation()", in favor of the more-
-        # accurate "operator_matrix()". But we need to override matrix()
-        # to keep parent class methods happy in the meantime.
-        #
-        operator_matrix = matrix
-
-
-        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):
             """
@@ -616,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()
@@ -779,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()
 
@@ -817,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)
 
 
 
@@ -872,7 +1062,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):
@@ -909,9 +1099,15 @@ 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 _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
@@ -919,7 +1115,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
@@ -929,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 = [ vec2mat(b) for b in W.basis() ]
+    S = tuple( _vec2mat(b) for b in W.basis() )
 
     Qs = []
     for s in S:
@@ -954,12 +1147,12 @@ 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)
 
-    return Qs
+    return (Qs, S)
 
 
 def _embed_complex_matrix(M):
@@ -986,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():
@@ -1014,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()
@@ -1028,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):
     """
@@ -1066,9 +1271,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):
@@ -1090,8 +1299,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):
@@ -1140,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 = []
@@ -1162,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)