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