mjo/eja/eja_element_subalgebra.py

   1 from sage.matrix.constructor import matrix
   2 from sage.rings.all import QQ
   3
   4 from mjo.eja.eja_subalgebra import FiniteDimensionalEuclideanJordanSubalgebra
   5
   6
   7 class FiniteDimensionalEuclideanJordanElementSubalgebra(FiniteDimensionalEuclideanJordanSubalgebra):
   8     def __init__(self, elt, orthonormalize_basis):
   9         self._superalgebra = elt.parent()
  10         category = self._superalgebra.category().Associative()
  11         V = self._superalgebra.vector_space()
  12         field = self._superalgebra.base_ring()
  13
  14         # This list is guaranteed to contain all independent powers,
  15         # because it's the maximal set of powers that could possibly
  16         # be independent (by a dimension argument).
  17         powers = [ elt**k for k in range(V.dimension()) ]
  18         power_vectors = [ p.to_vector() for p in powers ]
  19         P = matrix(field, power_vectors)
  20
  21         if orthonormalize_basis == False:
  22             # Echelonize the matrix ourselves, because otherwise the
  23             # call to P.pivot_rows() below can choose a non-optimal
  24             # row-reduction algorithm. In particular, scaling can
  25             # help over AA because it avoids the RecursionError that
  26             # gets thrown when we have to look too hard for a root.
  27             #
  28             # Beware: QQ supports an entirely different set of "algorithm"
  29             # keywords than do AA and RR.
  30             algo = None
  31             if field is not QQ:
  32                 algo = "scaled_partial_pivoting"
  33             P.echelonize(algorithm=algo)
  34
  35             # In this case, we just need to figure out which elements
  36             # of the "powers" list are redundant... First compute the
  37             # vector subspace spanned by the powers of the given
  38             # element.
  39
  40             # Figure out which powers form a linearly-independent set.
  41             ind_rows = P.pivot_rows()
  42
  43             # Pick those out of the list of all powers.
  44             superalgebra_basis = tuple(map(powers.__getitem__, ind_rows))
  45             basis_vectors = map(power_vectors.__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     def one(self):
  76         """
  77         Return the multiplicative identity element of this algebra.
  78
  79         The superclass method computes the identity element, which is
  80         beyond overkill in this case: the superalgebra identity
  81         restricted to this algebra is its identity. Note that we can't
  82         count on the first basis element being the identity -- it migth
  83         have been scaled if we orthonormalized the basis.
  84
  85         SETUP::
  86
  87             sage: from mjo.eja.eja_algebra import (HadamardEJA,
  88             ....:                                  random_eja)
  89
  90         EXAMPLES::
  91
  92             sage: J = HadamardEJA(5)
  93             sage: J.one()
  94             e0 + e1 + e2 + e3 + e4
  95             sage: x = sum(J.gens())
  96             sage: A = x.subalgebra_generated_by()
  97             sage: A.one()
  98             f0
  99             sage: A.one().superalgebra_element()
 100             e0 + e1 + e2 + e3 + e4
 101
 102         TESTS:
 103
 104         The identity element acts like the identity over the rationals::
 105
 106             sage: set_random_seed()
 107             sage: x = random_eja(field=QQ).random_element()
 108             sage: A = x.subalgebra_generated_by()
 109             sage: x = A.random_element()
 110             sage: A.one()*x == x and x*A.one() == x
 111             True
 112
 113         The identity element acts like the identity over the algebraic
 114         reals with an orthonormal basis::
 115
 116             sage: set_random_seed()
 117             sage: x = random_eja().random_element()
 118             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 119             sage: x = A.random_element()
 120             sage: A.one()*x == x and x*A.one() == x
 121             True
 122
 123         The matrix of the unit element's operator is the identity over
 124         the rationals::
 125
 126             sage: set_random_seed()
 127             sage: x = random_eja(field=QQ).random_element()
 128             sage: A = x.subalgebra_generated_by()
 129             sage: actual = A.one().operator().matrix()
 130             sage: expected = matrix.identity(A.base_ring(), A.dimension())
 131             sage: actual == expected
 132             True
 133
 134         The matrix of the unit element's operator is the identity over
 135         the algebraic reals with an orthonormal basis::
 136
 137             sage: set_random_seed()
 138             sage: x = random_eja().random_element()
 139             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 140             sage: actual = A.one().operator().matrix()
 141             sage: expected = matrix.identity(A.base_ring(), A.dimension())
 142             sage: actual == expected
 143             True
 144
 145         """
 146         if self.dimension() == 0:
 147             return self.zero()
 148         else:
 149             sa_one = self.superalgebra().one().to_vector()
 150             # The extra hackery is because foo.to_vector() might not
 151             # live in foo.parent().vector_space()!
 152             coords = sum( a*b for (a,b)
 153                           in zip(sa_one,
 154                                  self.superalgebra().vector_space().basis()) )
 155             return self.from_vector(self.vector_space().coordinate_vector(coords))
 156