eja: add the unique spectral_decomposition() for elements.

[sage.d.git] / mjo / eja / eja_element.py

diff --git a/mjo/eja/eja_element.py b/mjo/eja/eja_element.py
index 070501852ae4efa8d6845d74240d4495695f42a9..a4af4eaedbb4ce96c16aa31ab0e98c2fa4c5b6c7 100644 (file)
--- a/mjo/eja/eja_element.py
+++ b/mjo/eja/eja_element.py
@@ -1,3 +1,5 @@
+# -*- coding: utf-8 -*-
+

 from itertools import izip
 
 from sage.matrix.constructor import matrix
 from itertools import izip
 
 from sage.matrix.constructor import matrix


@@ -34,7 +36,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
         Return ``self`` raised to the power ``n``.
 
         Jordan algebras are always power-associative; see for
         Return ``self`` raised to the power ``n``.
 
         Jordan algebras are always power-associative; see for

-        example Faraut and Koranyi, Proposition II.1.2 (ii).
+        example Faraut and Korányi, Proposition II.1.2 (ii).

 
         We have to override this because our superclass uses row
         vectors instead of column vectors! We, on the other hand,
 
         We have to override this because our superclass uses row
         vectors instead of column vectors! We, on the other hand,


@@ -375,7 +377,7 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             True
 
         Ensure that the determinant is multiplicative on an associative
             True
 
         Ensure that the determinant is multiplicative on an associative

-        subalgebra as in Faraut and Koranyi's Proposition II.2.2::
+        subalgebra as in Faraut and Korányi's Proposition II.2.2::

 
             sage: set_random_seed()
             sage: J = random_eja().random_element().subalgebra_generated_by()
 
             sage: set_random_seed()
             sage: J = random_eja().random_element().subalgebra_generated_by()


@@ -405,7 +407,8 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
         SETUP::
 
 
         SETUP::
 

-            sage: from mjo.eja.eja_algebra import (JordanSpinEJA,
+            sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA,
+            ....:                                  JordanSpinEJA,

             ....:                                  random_eja)
 
         EXAMPLES:
             ....:                                  random_eja)
 
         EXAMPLES:


@@ -460,6 +463,33 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             ...
             ValueError: element is not invertible
 
             ...
             ValueError: element is not invertible
 

+        Proposition II.2.3 in Faraut and Korányi says that the inverse
+        of an element is the inverse of its left-multiplication operator
+        applied to the algebra's identity, when that inverse exists::
+
+            sage: set_random_seed()
+            sage: J = random_eja()
+            sage: x = J.random_element()
+            sage: (not x.operator().is_invertible()) or (
+            ....:    x.operator().inverse()(J.one()) == x.inverse() )
+            True
+
+        Proposition II.2.4 in Faraut and Korányi gives a formula for
+        the inverse based on the characteristic polynomial and the
+        Cayley-Hamilton theorem for Euclidean Jordan algebras::
+
+            sage: set_random_seed()
+            sage: J = ComplexHermitianEJA(3)
+            sage: x = J.random_element()
+            sage: while not x.is_invertible():
+            ....:     x = J.random_element()
+            sage: r = J.rank()
+            sage: a = x.characteristic_polynomial().coefficients(sparse=False)
+            sage: expected  = (-1)^(r+1)/x.det()
+            sage: expected *= sum( a[i+1]*x^i for i in range(r) )
+            sage: x.inverse() == expected
+            True
+

         """
         if not self.is_invertible():
             raise ValueError("element is not invertible")
         """
         if not self.is_invertible():
             raise ValueError("element is not invertible")


@@ -1010,12 +1040,82 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
 
 
 
 
 
 

+    def spectral_decomposition(self):
+        """
+        Return the unique spectral decomposition of this element.
+
+        ALGORITHM:
+
+        Following Faraut and Korányi's Theorem III.1.1, we restrict this
+        element's left-multiplication-by operator to the subalgebra it
+        generates. We then compute the spectral decomposition of that
+        operator, and the spectral projectors we get back must be the
+        left-multiplication-by operators for the idempotents we
+        seek. Thus applying them to the identity element gives us those
+        idempotents.
+
+        Since the eigenvalues are required to be distinct, we take
+        the spectral decomposition of the zero element to be zero
+        times the identity element of the algebra (which is idempotent,
+        obviously).
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealSymmetricEJA
+
+        EXAMPLES:
+
+        The spectral decomposition of the identity is ``1`` times itself,
+        and the spectral decomposition of zero is ``0`` times the identity::
+
+            sage: J = RealSymmetricEJA(3,AA)
+            sage: J.one()
+            e0 + e2 + e5
+            sage: J.one().spectral_decomposition()
+            [(1, e0 + e2 + e5)]
+            sage: J.zero().spectral_decomposition()
+            [(0, e0 + e2 + e5)]
+
+        TESTS::
+
+            sage: J = RealSymmetricEJA(4,AA)
+            sage: x = sum(J.gens())
+            sage: sd = x.spectral_decomposition()
+            sage: l0 = sd[0][0]
+            sage: l1 = sd[1][0]
+            sage: c0 = sd[0][1]
+            sage: c1 = sd[1][1]
+            sage: c0.inner_product(c1) == 0
+            True
+            sage: c0.is_idempotent()
+            True
+            sage: c1.is_idempotent()
+            True
+            sage: c0 + c1 == J.one()
+            True
+            sage: l0*c0 + l1*c1 == x
+            True

 
 

-    def subalgebra_generated_by(self):
+        """
+        P = self.parent()
+        A = self.subalgebra_generated_by(orthonormalize_basis=True)
+        result = []
+        for (evalue, proj) in A(self).operator().spectral_decomposition():
+            result.append( (evalue, proj(A.one()).superalgebra_element()) )
+        return result
+
+    def subalgebra_generated_by(self, orthonormalize_basis=False):

         """
         Return the associative subalgebra of the parent EJA generated
         by this element.
 
         """
         Return the associative subalgebra of the parent EJA generated
         by this element.
 

+        Since our parent algebra is unital, we want "subalgebra" to mean
+        "unital subalgebra" as well; thus the subalgebra that an element
+        generates will itself be a Euclidean Jordan algebra after
+        restricting the algebra operations appropriately. This is the
+        subalgebra that Faraut and Korányi work with in section II.2, for
+        example.
+

         SETUP::
 
             sage: from mjo.eja.eja_algebra import random_eja
         SETUP::
 
             sage: from mjo.eja.eja_algebra import random_eja


@@ -1040,17 +1140,16 @@ class FiniteDimensionalEuclideanJordanAlgebraElement(IndexedFreeModuleElement):
             sage: A(x^2) == A(x)*A(x)
             True
 
             sage: A(x^2) == A(x)*A(x)
             True
 

-        The subalgebra generated by the zero element is trivial::
+        By definition, the subalgebra generated by the zero element is the
+        one-dimensional algebra generated by the identity element::

 
             sage: set_random_seed()
             sage: A = random_eja().zero().subalgebra_generated_by()
 
             sage: set_random_seed()
             sage: A = random_eja().zero().subalgebra_generated_by()

-            sage: A
-            Euclidean Jordan algebra of dimension 0 over...
-            sage: A.one()
-            0
+            sage: A.dimension()
+            1

 
         """
 
         """

-        return FiniteDimensionalEuclideanJordanElementSubalgebra(self)
+        return FiniteDimensionalEuclideanJordanElementSubalgebra(self, orthonormalize_basis)

 
 
     def subalgebra_idempotent(self):
 
 
     def subalgebra_idempotent(self):