eja: get a rudimentary spectral decomposition for operators working.

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

diff --git a/mjo/eja/eja_operator.py b/mjo/eja/eja_operator.py
index 4b72e5783210b51cd6d8d658f0f5d4e8c40a7f57..c32ff1ed7c2ba0aa87c842e5545f9bf204f43fac 100644 (file)
--- a/mjo/eja/eja_operator.py
+++ b/mjo/eja/eja_operator.py
@@ -13,6 +13,8 @@ class FiniteDimensionalEuclideanJordanAlgebraOperator(Map):
         F = domain_eja.base_ring()
         if not (F == codomain_eja.base_ring()):
             raise ValueError("domain and codomain must have the same base ring")
+        if not (F == mat.base_ring()):
+            raise ValueError("domain and matrix must have the same base ring")
 
         # We need to supply something here to avoid getting the
         # default Homset of the parent FiniteDimensionalAlgebra class,
@@ -88,10 +90,11 @@ class FiniteDimensionalEuclideanJordanAlgebraOperator(Map):
         different EJAs, that should blow up::
 
             sage: J1 = RealSymmetricEJA(2)
+            sage: id1 = identity_matrix(J1.base_ring(), 3)
             sage: J2 = JordanSpinEJA(3)
-            sage: id = identity_matrix(QQ, 3)
-            sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,J1,id)
-            sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,J2,id)
+            sage: id2 = identity_matrix(J2.base_ring(), 3)
+            sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,J1,id1)
+            sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,J2,id2)
             sage: f + g
             Traceback (most recent call last):
             ...
@@ -423,3 +426,43 @@ class FiniteDimensionalEuclideanJordanAlgebraOperator(Map):
         """
         # The matrix method returns a polynomial in 'x' but want one in 't'.
         return self.matrix().minimal_polynomial().change_variable_name('t')
+
+
+    def spectral_decomposition(self):
+        """
+        Return the spectral decomposition of this operator as a list of
+        (eigenvalue, orthogonal projector) pairs.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import RealSymmetricEJA
+
+        EXAMPLES::
+
+            sage: J = RealSymmetricEJA(4,AA)
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
+            sage: L0x = A(x).operator()
+            sage: Ps = [ P*l for (l,P) in L0x.spectral_decomposition() ]
+            sage: Ps[0] + Ps[1] == L0x
+            True
+
+        """
+        if not self.matrix().is_symmetric():
+            raise ValueError('algebra basis is not orthonormal')
+
+        D,P = self.matrix().jordan_form(subdivide=False,transformation=True)
+        eigenvalues = D.diagonal()
+        us = P.columns()
+        projectors = []
+        for i in range(len(us)):
+            # they won't be normalized, but they have to be
+            # for the spectral theorem to work.
+            us[i] = us[i]/us[i].norm()
+            mat = us[i].column()*us[i].row()
+            Pi = FiniteDimensionalEuclideanJordanAlgebraOperator(
+                   self.domain(),
+                   self.codomain(),
+                   mat)
+            projectors.append(Pi)
+        return zip(eigenvalues, projectors)