mjo/eja/euclidean_jordan_algebra.py

   1 """
   2 Euclidean Jordan Algebras. These are formally-real Jordan Algebras;
   3 specifically those where u^2 + v^2 = 0 implies that u = v = 0. They
   4 are used in optimization, and have some additional nice methods beyond
   5 what can be supported in a general Jordan Algebra.
   6 """
   7
   8 from sage.categories.magmatic_algebras import MagmaticAlgebras
   9 from sage.structure.element import is_Matrix
  10 from sage.structure.category_object import normalize_names
  11
  12 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra import FiniteDimensionalAlgebra
  13 from sage.algebras.finite_dimensional_algebras.finite_dimensional_algebra_element import FiniteDimensionalAlgebraElement
  14
  15 class FiniteDimensionalEuclideanJordanAlgebra(FiniteDimensionalAlgebra):
  16     @staticmethod
  17     def __classcall_private__(cls,
  18                               field,
  19                               mult_table,
  20                               names='e',
  21                               assume_associative=False,
  22                               category=None,
  23                               rank=None):
  24         n = len(mult_table)
  25         mult_table = [b.base_extend(field) for b in mult_table]
  26         for b in mult_table:
  27             b.set_immutable()
  28             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  29                 raise ValueError("input is not a multiplication table")
  30             if not (b.is_symmetric()):
  31                 # Euclidean jordan algebras are commutative, so left/right
  32                 # multiplication is the same.
  33                 raise ValueError("multiplication table must be symmetric")
  34         mult_table = tuple(mult_table)
  35
  36         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  37         cat.or_subcategory(category)
  38         if assume_associative:
  39             cat = cat.Associative()
  40
  41         names = normalize_names(n, names)
  42
  43         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  44         return fda.__classcall__(cls,
  45                                  field,
  46                                  mult_table,
  47                                  assume_associative=assume_associative,
  48                                  names=names,
  49                                  category=cat,
  50                                  rank=rank)
  51
  52
  53     def __init__(self, field,
  54                  mult_table,
  55                  names='e',
  56                  assume_associative=False,
  57                  category=None,
  58                  rank=None):
  59         self._rank = rank
  60         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  61         fda.__init__(field,
  62                      mult_table,
  63                      names=names,
  64                      category=category)
  65
  66
  67     def _repr_(self):
  68         """
  69         Return a string representation of ``self``.
  70         """
  71         fmt = "Euclidean Jordan algebra of degree {} over {}"
  72         return fmt.format(self.degree(), self.base_ring())
  73
  74     def rank(self):
  75         """
  76         Return the rank of this EJA.
  77         """
  78         if self._rank is None:
  79             raise ValueError("no rank specified at genesis")
  80         else:
  81             return self._rank
  82
  83
  84     class Element(FiniteDimensionalAlgebraElement):
  85         """
  86         An element of a Euclidean Jordan algebra.
  87
  88         Since EJAs are commutative, the "right multiplication" matrix is
  89         also the left multiplication matrix and must be symmetric::
  90
  91             sage: set_random_seed()
  92             sage: J = eja_ln(5)
  93             sage: J.random_element().matrix().is_symmetric()
  94             True
  95
  96         """
  97
  98         def __pow__(self, n):
  99             """
 100             Return ``self`` raised to the power ``n``.
 101
 102             Jordan algebras are always power-associative; see for
 103             example Faraut and Koranyi, Proposition II.1.2 (ii).
 104             """
 105             A = self.parent()
 106             if n == 0:
 107                 return A.one()
 108             elif n == 1:
 109                 return self
 110             else:
 111                 return A.element_class(A, self.vector()*(self.matrix()**(n-1)))
 112
 113
 114         def span_of_powers(self):
 115             """
 116             Return the vector space spanned by successive powers of
 117             this element.
 118             """
 119             # The dimension of the subalgebra can't be greater than
 120             # the big algebra, so just put everything into a list
 121             # and let span() get rid of the excess.
 122             V = self.vector().parent()
 123             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 124
 125
 126         def degree(self):
 127             """
 128             Compute the degree of this element the straightforward way
 129             according to the definition; by appending powers to a list
 130             and figuring out its dimension (that is, whether or not
 131             they're linearly dependent).
 132
 133             EXAMPLES::
 134
 135                 sage: J = eja_ln(4)
 136                 sage: J.one().degree()
 137                 1
 138                 sage: e0,e1,e2,e3 = J.gens()
 139                 sage: (e0 - e1).degree()
 140                 2
 141
 142             In the spin factor algebra (of rank two), all elements that
 143             aren't multiples of the identity are regular::
 144
 145                 sage: set_random_seed()
 146                 sage: n = ZZ.random_element(1,10).abs()
 147                 sage: J = eja_ln(n)
 148                 sage: x = J.random_element()
 149                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 150                 True
 151
 152             """
 153             return self.span_of_powers().dimension()
 154
 155
 156         def subalgebra_generated_by(self):
 157             """
 158             Return the subalgebra of the parent EJA generated by this element.
 159             """
 160             # First get the subspace spanned by the powers of myself...
 161             V = self.span_of_powers()
 162             F = self.base_ring()
 163
 164             # Now figure out the entries of the right-multiplication
 165             # matrix for the successive basis elements b0, b1,... of
 166             # that subspace.
 167             mats = []
 168             for b_right in V.basis():
 169                 eja_b_right = self.parent()(b_right)
 170                 b_right_rows = []
 171                 # The first row of the right-multiplication matrix by
 172                 # b1 is what we get if we apply that matrix to b1. The
 173                 # second row of the right multiplication matrix by b1
 174                 # is what we get when we apply that matrix to b2...
 175                 for b_left in V.basis():
 176                     eja_b_left = self.parent()(b_left)
 177                     # Multiply in the original EJA, but then get the
 178                     # coordinates from the subalgebra in terms of its
 179                     # basis.
 180                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 181                     b_right_rows.append(this_row)
 182                 b_right_matrix = matrix(F, b_right_rows)
 183                 mats.append(b_right_matrix)
 184
 185             return FiniteDimensionalEuclideanJordanAlgebra(F, mats)
 186
 187
 188         def minimal_polynomial(self):
 189             """
 190             EXAMPLES::
 191
 192                 sage: set_random_seed()
 193                 sage: n = ZZ.random_element(1,10).abs()
 194                 sage: J = eja_rn(n)
 195                 sage: x = J.random_element()
 196                 sage: x.degree() == x.minimal_polynomial().degree()
 197                 True
 198
 199             ::
 200
 201                 sage: set_random_seed()
 202                 sage: n = ZZ.random_element(1,10).abs()
 203                 sage: J = eja_ln(n)
 204                 sage: x = J.random_element()
 205                 sage: x.degree() == x.minimal_polynomial().degree()
 206                 True
 207
 208             The minimal polynomial and the characteristic polynomial coincide
 209             and are known (see Alizadeh, Example 11.11) for all elements of
 210             the spin factor algebra that aren't scalar multiples of the
 211             identity::
 212
 213                 sage: set_random_seed()
 214                 sage: n = ZZ.random_element(2,10).abs()
 215                 sage: J = eja_ln(n)
 216                 sage: y = J.random_element()
 217                 sage: while y == y.coefficient(0)*J.one():
 218                 ....:     y = J.random_element()
 219                 sage: y0 = y.vector()[0]
 220                 sage: y_bar = y.vector()[1:]
 221                 sage: actual = y.minimal_polynomial()
 222                 sage: x = SR.symbol('x', domain='real')
 223                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 224                 sage: bool(actual == expected)
 225                 True
 226
 227             """
 228             V = self.span_of_powers()
 229             assoc_subalg = self.subalgebra_generated_by()
 230             # Mis-design warning: the basis used for span_of_powers()
 231             # and subalgebra_generated_by() must be the same, and in
 232             # the same order!
 233             subalg_self = assoc_subalg(V.coordinates(self.vector()))
 234             return subalg_self.matrix().minimal_polynomial()
 235
 236
 237         def characteristic_polynomial(self):
 238             return self.matrix().characteristic_polynomial()
 239
 240
 241 def eja_rn(dimension, field=QQ):
 242     """
 243     Return the Euclidean Jordan Algebra corresponding to the set
 244     `R^n` under the Hadamard product.
 245
 246     EXAMPLES:
 247
 248     This multiplication table can be verified by hand::
 249
 250         sage: J = eja_rn(3)
 251         sage: e0,e1,e2 = J.gens()
 252         sage: e0*e0
 253         e0
 254         sage: e0*e1
 255         0
 256         sage: e0*e2
 257         0
 258         sage: e1*e1
 259         e1
 260         sage: e1*e2
 261         0
 262         sage: e2*e2
 263         e2
 264
 265     """
 266     # The FiniteDimensionalAlgebra constructor takes a list of
 267     # matrices, the ith representing right multiplication by the ith
 268     # basis element in the vector space. So if e_1 = (1,0,0), then
 269     # right (Hadamard) multiplication of x by e_1 picks out the first
 270     # component of x; and likewise for the ith basis element e_i.
 271     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 272            for i in xrange(dimension) ]
 273
 274     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 275
 276
 277 def eja_ln(dimension, field=QQ):
 278     """
 279     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 280     cone of the given ``dimension``.
 281
 282     EXAMPLES:
 283
 284     This multiplication table can be verified by hand::
 285
 286         sage: J = eja_ln(4)
 287         sage: e0,e1,e2,e3 = J.gens()
 288         sage: e0*e0
 289         e0
 290         sage: e0*e1
 291         e1
 292         sage: e0*e2
 293         e2
 294         sage: e0*e3
 295         e3
 296         sage: e1*e2
 297         0
 298         sage: e1*e3
 299         0
 300         sage: e2*e3
 301         0
 302
 303     In one dimension, this is the reals under multiplication::
 304
 305       sage: J1 = eja_ln(1)
 306       sage: J2 = eja_rn(1)
 307       sage: J1 == J2
 308       True
 309
 310     """
 311     Qs = []
 312     id_matrix = identity_matrix(field,dimension)
 313     for i in xrange(dimension):
 314         ei = id_matrix.column(i)
 315         Qi = zero_matrix(field,dimension)
 316         Qi.set_row(0, ei)
 317         Qi.set_column(0, ei)
 318         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 319         # The addition of the diagonal matrix adds an extra ei[0] in the
 320         # upper-left corner of the matrix.
 321         Qi[0,0] = Qi[0,0] * ~field(2)
 322         Qs.append(Qi)
 323
 324     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=2)