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         else:
  46             # If we're going to orthonormalize the basis anyway, we
  47             # might as well just do Gram-Schmidt on the whole list of
  48             # powers. The redundant ones will get zero'd out. If this
  49             # looks like a roundabout way to orthonormalize, it is.
  50             # But converting everything from algebra elements to vectors
  51             # to matrices and then back again turns out to be about
  52             # as fast as reimplementing our own Gram-Schmidt that
  53             # works in an EJA.
  54             G,_ = P.gram_schmidt(orthonormal=True)
  55             basis_vectors = [ g for g in G.rows() if not g.is_zero() ]
  56             superalgebra_basis = [ self._superalgebra.from_vector(b)
  57                                    for b in basis_vectors ]
  58
  59         fdeja = super(FiniteDimensionalEuclideanJordanElementSubalgebra, self)
  60         fdeja.__init__(self._superalgebra,
  61                        superalgebra_basis,
  62                        category=category,
  63                        check_axioms=False)
  64
  65         # The rank is the highest possible degree of a minimal
  66         # polynomial, and is bounded above by the dimension. We know
  67         # in this case that there's an element whose minimal
  68         # polynomial has the same degree as the space's dimension
  69         # (remember how we constructed the space?), so that must be
  70         # its rank too.
  71         self.rank.set_cache(self.dimension())
  72
  73
  74     def one(self):
  75         """
  76         Return the multiplicative identity element of this algebra.
  77
  78         The superclass method computes the identity element, which is
  79         beyond overkill in this case: the superalgebra identity
  80         restricted to this algebra is its identity. Note that we can't
  81         count on the first basis element being the identity -- it migth
  82         have been scaled if we orthonormalized the basis.
  83
  84         SETUP::
  85
  86             sage: from mjo.eja.eja_algebra import (HadamardEJA,
  87             ....:                                  random_eja)
  88
  89         EXAMPLES::
  90
  91             sage: J = HadamardEJA(5)
  92             sage: J.one()
  93             e0 + e1 + e2 + e3 + e4
  94             sage: x = sum(J.gens())
  95             sage: A = x.subalgebra_generated_by()
  96             sage: A.one()
  97             f0
  98             sage: A.one().superalgebra_element()
  99             e0 + e1 + e2 + e3 + e4
 100
 101         TESTS:
 102
 103         The identity element acts like the identity over the rationals::
 104
 105             sage: set_random_seed()
 106             sage: x = random_eja(field=QQ).random_element()
 107             sage: A = x.subalgebra_generated_by()
 108             sage: x = A.random_element()
 109             sage: A.one()*x == x and x*A.one() == x
 110             True
 111
 112         The identity element acts like the identity over the algebraic
 113         reals with an orthonormal basis::
 114
 115             sage: set_random_seed()
 116             sage: x = random_eja().random_element()
 117             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 118             sage: x = A.random_element()
 119             sage: A.one()*x == x and x*A.one() == x
 120             True
 121
 122         The matrix of the unit element's operator is the identity over
 123         the rationals::
 124
 125             sage: set_random_seed()
 126             sage: x = random_eja(field=QQ).random_element()
 127             sage: A = x.subalgebra_generated_by()
 128             sage: actual = A.one().operator().matrix()
 129             sage: expected = matrix.identity(A.base_ring(), A.dimension())
 130             sage: actual == expected
 131             True
 132
 133         The matrix of the unit element's operator is the identity over
 134         the algebraic reals with an orthonormal basis::
 135
 136             sage: set_random_seed()
 137             sage: x = random_eja().random_element()
 138             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 139             sage: actual = A.one().operator().matrix()
 140             sage: expected = matrix.identity(A.base_ring(), A.dimension())
 141             sage: actual == expected
 142             True
 143
 144         """
 145         if self.dimension() == 0:
 146             return self.zero()
 147         else:
 148             sa_one = self.superalgebra().one().to_vector()
 149             # The extra hackery is because foo.to_vector() might not
 150             # live in foo.parent().vector_space()!
 151             coords = sum( a*b for (a,b)
 152                           in zip(sa_one,
 153                                  self.superalgebra().vector_space().basis()) )
 154             return self.from_vector(self.vector_space().coordinate_vector(coords))
 155