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             .. WARNING:
 106
 107                 We have to override this because our superclass uses row vectors
 108                 instead of column vectors! We, on the other hand, assume column
 109                 vectors everywhere.
 110
 111             EXAMPLES:
 112
 113                 sage: set_random_seed()
 114                 sage: J = eja_ln(5)
 115                 sage: x = J.random_element()
 116                 sage: x.matrix()*x.vector() == (x**2).vector()
 117                 True
 118
 119             """
 120             A = self.parent()
 121             if n == 0:
 122                 return A.one()
 123             elif n == 1:
 124                 return self
 125             else:
 126                 return A.element_class(A, (self.matrix()**(n-1))*self.vector())
 127
 128
 129         def characteristic_polynomial(self):
 130             """
 131             Return my characteristic polynomial (if I'm a regular
 132             element).
 133
 134             Eventually this should be implemented in terms of the parent
 135             algebra's characteristic polynomial that works for ALL
 136             elements.
 137             """
 138             if self.is_regular():
 139                 return self.minimal_polynomial()
 140             else:
 141                 return NotImplementedError('irregular element')
 142
 143
 144         def is_nilpotent(self):
 145             """
 146             Return whether or not some power of this element is zero.
 147
 148             The superclass method won't work unless we're in an
 149             associative algebra, and we aren't. However, we generate
 150             an assocoative subalgebra and we're nilpotent there if and
 151             only if we're nilpotent here (probably).
 152
 153             TESTS:
 154
 155             The identity element is never nilpotent::
 156
 157                 sage: set_random_seed()
 158                 sage: n = ZZ.random_element(2,10).abs()
 159                 sage: J = eja_rn(n)
 160                 sage: J.one().is_nilpotent()
 161                 False
 162                 sage: J = eja_ln(n)
 163                 sage: J.one().is_nilpotent()
 164                 False
 165
 166             The additive identity is always nilpotent::
 167
 168                 sage: set_random_seed()
 169                 sage: n = ZZ.random_element(2,10).abs()
 170                 sage: J = eja_rn(n)
 171                 sage: J.zero().is_nilpotent()
 172                 True
 173                 sage: J = eja_ln(n)
 174                 sage: J.zero().is_nilpotent()
 175                 True
 176
 177             """
 178             # The element we're going to call "is_nilpotent()" on.
 179             # Either myself, interpreted as an element of a finite-
 180             # dimensional algebra, or an element of an associative
 181             # subalgebra.
 182             elt = None
 183
 184             if self.parent().is_associative():
 185                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 186             else:
 187                 V = self.span_of_powers()
 188                 assoc_subalg = self.subalgebra_generated_by()
 189                 # Mis-design warning: the basis used for span_of_powers()
 190                 # and subalgebra_generated_by() must be the same, and in
 191                 # the same order!
 192                 elt = assoc_subalg(V.coordinates(self.vector()))
 193
 194             # Recursive call, but should work since elt lives in an
 195             # associative algebra.
 196             return elt.is_nilpotent()
 197
 198
 199         def is_regular(self):
 200             """
 201             Return whether or not this is a regular element.
 202
 203             EXAMPLES:
 204
 205             The identity element always has degree one, but any element
 206             linearly-independent from it is regular::
 207
 208                 sage: J = eja_ln(5)
 209                 sage: J.one().is_regular()
 210                 False
 211                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 212                 sage: for x in J.gens():
 213                 ....:     (J.one() + x).is_regular()
 214                 False
 215                 True
 216                 True
 217                 True
 218                 True
 219
 220             """
 221             return self.degree() == self.parent().rank()
 222
 223
 224         def degree(self):
 225             """
 226             Compute the degree of this element the straightforward way
 227             according to the definition; by appending powers to a list
 228             and figuring out its dimension (that is, whether or not
 229             they're linearly dependent).
 230
 231             EXAMPLES::
 232
 233                 sage: J = eja_ln(4)
 234                 sage: J.one().degree()
 235                 1
 236                 sage: e0,e1,e2,e3 = J.gens()
 237                 sage: (e0 - e1).degree()
 238                 2
 239
 240             In the spin factor algebra (of rank two), all elements that
 241             aren't multiples of the identity are regular::
 242
 243                 sage: set_random_seed()
 244                 sage: n = ZZ.random_element(1,10).abs()
 245                 sage: J = eja_ln(n)
 246                 sage: x = J.random_element()
 247                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 248                 True
 249
 250             """
 251             return self.span_of_powers().dimension()
 252
 253
 254         def matrix(self):
 255             """
 256             Return the matrix that represents left- (or right-)
 257             multiplication by this element in the parent algebra.
 258
 259             We have to override this because the superclass method
 260             returns a matrix that acts on row vectors (that is, on
 261             the right).
 262             """
 263             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 264             return fda_elt.matrix().transpose()
 265
 266
 267         def minimal_polynomial(self):
 268             """
 269             EXAMPLES::
 270
 271                 sage: set_random_seed()
 272                 sage: n = ZZ.random_element(1,10).abs()
 273                 sage: J = eja_rn(n)
 274                 sage: x = J.random_element()
 275                 sage: x.degree() == x.minimal_polynomial().degree()
 276                 True
 277
 278             ::
 279
 280                 sage: set_random_seed()
 281                 sage: n = ZZ.random_element(1,10).abs()
 282                 sage: J = eja_ln(n)
 283                 sage: x = J.random_element()
 284                 sage: x.degree() == x.minimal_polynomial().degree()
 285                 True
 286
 287             The minimal polynomial and the characteristic polynomial coincide
 288             and are known (see Alizadeh, Example 11.11) for all elements of
 289             the spin factor algebra that aren't scalar multiples of the
 290             identity::
 291
 292                 sage: set_random_seed()
 293                 sage: n = ZZ.random_element(2,10).abs()
 294                 sage: J = eja_ln(n)
 295                 sage: y = J.random_element()
 296                 sage: while y == y.coefficient(0)*J.one():
 297                 ....:     y = J.random_element()
 298                 sage: y0 = y.vector()[0]
 299                 sage: y_bar = y.vector()[1:]
 300                 sage: actual = y.minimal_polynomial()
 301                 sage: x = SR.symbol('x', domain='real')
 302                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 303                 sage: bool(actual == expected)
 304                 True
 305
 306             """
 307             # The element we're going to call "minimal_polynomial()" on.
 308             # Either myself, interpreted as an element of a finite-
 309             # dimensional algebra, or an element of an associative
 310             # subalgebra.
 311             elt = None
 312
 313             if self.parent().is_associative():
 314                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 315             else:
 316                 V = self.span_of_powers()
 317                 assoc_subalg = self.subalgebra_generated_by()
 318                 # Mis-design warning: the basis used for span_of_powers()
 319                 # and subalgebra_generated_by() must be the same, and in
 320                 # the same order!
 321                 elt = assoc_subalg(V.coordinates(self.vector()))
 322
 323             # Recursive call, but should work since elt lives in an
 324             # associative algebra.
 325             return elt.minimal_polynomial()
 326
 327
 328         def span_of_powers(self):
 329             """
 330             Return the vector space spanned by successive powers of
 331             this element.
 332             """
 333             # The dimension of the subalgebra can't be greater than
 334             # the big algebra, so just put everything into a list
 335             # and let span() get rid of the excess.
 336             V = self.vector().parent()
 337             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 338
 339
 340         def subalgebra_generated_by(self):
 341             """
 342             Return the associative subalgebra of the parent EJA generated
 343             by this element.
 344
 345             TESTS::
 346
 347                 sage: set_random_seed()
 348                 sage: n = ZZ.random_element(1,10).abs()
 349                 sage: J = eja_rn(n)
 350                 sage: x = J.random_element()
 351                 sage: x.subalgebra_generated_by().is_associative()
 352                 True
 353                 sage: J = eja_ln(n)
 354                 sage: x = J.random_element()
 355                 sage: x.subalgebra_generated_by().is_associative()
 356                 True
 357
 358             Squaring in the subalgebra should be the same thing as
 359             squaring in the superalgebra::
 360
 361                 sage: J = eja_ln(5)
 362                 sage: x = J.random_element()
 363                 sage: u = x.subalgebra_generated_by().random_element()
 364                 sage: u.matrix()*u.vector() == (u**2).vector()
 365                 True
 366
 367             """
 368             # First get the subspace spanned by the powers of myself...
 369             V = self.span_of_powers()
 370             F = self.base_ring()
 371
 372             # Now figure out the entries of the right-multiplication
 373             # matrix for the successive basis elements b0, b1,... of
 374             # that subspace.
 375             mats = []
 376             for b_right in V.basis():
 377                 eja_b_right = self.parent()(b_right)
 378                 b_right_rows = []
 379                 # The first row of the right-multiplication matrix by
 380                 # b1 is what we get if we apply that matrix to b1. The
 381                 # second row of the right multiplication matrix by b1
 382                 # is what we get when we apply that matrix to b2...
 383                 #
 384                 # IMPORTANT: this assumes that all vectors are COLUMN
 385                 # vectors, unlike our superclass (which uses row vectors).
 386                 for b_left in V.basis():
 387                     eja_b_left = self.parent()(b_left)
 388                     # Multiply in the original EJA, but then get the
 389                     # coordinates from the subalgebra in terms of its
 390                     # basis.
 391                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 392                     b_right_rows.append(this_row)
 393                 b_right_matrix = matrix(F, b_right_rows)
 394                 mats.append(b_right_matrix)
 395
 396             # It's an algebra of polynomials in one element, and EJAs
 397             # are power-associative.
 398             #
 399             # TODO: choose generator names intelligently.
 400             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 401
 402
 403         def subalgebra_idempotent(self):
 404             """
 405             Find an idempotent in the associative subalgebra I generate
 406             using Proposition 2.3.5 in Baes.
 407
 408             TESTS::
 409
 410                 sage: set_random_seed()
 411                 sage: J = eja_rn(5)
 412                 sage: c = J.random_element().subalgebra_idempotent()
 413                 sage: c^2 == c
 414                 True
 415                 sage: J = eja_ln(5)
 416                 sage: c = J.random_element().subalgebra_idempotent()
 417                 sage: c^2 == c
 418                 True
 419
 420             """
 421             if self.is_nilpotent():
 422                 raise ValueError("this only works with non-nilpotent elements!")
 423
 424             V = self.span_of_powers()
 425             J = self.subalgebra_generated_by()
 426             # Mis-design warning: the basis used for span_of_powers()
 427             # and subalgebra_generated_by() must be the same, and in
 428             # the same order!
 429             u = J(V.coordinates(self.vector()))
 430
 431             # The image of the matrix of left-u^m-multiplication
 432             # will be minimal for some natural number s...
 433             s = 0
 434             minimal_dim = V.dimension()
 435             for i in xrange(1, V.dimension()):
 436                 this_dim = (u**i).matrix().image().dimension()
 437                 if this_dim < minimal_dim:
 438                     minimal_dim = this_dim
 439                     s = i
 440
 441             # Now minimal_matrix should correspond to the smallest
 442             # non-zero subspace in Baes's (or really, Koecher's)
 443             # proposition.
 444             #
 445             # However, we need to restrict the matrix to work on the
 446             # subspace... or do we? Can't we just solve, knowing that
 447             # A(c) = u^(s+1) should have a solution in the big space,
 448             # too?
 449             #
 450             # Beware, solve_right() means that we're using COLUMN vectors.
 451             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 452             u_next = u**(s+1)
 453             A = u_next.matrix()
 454             c_coordinates = A.solve_right(u_next.vector())
 455
 456             # Now c_coordinates is the idempotent we want, but it's in
 457             # the coordinate system of the subalgebra.
 458             #
 459             # We need the basis for J, but as elements of the parent algebra.
 460             #
 461             basis = [self.parent(v) for v in V.basis()]
 462             return self.parent().linear_combination(zip(c_coordinates, basis))
 463
 464
 465
 466 def eja_rn(dimension, field=QQ):
 467     """
 468     Return the Euclidean Jordan Algebra corresponding to the set
 469     `R^n` under the Hadamard product.
 470
 471     EXAMPLES:
 472
 473     This multiplication table can be verified by hand::
 474
 475         sage: J = eja_rn(3)
 476         sage: e0,e1,e2 = J.gens()
 477         sage: e0*e0
 478         e0
 479         sage: e0*e1
 480         0
 481         sage: e0*e2
 482         0
 483         sage: e1*e1
 484         e1
 485         sage: e1*e2
 486         0
 487         sage: e2*e2
 488         e2
 489
 490     """
 491     # The FiniteDimensionalAlgebra constructor takes a list of
 492     # matrices, the ith representing right multiplication by the ith
 493     # basis element in the vector space. So if e_1 = (1,0,0), then
 494     # right (Hadamard) multiplication of x by e_1 picks out the first
 495     # component of x; and likewise for the ith basis element e_i.
 496     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 497            for i in xrange(dimension) ]
 498
 499     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 500
 501
 502 def eja_ln(dimension, field=QQ):
 503     """
 504     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 505     cone of the given ``dimension``.
 506
 507     EXAMPLES:
 508
 509     This multiplication table can be verified by hand::
 510
 511         sage: J = eja_ln(4)
 512         sage: e0,e1,e2,e3 = J.gens()
 513         sage: e0*e0
 514         e0
 515         sage: e0*e1
 516         e1
 517         sage: e0*e2
 518         e2
 519         sage: e0*e3
 520         e3
 521         sage: e1*e2
 522         0
 523         sage: e1*e3
 524         0
 525         sage: e2*e3
 526         0
 527
 528     In one dimension, this is the reals under multiplication::
 529
 530       sage: J1 = eja_ln(1)
 531       sage: J2 = eja_rn(1)
 532       sage: J1 == J2
 533       True
 534
 535     """
 536     Qs = []
 537     id_matrix = identity_matrix(field,dimension)
 538     for i in xrange(dimension):
 539         ei = id_matrix.column(i)
 540         Qi = zero_matrix(field,dimension)
 541         Qi.set_row(0, ei)
 542         Qi.set_column(0, ei)
 543         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 544         # The addition of the diagonal matrix adds an extra ei[0] in the
 545         # upper-left corner of the matrix.
 546         Qi[0,0] = Qi[0,0] * ~field(2)
 547         Qs.append(Qi)
 548
 549     # The rank of the spin factor algebra is two, UNLESS we're in a
 550     # one-dimensional ambient space (the rank is bounded by the
 551     # ambient dimension).
 552     rank = min(dimension,2)
 553     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)