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 0e898b59fb9d8a99090da6f4baec3e49df51d9cd..c32ff1ed7c2ba0aa87c842e5545f9bf204f43fac 100644 (file)
--- a/mjo/eja/eja_operator.py
+++ b/mjo/eja/eja_operator.py
@@ -426,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)