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
  85         def __pow__(self, n):
  86             """
  87             Return ``self`` raised to the power ``n``.
  88
  89             Jordan algebras are always power-associative; see for
  90             example Faraut and Koranyi, Proposition II.1.2 (ii).
  91
  92             .. WARNING:
  93
  94                 We have to override this because our superclass uses row vectors
  95                 instead of column vectors! We, on the other hand, assume column
  96                 vectors everywhere.
  97
  98             EXAMPLES:
  99
 100                 sage: set_random_seed()
 101                 sage: x = random_eja().random_element()
 102                 sage: x.matrix()*x.vector() == (x**2).vector()
 103                 True
 104
 105             """
 106             A = self.parent()
 107             if n == 0:
 108                 return A.one()
 109             elif n == 1:
 110                 return self
 111             else:
 112                 return A.element_class(A, (self.matrix()**(n-1))*self.vector())
 113
 114
 115         def characteristic_polynomial(self):
 116             """
 117             Return my characteristic polynomial (if I'm a regular
 118             element).
 119
 120             Eventually this should be implemented in terms of the parent
 121             algebra's characteristic polynomial that works for ALL
 122             elements.
 123             """
 124             if self.is_regular():
 125                 return self.minimal_polynomial()
 126             else:
 127                 raise NotImplementedError('irregular element')
 128
 129
 130         def det(self):
 131             """
 132             Return my determinant, the product of my eigenvalues.
 133
 134             EXAMPLES::
 135
 136                 sage: J = eja_ln(2)
 137                 sage: e0,e1 = J.gens()
 138                 sage: x = e0 + e1
 139                 sage: x.det()
 140                 0
 141                 sage: J = eja_ln(3)
 142                 sage: e0,e1,e2 = J.gens()
 143                 sage: x = e0 + e1 + e2
 144                 sage: x.det()
 145                 -1
 146
 147             """
 148             cs = self.characteristic_polynomial().coefficients(sparse=False)
 149             r = len(cs) - 1
 150             if r >= 0:
 151                 return cs[0] * (-1)**r
 152             else:
 153                 raise ValueError('charpoly had no coefficients')
 154
 155
 156         def is_nilpotent(self):
 157             """
 158             Return whether or not some power of this element is zero.
 159
 160             The superclass method won't work unless we're in an
 161             associative algebra, and we aren't. However, we generate
 162             an assocoative subalgebra and we're nilpotent there if and
 163             only if we're nilpotent here (probably).
 164
 165             TESTS:
 166
 167             The identity element is never nilpotent::
 168
 169                 sage: set_random_seed()
 170                 sage: random_eja().one().is_nilpotent()
 171                 False
 172
 173             The additive identity is always nilpotent::
 174
 175                 sage: set_random_seed()
 176                 sage: random_eja().zero().is_nilpotent()
 177                 True
 178
 179             """
 180             # The element we're going to call "is_nilpotent()" on.
 181             # Either myself, interpreted as an element of a finite-
 182             # dimensional algebra, or an element of an associative
 183             # subalgebra.
 184             elt = None
 185
 186             if self.parent().is_associative():
 187                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 188             else:
 189                 V = self.span_of_powers()
 190                 assoc_subalg = self.subalgebra_generated_by()
 191                 # Mis-design warning: the basis used for span_of_powers()
 192                 # and subalgebra_generated_by() must be the same, and in
 193                 # the same order!
 194                 elt = assoc_subalg(V.coordinates(self.vector()))
 195
 196             # Recursive call, but should work since elt lives in an
 197             # associative algebra.
 198             return elt.is_nilpotent()
 199
 200
 201         def is_regular(self):
 202             """
 203             Return whether or not this is a regular element.
 204
 205             EXAMPLES:
 206
 207             The identity element always has degree one, but any element
 208             linearly-independent from it is regular::
 209
 210                 sage: J = eja_ln(5)
 211                 sage: J.one().is_regular()
 212                 False
 213                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 214                 sage: for x in J.gens():
 215                 ....:     (J.one() + x).is_regular()
 216                 False
 217                 True
 218                 True
 219                 True
 220                 True
 221
 222             """
 223             return self.degree() == self.parent().rank()
 224
 225
 226         def degree(self):
 227             """
 228             Compute the degree of this element the straightforward way
 229             according to the definition; by appending powers to a list
 230             and figuring out its dimension (that is, whether or not
 231             they're linearly dependent).
 232
 233             EXAMPLES::
 234
 235                 sage: J = eja_ln(4)
 236                 sage: J.one().degree()
 237                 1
 238                 sage: e0,e1,e2,e3 = J.gens()
 239                 sage: (e0 - e1).degree()
 240                 2
 241
 242             In the spin factor algebra (of rank two), all elements that
 243             aren't multiples of the identity are regular::
 244
 245                 sage: set_random_seed()
 246                 sage: n = ZZ.random_element(1,10).abs()
 247                 sage: J = eja_ln(n)
 248                 sage: x = J.random_element()
 249                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 250                 True
 251
 252             """
 253             return self.span_of_powers().dimension()
 254
 255
 256         def matrix(self):
 257             """
 258             Return the matrix that represents left- (or right-)
 259             multiplication by this element in the parent algebra.
 260
 261             We have to override this because the superclass method
 262             returns a matrix that acts on row vectors (that is, on
 263             the right).
 264             """
 265             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 266             return fda_elt.matrix().transpose()
 267
 268
 269         def minimal_polynomial(self):
 270             """
 271             EXAMPLES::
 272
 273                 sage: set_random_seed()
 274                 sage: x = random_eja().random_element()
 275                 sage: x.degree() == x.minimal_polynomial().degree()
 276                 True
 277
 278             ::
 279
 280                 sage: set_random_seed()
 281                 sage: x = random_eja().random_element()
 282                 sage: x.degree() == x.minimal_polynomial().degree()
 283                 True
 284
 285             The minimal polynomial and the characteristic polynomial coincide
 286             and are known (see Alizadeh, Example 11.11) for all elements of
 287             the spin factor algebra that aren't scalar multiples of the
 288             identity::
 289
 290                 sage: set_random_seed()
 291                 sage: n = ZZ.random_element(2,10).abs()
 292                 sage: J = eja_ln(n)
 293                 sage: y = J.random_element()
 294                 sage: while y == y.coefficient(0)*J.one():
 295                 ....:     y = J.random_element()
 296                 sage: y0 = y.vector()[0]
 297                 sage: y_bar = y.vector()[1:]
 298                 sage: actual = y.minimal_polynomial()
 299                 sage: x = SR.symbol('x', domain='real')
 300                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 301                 sage: bool(actual == expected)
 302                 True
 303
 304             """
 305             # The element we're going to call "minimal_polynomial()" on.
 306             # Either myself, interpreted as an element of a finite-
 307             # dimensional algebra, or an element of an associative
 308             # subalgebra.
 309             elt = None
 310
 311             if self.parent().is_associative():
 312                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 313             else:
 314                 V = self.span_of_powers()
 315                 assoc_subalg = self.subalgebra_generated_by()
 316                 # Mis-design warning: the basis used for span_of_powers()
 317                 # and subalgebra_generated_by() must be the same, and in
 318                 # the same order!
 319                 elt = assoc_subalg(V.coordinates(self.vector()))
 320
 321             # Recursive call, but should work since elt lives in an
 322             # associative algebra.
 323             return elt.minimal_polynomial()
 324
 325
 326         def span_of_powers(self):
 327             """
 328             Return the vector space spanned by successive powers of
 329             this element.
 330             """
 331             # The dimension of the subalgebra can't be greater than
 332             # the big algebra, so just put everything into a list
 333             # and let span() get rid of the excess.
 334             V = self.vector().parent()
 335             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 336
 337
 338         def subalgebra_generated_by(self):
 339             """
 340             Return the associative subalgebra of the parent EJA generated
 341             by this element.
 342
 343             TESTS::
 344
 345                 sage: set_random_seed()
 346                 sage: x = random_eja().random_element()
 347                 sage: x.subalgebra_generated_by().is_associative()
 348                 True
 349
 350             Squaring in the subalgebra should be the same thing as
 351             squaring in the superalgebra::
 352
 353                 sage: set_random_seed()
 354                 sage: x = random_eja().random_element()
 355                 sage: u = x.subalgebra_generated_by().random_element()
 356                 sage: u.matrix()*u.vector() == (u**2).vector()
 357                 True
 358
 359             """
 360             # First get the subspace spanned by the powers of myself...
 361             V = self.span_of_powers()
 362             F = self.base_ring()
 363
 364             # Now figure out the entries of the right-multiplication
 365             # matrix for the successive basis elements b0, b1,... of
 366             # that subspace.
 367             mats = []
 368             for b_right in V.basis():
 369                 eja_b_right = self.parent()(b_right)
 370                 b_right_rows = []
 371                 # The first row of the right-multiplication matrix by
 372                 # b1 is what we get if we apply that matrix to b1. The
 373                 # second row of the right multiplication matrix by b1
 374                 # is what we get when we apply that matrix to b2...
 375                 #
 376                 # IMPORTANT: this assumes that all vectors are COLUMN
 377                 # vectors, unlike our superclass (which uses row vectors).
 378                 for b_left in V.basis():
 379                     eja_b_left = self.parent()(b_left)
 380                     # Multiply in the original EJA, but then get the
 381                     # coordinates from the subalgebra in terms of its
 382                     # basis.
 383                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 384                     b_right_rows.append(this_row)
 385                 b_right_matrix = matrix(F, b_right_rows)
 386                 mats.append(b_right_matrix)
 387
 388             # It's an algebra of polynomials in one element, and EJAs
 389             # are power-associative.
 390             #
 391             # TODO: choose generator names intelligently.
 392             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 393
 394
 395         def subalgebra_idempotent(self):
 396             """
 397             Find an idempotent in the associative subalgebra I generate
 398             using Proposition 2.3.5 in Baes.
 399
 400             TESTS::
 401
 402                 sage: set_random_seed()
 403                 sage: J = eja_rn(5)
 404                 sage: c = J.random_element().subalgebra_idempotent()
 405                 sage: c^2 == c
 406                 True
 407                 sage: J = eja_ln(5)
 408                 sage: c = J.random_element().subalgebra_idempotent()
 409                 sage: c^2 == c
 410                 True
 411
 412             """
 413             if self.is_nilpotent():
 414                 raise ValueError("this only works with non-nilpotent elements!")
 415
 416             V = self.span_of_powers()
 417             J = self.subalgebra_generated_by()
 418             # Mis-design warning: the basis used for span_of_powers()
 419             # and subalgebra_generated_by() must be the same, and in
 420             # the same order!
 421             u = J(V.coordinates(self.vector()))
 422
 423             # The image of the matrix of left-u^m-multiplication
 424             # will be minimal for some natural number s...
 425             s = 0
 426             minimal_dim = V.dimension()
 427             for i in xrange(1, V.dimension()):
 428                 this_dim = (u**i).matrix().image().dimension()
 429                 if this_dim < minimal_dim:
 430                     minimal_dim = this_dim
 431                     s = i
 432
 433             # Now minimal_matrix should correspond to the smallest
 434             # non-zero subspace in Baes's (or really, Koecher's)
 435             # proposition.
 436             #
 437             # However, we need to restrict the matrix to work on the
 438             # subspace... or do we? Can't we just solve, knowing that
 439             # A(c) = u^(s+1) should have a solution in the big space,
 440             # too?
 441             #
 442             # Beware, solve_right() means that we're using COLUMN vectors.
 443             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 444             u_next = u**(s+1)
 445             A = u_next.matrix()
 446             c_coordinates = A.solve_right(u_next.vector())
 447
 448             # Now c_coordinates is the idempotent we want, but it's in
 449             # the coordinate system of the subalgebra.
 450             #
 451             # We need the basis for J, but as elements of the parent algebra.
 452             #
 453             basis = [self.parent(v) for v in V.basis()]
 454             return self.parent().linear_combination(zip(c_coordinates, basis))
 455
 456
 457         def trace(self):
 458             """
 459             Return my trace, the sum of my eigenvalues.
 460
 461             EXAMPLES::
 462
 463                 sage: J = eja_ln(3)
 464                 sage: e0,e1,e2 = J.gens()
 465                 sage: x = e0 + e1 + e2
 466                 sage: x.trace()
 467                 2
 468
 469             """
 470             cs = self.characteristic_polynomial().coefficients(sparse=False)
 471             if len(cs) >= 2:
 472                 return -1*cs[-2]
 473             else:
 474                 raise ValueError('charpoly had fewer than 2 coefficients')
 475
 476
 477 def eja_rn(dimension, field=QQ):
 478     """
 479     Return the Euclidean Jordan Algebra corresponding to the set
 480     `R^n` under the Hadamard product.
 481
 482     EXAMPLES:
 483
 484     This multiplication table can be verified by hand::
 485
 486         sage: J = eja_rn(3)
 487         sage: e0,e1,e2 = J.gens()
 488         sage: e0*e0
 489         e0
 490         sage: e0*e1
 491         0
 492         sage: e0*e2
 493         0
 494         sage: e1*e1
 495         e1
 496         sage: e1*e2
 497         0
 498         sage: e2*e2
 499         e2
 500
 501     """
 502     # The FiniteDimensionalAlgebra constructor takes a list of
 503     # matrices, the ith representing right multiplication by the ith
 504     # basis element in the vector space. So if e_1 = (1,0,0), then
 505     # right (Hadamard) multiplication of x by e_1 picks out the first
 506     # component of x; and likewise for the ith basis element e_i.
 507     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 508            for i in xrange(dimension) ]
 509
 510     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 511
 512
 513 def eja_ln(dimension, field=QQ):
 514     """
 515     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 516     cone of the given ``dimension``.
 517
 518     EXAMPLES:
 519
 520     This multiplication table can be verified by hand::
 521
 522         sage: J = eja_ln(4)
 523         sage: e0,e1,e2,e3 = J.gens()
 524         sage: e0*e0
 525         e0
 526         sage: e0*e1
 527         e1
 528         sage: e0*e2
 529         e2
 530         sage: e0*e3
 531         e3
 532         sage: e1*e2
 533         0
 534         sage: e1*e3
 535         0
 536         sage: e2*e3
 537         0
 538
 539     In one dimension, this is the reals under multiplication::
 540
 541       sage: J1 = eja_ln(1)
 542       sage: J2 = eja_rn(1)
 543       sage: J1 == J2
 544       True
 545
 546     """
 547     Qs = []
 548     id_matrix = identity_matrix(field,dimension)
 549     for i in xrange(dimension):
 550         ei = id_matrix.column(i)
 551         Qi = zero_matrix(field,dimension)
 552         Qi.set_row(0, ei)
 553         Qi.set_column(0, ei)
 554         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 555         # The addition of the diagonal matrix adds an extra ei[0] in the
 556         # upper-left corner of the matrix.
 557         Qi[0,0] = Qi[0,0] * ~field(2)
 558         Qs.append(Qi)
 559
 560     # The rank of the spin factor algebra is two, UNLESS we're in a
 561     # one-dimensional ambient space (the rank is bounded by the
 562     # ambient dimension).
 563     rank = min(dimension,2)
 564     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)
 565
 566
 567 def eja_sn(dimension, field=QQ):
 568     """
 569     Return the simple Jordan algebra of ``dimension``-by-``dimension``
 570     symmetric matrices over ``field``.
 571
 572     EXAMPLES::
 573
 574         sage: J = eja_sn(2)
 575         sage: e0, e1, e2 = J.gens()
 576         sage: e0*e0
 577         e0
 578         sage: e1*e1
 579         e0 + e2
 580         sage: e2*e2
 581         e2
 582
 583     """
 584     Qs = []
 585
 586     # In S^2, for example, we nominally have four coordinates even
 587     # though the space is of dimension three only. The vector space V
 588     # is supposed to hold the entire long vector, and the subspace W
 589     # of V will be spanned by the vectors that arise from symmetric
 590     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
 591     V = VectorSpace(field, dimension**2)
 592
 593     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
 594     # coordinates.
 595     S = []
 596
 597     for i in xrange(dimension):
 598         for j in xrange(i+1):
 599             Eij = matrix(field, dimension, lambda k,l: k==i and l==j)
 600             if i == j:
 601                 Sij = Eij
 602             else:
 603                 Sij = Eij + Eij.transpose()
 604             S.append(Sij)
 605
 606     def mat2vec(m):
 607         return vector(field, m.list())
 608
 609     def vec2mat(v):
 610         return matrix(field, dimension, v.list())
 611
 612     W = V.span( mat2vec(s) for s in S )
 613
 614     # Taking the span above reorders our basis (thanks, jerk!) so we
 615     # need to put our "matrix basis" in the same order as the
 616     # (reordered) vector basis.
 617     S = [ vec2mat(b) for b in W.basis() ]
 618
 619     for s in S:
 620         # Brute force the multiplication-by-s matrix by looping
 621         # through all elements of the basis and doing the computation
 622         # to find out what the corresponding row should be. BEWARE:
 623         # these multiplication tables won't be symmetric! It therefore
 624         # becomes REALLY IMPORTANT that the underlying algebra
 625         # constructor uses ROW vectors and not COLUMN vectors. That's
 626         # why we're computing rows here and not columns.
 627         Q_rows = []
 628         for t in S:
 629             this_row = mat2vec((s*t + t*s)/2)
 630             Q_rows.append(W.coordinates(this_row))
 631         Q = matrix(field,Q_rows)
 632         Qs.append(Q)
 633
 634     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 635
 636
 637 def random_eja():
 638     """
 639     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 640
 641     ALGORITHM:
 642
 643     For now, we choose a random natural number ``n`` (greater than zero)
 644     and then give you back one of the following:
 645
 646       * The cartesian product of the rational numbers ``n`` times; this is
 647         ``QQ^n`` with the Hadamard product.
 648
 649       * The Jordan spin algebra on ``QQ^n``.
 650
 651       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
 652         product.
 653
 654     Later this might be extended to return Cartesian products of the
 655     EJAs above.
 656
 657     TESTS::
 658
 659         sage: random_eja()
 660         Euclidean Jordan algebra of degree...
 661
 662     """
 663     n = ZZ.random_element(1,10).abs()
 664     constructor = choice([eja_rn, eja_ln, eja_sn])
 665     return constructor(dimension=n, field=QQ)