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: n = ZZ.random_element(1,10).abs()
  89             sage: J = eja_rn(5)
  90             sage: J.random_element().matrix().is_symmetric()
  91             True
  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 associative subalgebra of the parent EJA generated
 159             by this element.
 160
 161             TESTS::
 162
 163                 sage: set_random_seed()
 164                 sage: n = ZZ.random_element(1,10).abs()
 165                 sage: J = eja_rn(n)
 166                 sage: x = J.random_element()
 167                 sage: x.subalgebra_generated_by().is_associative()
 168                 True
 169                 sage: J = eja_ln(n)
 170                 sage: x = J.random_element()
 171                 sage: x.subalgebra_generated_by().is_associative()
 172                 True
 173
 174             This is buggy right now::
 175
 176                 sage: J = eja_ln(5)
 177                 sage: x = J.random_element()
 178                 sage: x.matrix()*x.vector() == (x**2).vector() # works
 179                 True
 180                 sage: u = x.subalgebra_generated_by().random_element()
 181                 sage: u.matrix()*u.vector() == (u**2).vector() # busted
 182                 True
 183
 184             """
 185             # First get the subspace spanned by the powers of myself...
 186             V = self.span_of_powers()
 187             F = self.base_ring()
 188
 189             # Now figure out the entries of the right-multiplication
 190             # matrix for the successive basis elements b0, b1,... of
 191             # that subspace.
 192             mats = []
 193             for b_right in V.basis():
 194                 eja_b_right = self.parent()(b_right)
 195                 b_right_rows = []
 196                 # The first row of the right-multiplication matrix by
 197                 # b1 is what we get if we apply that matrix to b1. The
 198                 # second row of the right multiplication matrix by b1
 199                 # is what we get when we apply that matrix to b2...
 200                 for b_left in V.basis():
 201                     eja_b_left = self.parent()(b_left)
 202                     # Multiply in the original EJA, but then get the
 203                     # coordinates from the subalgebra in terms of its
 204                     # basis.
 205                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 206                     b_right_rows.append(this_row)
 207                 b_right_matrix = matrix(F, b_right_rows)
 208                 mats.append(b_right_matrix)
 209
 210             # It's an algebra of polynomials in one element, and EJAs
 211             # are power-associative.
 212             #
 213             # TODO: choose generator names intelligently.
 214             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 215
 216
 217         def minimal_polynomial(self):
 218             """
 219             EXAMPLES::
 220
 221                 sage: set_random_seed()
 222                 sage: n = ZZ.random_element(1,10).abs()
 223                 sage: J = eja_rn(n)
 224                 sage: x = J.random_element()
 225                 sage: x.degree() == x.minimal_polynomial().degree()
 226                 True
 227
 228             ::
 229
 230                 sage: set_random_seed()
 231                 sage: n = ZZ.random_element(1,10).abs()
 232                 sage: J = eja_ln(n)
 233                 sage: x = J.random_element()
 234                 sage: x.degree() == x.minimal_polynomial().degree()
 235                 True
 236
 237             The minimal polynomial and the characteristic polynomial coincide
 238             and are known (see Alizadeh, Example 11.11) for all elements of
 239             the spin factor algebra that aren't scalar multiples of the
 240             identity::
 241
 242                 sage: set_random_seed()
 243                 sage: n = ZZ.random_element(2,10).abs()
 244                 sage: J = eja_ln(n)
 245                 sage: y = J.random_element()
 246                 sage: while y == y.coefficient(0)*J.one():
 247                 ....:     y = J.random_element()
 248                 sage: y0 = y.vector()[0]
 249                 sage: y_bar = y.vector()[1:]
 250                 sage: actual = y.minimal_polynomial()
 251                 sage: x = SR.symbol('x', domain='real')
 252                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 253                 sage: bool(actual == expected)
 254                 True
 255
 256             """
 257             # The element we're going to call "minimal_polynomial()" on.
 258             # Either myself, interpreted as an element of a finite-
 259             # dimensional algebra, or an element of an associative
 260             # subalgebra.
 261             elt = None
 262
 263             if self.parent().is_associative():
 264                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 265             else:
 266                 V = self.span_of_powers()
 267                 assoc_subalg = self.subalgebra_generated_by()
 268                 # Mis-design warning: the basis used for span_of_powers()
 269                 # and subalgebra_generated_by() must be the same, and in
 270                 # the same order!
 271                 elt = assoc_subalg(V.coordinates(self.vector()))
 272
 273             # Recursive call, but should work since elt lives in an
 274             # associative algebra.
 275             return elt.minimal_polynomial()
 276
 277
 278         def is_nilpotent(self):
 279             """
 280             Return whether or not some power of this element is zero.
 281
 282             The superclass method won't work unless we're in an
 283             associative algebra, and we aren't. However, we generate
 284             an assocoative subalgebra and we're nilpotent there if and
 285             only if we're nilpotent here (probably).
 286
 287             TESTS:
 288
 289             The identity element is never nilpotent::
 290
 291                 sage: set_random_seed()
 292                 sage: n = ZZ.random_element(2,10).abs()
 293                 sage: J = eja_rn(n)
 294                 sage: J.one().is_nilpotent()
 295                 False
 296                 sage: J = eja_ln(n)
 297                 sage: J.one().is_nilpotent()
 298                 False
 299
 300             The additive identity is always nilpotent::
 301
 302                 sage: set_random_seed()
 303                 sage: n = ZZ.random_element(2,10).abs()
 304                 sage: J = eja_rn(n)
 305                 sage: J.zero().is_nilpotent()
 306                 True
 307                 sage: J = eja_ln(n)
 308                 sage: J.zero().is_nilpotent()
 309                 True
 310
 311             """
 312             # The element we're going to call "is_nilpotent()" on.
 313             # Either myself, interpreted as an element of a finite-
 314             # dimensional algebra, or an element of an associative
 315             # subalgebra.
 316             elt = None
 317
 318             if self.parent().is_associative():
 319                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 320             else:
 321                 V = self.span_of_powers()
 322                 assoc_subalg = self.subalgebra_generated_by()
 323                 # Mis-design warning: the basis used for span_of_powers()
 324                 # and subalgebra_generated_by() must be the same, and in
 325                 # the same order!
 326                 elt = assoc_subalg(V.coordinates(self.vector()))
 327
 328             # Recursive call, but should work since elt lives in an
 329             # associative algebra.
 330             return elt.is_nilpotent()
 331
 332
 333         def subalgebra_idempotent(self):
 334             """
 335             Find an idempotent in the associative subalgebra I generate
 336             using Proposition 2.3.5 in Baes.
 337             """
 338             if self.is_nilpotent():
 339                 raise ValueError("this only works with non-nilpotent elements!")
 340
 341             V = self.span_of_powers()
 342             J = self.subalgebra_generated_by()
 343             # Mis-design warning: the basis used for span_of_powers()
 344             # and subalgebra_generated_by() must be the same, and in
 345             # the same order!
 346             u = J(V.coordinates(self.vector()))
 347
 348             # The image of the matrix of left-u^m-multiplication
 349             # will be minimal for some natural number s...
 350             s = 0
 351             minimal_dim = V.dimension()
 352             for i in xrange(1, V.dimension()):
 353                 this_dim = (u**i).matrix().image().dimension()
 354                 if this_dim < minimal_dim:
 355                     minimal_dim = this_dim
 356                     s = i
 357
 358             # Now minimal_matrix should correspond to the smallest
 359             # non-zero subspace in Baes's (or really, Koecher's)
 360             # proposition.
 361             #
 362             # However, we need to restrict the matrix to work on the
 363             # subspace... or do we? Can't we just solve, knowing that
 364             # A(c) = u^(s+1) should have a solution in the big space,
 365             # too?
 366             u_next = u**(s+1)
 367             A = u_next.matrix()
 368             c_coordinates = A.solve_right(u_next.vector())
 369
 370             # Now c_coordinates is the idempotent we want, but it's in
 371             # the coordinate system of the subalgebra.
 372             #
 373             # We need the basis for J, but as elements of the parent algebra.
 374             #
 375             #
 376             # TODO: this is buggy, but it's probably because the
 377             # multiplication table for the subalgebra is wrong! The
 378             # matrices should be symmetric I bet.
 379             basis = [self.parent(v) for v in V.basis()]
 380             return self.parent().linear_combination(zip(c_coordinates, basis))
 381
 382
 383
 384         def characteristic_polynomial(self):
 385             return self.matrix().characteristic_polynomial()
 386
 387
 388 def eja_rn(dimension, field=QQ):
 389     """
 390     Return the Euclidean Jordan Algebra corresponding to the set
 391     `R^n` under the Hadamard product.
 392
 393     EXAMPLES:
 394
 395     This multiplication table can be verified by hand::
 396
 397         sage: J = eja_rn(3)
 398         sage: e0,e1,e2 = J.gens()
 399         sage: e0*e0
 400         e0
 401         sage: e0*e1
 402         0
 403         sage: e0*e2
 404         0
 405         sage: e1*e1
 406         e1
 407         sage: e1*e2
 408         0
 409         sage: e2*e2
 410         e2
 411
 412     """
 413     # The FiniteDimensionalAlgebra constructor takes a list of
 414     # matrices, the ith representing right multiplication by the ith
 415     # basis element in the vector space. So if e_1 = (1,0,0), then
 416     # right (Hadamard) multiplication of x by e_1 picks out the first
 417     # component of x; and likewise for the ith basis element e_i.
 418     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 419            for i in xrange(dimension) ]
 420
 421     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 422
 423
 424 def eja_ln(dimension, field=QQ):
 425     """
 426     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 427     cone of the given ``dimension``.
 428
 429     EXAMPLES:
 430
 431     This multiplication table can be verified by hand::
 432
 433         sage: J = eja_ln(4)
 434         sage: e0,e1,e2,e3 = J.gens()
 435         sage: e0*e0
 436         e0
 437         sage: e0*e1
 438         e1
 439         sage: e0*e2
 440         e2
 441         sage: e0*e3
 442         e3
 443         sage: e1*e2
 444         0
 445         sage: e1*e3
 446         0
 447         sage: e2*e3
 448         0
 449
 450     In one dimension, this is the reals under multiplication::
 451
 452       sage: J1 = eja_ln(1)
 453       sage: J2 = eja_rn(1)
 454       sage: J1 == J2
 455       True
 456
 457     """
 458     Qs = []
 459     id_matrix = identity_matrix(field,dimension)
 460     for i in xrange(dimension):
 461         ei = id_matrix.column(i)
 462         Qi = zero_matrix(field,dimension)
 463         Qi.set_row(0, ei)
 464         Qi.set_column(0, ei)
 465         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 466         # The addition of the diagonal matrix adds an extra ei[0] in the
 467         # upper-left corner of the matrix.
 468         Qi[0,0] = Qi[0,0] * ~field(2)
 469         Qs.append(Qi)
 470
 471     # The rank of the spin factor algebra is two, UNLESS we're in a
 472     # one-dimensional ambient space (the rank is bounded by the
 473     # ambient dimension).
 474     rank = min(dimension,2)
 475     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)