From e529e0e2775cf50207c7d01d5907214d03cdff5c Mon Sep 17 00:00:00 2001
From: Michael Orlitzky <michael@orlitzky.com>
Date: Thu, 29 Aug 2019 22:48:51 -0400
Subject: [PATCH] eja: get a rudimentary spectral decomposition for operators
 working.

---
 mjo/eja/TODO            |  9 ---------
 mjo/eja/eja_operator.py | 40 ++++++++++++++++++++++++++++++++++++++++
 mjo/eja/eja_utils.py    | 33 +++++++--------------------------
 3 files changed, 47 insertions(+), 35 deletions(-)

diff --git a/mjo/eja/TODO b/mjo/eja/TODO
index b2495b5..985e50a 100644
--- a/mjo/eja/TODO
+++ b/mjo/eja/TODO
@@ -8,12 +8,3 @@
 
 5. Factor out the unit-norm basis (and operator symmetry) tests once
    all of the algebras pass.
-
-6. Implement spectral projector decomposition for EJA operators
-   using jordan_form() or eigenmatrix_right(). I suppose we can
-   ignore the problem of base rings for now and just let it crash
-   if we're not using AA as our base field.
-
-7. Do we really need to orthonormalize the basis in a subalgebra?
-   So long as we can decompose the operator (which is invariant
-   under changes of basis), who cares?
diff --git a/mjo/eja/eja_operator.py b/mjo/eja/eja_operator.py
index 0e898b5..c32ff1e 100644
--- 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)
diff --git a/mjo/eja/eja_utils.py b/mjo/eja/eja_utils.py
index d486b4c..cf75e32 100644
--- a/mjo/eja/eja_utils.py
+++ b/mjo/eja/eja_utils.py
@@ -26,11 +26,11 @@ def gram_schmidt(v):
         sage: u = gram_schmidt(v)
         sage: all( u_i.inner_product(u_i).sqrt() == 1 for u_i in u )
         True
-        sage: u[0].inner_product(u[1]) == 0
+        sage: bool(u[0].inner_product(u[1]) == 0)
         True
-        sage: u[0].inner_product(u[2]) == 0
+        sage: bool(u[0].inner_product(u[2]) == 0)
         True
-        sage: u[1].inner_product(u[2]) == 0
+        sage: bool(u[1].inner_product(u[2]) == 0)
         True
 
     TESTS:
@@ -67,29 +67,10 @@ def gram_schmidt(v):
     # And now drop all zero vectors again if they were "orthogonalized out."
     v = [ v_i for v_i in v if not v_i.is_zero() ]
 
-    # Now pretend to normalize, building a new ring R that contains
-    # all of the necessary square roots.
-    norms_squared = [0]*len(v)
-
-    for i in xrange(len(v)):
-        norms_squared[i] = v[i].inner_product(v[i])
-        ns = [norms_squared[i].numerator(), norms_squared[i].denominator()]
-
-        # Do the numerator and denominator separately so that we
-        # adjoin e.g. sqrt(2) and sqrt(3) instead of sqrt(2/3).
-        for j in xrange(len(ns)):
-            PR = PolynomialRing(R, 'z')
-            z = PR.gen()
-            p = z**2 - ns[j]
-            if p.is_irreducible():
-                R = NumberField(p,
-                                'sqrt' + str(ns[j]),
-                                embedding=RLF(ns[j]).sqrt())
-
-    # When we're done, we have to change every element's ring to the
-    # extension that we wound up with, and then normalize it (which
-    # should work, since "R" contains its norm now).
+    # Just normalize. If the algebra is missing the roots, we can't add
+    # them here because then our subalgebra would have a bigger field
+    # than the superalgebra.
     for i in xrange(len(v)):
-        v[i] = v[i].change_ring(R) / R(norms_squared[i]).sqrt()
+        v[i] = v[i] / v[i].norm()
 
     return v
-- 
2.43.2