mjo/eja/eja_element_subalgebra.py

   1 from sage.matrix.constructor import matrix
   2 from sage.misc.cachefunc import cached_method
   3 from sage.rings.all import QQ
   4
   5 from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
   6
   7
   8 class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclideanJordanSubalgebra):
   9     def __init__(self, elt, orthonormalize_basis):
  10         self._superalgebra = elt.parent()
  11         category = self._superalgebra.category().Associative()
  12         V = self._superalgebra.vector_space()
  13         field = self._superalgebra.base_ring()
  14
  15         # This list is guaranteed to contain all independent powers,
  16         # because it's the maximal set of powers that could possibly
  17         # be independent (by a dimension argument).
  18         powers = [ elt**k for k in range(V.dimension()) ]
  19         power_vectors = [ p.to_vector() for p in powers ]
  20         P = matrix(field, power_vectors)
  21
  22         if orthonormalize_basis == False:
  23             # Echelonize the matrix ourselves, because otherwise the
  24             # call to P.pivot_rows() below can choose a non-optimal
  25             # row-reduction algorithm. In particular, scaling can
  26             # help over AA because it avoids the RecursionError that
  27             # gets thrown when we have to look too hard for a root.
  28             #
  29             # Beware: QQ supports an entirely different set of "algorithm"
  30             # keywords than do AA and RR.
  31             algo = None
  32             if field is not QQ:
  33                 algo = "scaled_partial_pivoting"
  34             P.echelonize(algorithm=algo)
  35
  36             # In this case, we just need to figure out which elements
  37             # of the "powers" list are redundant... First compute the
  38             # vector subspace spanned by the powers of the given
  39             # element.
  40
  41             # Figure out which powers form a linearly-independent set.
  42             ind_rows = P.pivot_rows()
  43
  44             # Pick those out of the list of all powers.
  45             superalgebra_basis = tuple(map(powers.__getitem__, ind_rows))
  46         else:
  47             # If we're going to orthonormalize the basis anyway, we
  48             # might as well just do Gram-Schmidt on the whole list of
  49             # powers. The redundant ones will get zero'd out. If this
  50             # looks like a roundabout way to orthonormalize, it is.
  51             # But converting everything from algebra elements to vectors
  52             # to matrices and then back again turns out to be about
  53             # as fast as reimplementing our own Gram-Schmidt that
  54             # works in an EJA.
  55             G,_ = P.gram_schmidt(orthonormal=True)
  56             basis_vectors = [ g for g in G.rows() if not g.is_zero() ]
  57             superalgebra_basis = [ self._superalgebra.from_vector(b)
  58                                    for b in basis_vectors ]
  59
  60         fdeja = super(FiniteDimensionalEuclideanJordanElementSubalgebra, self)
  61         fdeja.__init__(self._superalgebra,
  62                        superalgebra_basis,
  63                        category=category,
  64                        check_axioms=False)
  65
  66         # The rank is the highest possible degree of a minimal
  67         # polynomial, and is bounded above by the dimension. We know
  68         # in this case that there's an element whose minimal
  69         # polynomial has the same degree as the space's dimension
  70         # (remember how we constructed the space?), so that must be
  71         # its rank too.
  72         self.rank.set_cache(self.dimension())
  73
  74
  75     @cached_method
  76     def one(self):
  77         """
  78         Return the multiplicative identity element of this algebra.
  79
  80         The superclass method computes the identity element, which is
  81         beyond overkill in this case: the superalgebra identity
  82         restricted to this algebra is its identity. Note that we can't
  83         count on the first basis element being the identity -- it might
  84         have been scaled if we orthonormalized the basis.
  85
  86         SETUP::
  87
  88             sage: from mjo.eja.eja_algebra import (HadamardEJA,
  89             ....:                                  random_eja)
  90
  91         EXAMPLES::
  92
  93             sage: J = HadamardEJA(5)
  94             sage: J.one()
  95             e0 + e1 + e2 + e3 + e4
  96             sage: x = sum(J.gens())
  97             sage: A = x.subalgebra_generated_by()
  98             sage: A.one()
  99             f0
 100             sage: A.one().superalgebra_element()
 101             e0 + e1 + e2 + e3 + e4
 102
 103         TESTS:
 104
 105         The identity element acts like the identity over the rationals::
 106
 107             sage: set_random_seed()
 108             sage: x = random_eja(field=QQ,orthonormalize=False).random_element()
 109             sage: A = x.subalgebra_generated_by()
 110             sage: x = A.random_element()
 111             sage: A.one()*x == x and x*A.one() == x
 112             True
 113
 114         The identity element acts like the identity over the algebraic
 115         reals with an orthonormal basis::
 116
 117             sage: set_random_seed()
 118             sage: x = random_eja().random_element()
 119             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 120             sage: x = A.random_element()
 121             sage: A.one()*x == x and x*A.one() == x
 122             True
 123
 124         The matrix of the unit element's operator is the identity over
 125         the rationals::
 126
 127             sage: set_random_seed()
 128             sage: x = random_eja(field=QQ,orthonormalize=False).random_element()
 129             sage: A = x.subalgebra_generated_by()
 130             sage: actual = A.one().operator().matrix()
 131             sage: expected = matrix.identity(A.base_ring(), A.dimension())
 132             sage: actual == expected
 133             True
 134
 135         The matrix of the unit element's operator is the identity over
 136         the algebraic reals with an orthonormal basis::
 137
 138             sage: set_random_seed()
 139             sage: x = random_eja().random_element()
 140             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 141             sage: actual = A.one().operator().matrix()
 142             sage: expected = matrix.identity(A.base_ring(), A.dimension())
 143             sage: actual == expected
 144             True
 145
 146         """
 147         if self.dimension() == 0:
 148             return self.zero()
 149
 150         return self(self.superalgebra().one())
 151