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                               natural_basis=None,
  25                               inner_product=None):
  26         n = len(mult_table)
  27         mult_table = [b.base_extend(field) for b in mult_table]
  28         for b in mult_table:
  29             b.set_immutable()
  30             if not (is_Matrix(b) and b.dimensions() == (n, n)):
  31                 raise ValueError("input is not a multiplication table")
  32         mult_table = tuple(mult_table)
  33
  34         cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
  35         cat.or_subcategory(category)
  36         if assume_associative:
  37             cat = cat.Associative()
  38
  39         names = normalize_names(n, names)
  40
  41         fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
  42         return fda.__classcall__(cls,
  43                                  field,
  44                                  mult_table,
  45                                  assume_associative=assume_associative,
  46                                  names=names,
  47                                  category=cat,
  48                                  rank=rank,
  49                                  natural_basis=natural_basis,
  50                                  inner_product=inner_product)
  51
  52
  53     def __init__(self, field,
  54                  mult_table,
  55                  names='e',
  56                  assume_associative=False,
  57                  category=None,
  58                  rank=None,
  59                  natural_basis=None,
  60                  inner_product=None):
  61         """
  62         EXAMPLES:
  63
  64         By definition, Jordan multiplication commutes::
  65
  66             sage: set_random_seed()
  67             sage: J = random_eja()
  68             sage: x = J.random_element()
  69             sage: y = J.random_element()
  70             sage: x*y == y*x
  71             True
  72
  73         """
  74         self._rank = rank
  75         self._natural_basis = natural_basis
  76         self._inner_product = inner_product
  77         fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
  78         fda.__init__(field,
  79                      mult_table,
  80                      names=names,
  81                      category=category)
  82
  83
  84     def _repr_(self):
  85         """
  86         Return a string representation of ``self``.
  87         """
  88         fmt = "Euclidean Jordan algebra of degree {} over {}"
  89         return fmt.format(self.degree(), self.base_ring())
  90
  91
  92     def inner_product(self, x, y):
  93         """
  94         The inner product associated with this Euclidean Jordan algebra.
  95
  96         Will default to the trace inner product if nothing else.
  97
  98         EXAMPLES:
  99
 100         The inner product must satisfy its axiom for this algebra to truly
 101         be a Euclidean Jordan Algebra::
 102
 103             sage: set_random_seed()
 104             sage: J = random_eja()
 105             sage: x = J.random_element()
 106             sage: y = J.random_element()
 107             sage: z = J.random_element()
 108             sage: (x*y).inner_product(z) == y.inner_product(x*z)
 109             True
 110
 111         """
 112         if (not x in self) or (not y in self):
 113             raise TypeError("arguments must live in this algebra")
 114         if self._inner_product is None:
 115             return x.trace_inner_product(y)
 116         else:
 117             return self._inner_product(x,y)
 118
 119
 120     def natural_basis(self):
 121         """
 122         Return a more-natural representation of this algebra's basis.
 123
 124         Every finite-dimensional Euclidean Jordan Algebra is a direct
 125         sum of five simple algebras, four of which comprise Hermitian
 126         matrices. This method returns the original "natural" basis
 127         for our underlying vector space. (Typically, the natural basis
 128         is used to construct the multiplication table in the first place.)
 129
 130         Note that this will always return a matrix. The standard basis
 131         in `R^n` will be returned as `n`-by-`1` column matrices.
 132
 133         EXAMPLES::
 134
 135             sage: J = RealSymmetricSimpleEJA(2)
 136             sage: J.basis()
 137             Family (e0, e1, e2)
 138             sage: J.natural_basis()
 139             (
 140             [1 0]  [0 1]  [0 0]
 141             [0 0], [1 0], [0 1]
 142             )
 143
 144         ::
 145
 146             sage: J = JordanSpinSimpleEJA(2)
 147             sage: J.basis()
 148             Family (e0, e1)
 149             sage: J.natural_basis()
 150             (
 151             [1]  [0]
 152             [0], [1]
 153             )
 154
 155         """
 156         if self._natural_basis is None:
 157             return tuple( b.vector().column() for b in self.basis() )
 158         else:
 159             return self._natural_basis
 160
 161
 162     def rank(self):
 163         """
 164         Return the rank of this EJA.
 165         """
 166         if self._rank is None:
 167             raise ValueError("no rank specified at genesis")
 168         else:
 169             return self._rank
 170
 171
 172     class Element(FiniteDimensionalAlgebraElement):
 173         """
 174         An element of a Euclidean Jordan algebra.
 175         """
 176
 177         def __init__(self, A, elt=None):
 178             """
 179             EXAMPLES:
 180
 181             The identity in `S^n` is converted to the identity in the EJA::
 182
 183                 sage: J = RealSymmetricSimpleEJA(3)
 184                 sage: I = identity_matrix(QQ,3)
 185                 sage: J(I) == J.one()
 186                 True
 187
 188             This skew-symmetric matrix can't be represented in the EJA::
 189
 190                 sage: J = RealSymmetricSimpleEJA(3)
 191                 sage: A = matrix(QQ,3, lambda i,j: i-j)
 192                 sage: J(A)
 193                 Traceback (most recent call last):
 194                 ...
 195                 ArithmeticError: vector is not in free module
 196
 197             """
 198             # Goal: if we're given a matrix, and if it lives in our
 199             # parent algebra's "natural ambient space," convert it
 200             # into an algebra element.
 201             #
 202             # The catch is, we make a recursive call after converting
 203             # the given matrix into a vector that lives in the algebra.
 204             # This we need to try the parent class initializer first,
 205             # to avoid recursing forever if we're given something that
 206             # already fits into the algebra, but also happens to live
 207             # in the parent's "natural ambient space" (this happens with
 208             # vectors in R^n).
 209             try:
 210                 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
 211             except ValueError:
 212                 natural_basis = A.natural_basis()
 213                 if elt in natural_basis[0].matrix_space():
 214                     # Thanks for nothing! Matrix spaces aren't vector
 215                     # spaces in Sage, so we have to figure out its
 216                     # natural-basis coordinates ourselves.
 217                     V = VectorSpace(elt.base_ring(), elt.nrows()**2)
 218                     W = V.span( _mat2vec(s) for s in natural_basis )
 219                     coords =  W.coordinates(_mat2vec(elt))
 220                     FiniteDimensionalAlgebraElement.__init__(self, A, coords)
 221
 222         def __pow__(self, n):
 223             """
 224             Return ``self`` raised to the power ``n``.
 225
 226             Jordan algebras are always power-associative; see for
 227             example Faraut and Koranyi, Proposition II.1.2 (ii).
 228
 229             .. WARNING:
 230
 231                 We have to override this because our superclass uses row vectors
 232                 instead of column vectors! We, on the other hand, assume column
 233                 vectors everywhere.
 234
 235             EXAMPLES::
 236
 237                 sage: set_random_seed()
 238                 sage: x = random_eja().random_element()
 239                 sage: x.operator_matrix()*x.vector() == (x^2).vector()
 240                 True
 241
 242             A few examples of power-associativity::
 243
 244                 sage: set_random_seed()
 245                 sage: x = random_eja().random_element()
 246                 sage: x*(x*x)*(x*x) == x^5
 247                 True
 248                 sage: (x*x)*(x*x*x) == x^5
 249                 True
 250
 251             We also know that powers operator-commute (Koecher, Chapter
 252             III, Corollary 1)::
 253
 254                 sage: set_random_seed()
 255                 sage: x = random_eja().random_element()
 256                 sage: m = ZZ.random_element(0,10)
 257                 sage: n = ZZ.random_element(0,10)
 258                 sage: Lxm = (x^m).operator_matrix()
 259                 sage: Lxn = (x^n).operator_matrix()
 260                 sage: Lxm*Lxn == Lxn*Lxm
 261                 True
 262
 263             """
 264             A = self.parent()
 265             if n == 0:
 266                 return A.one()
 267             elif n == 1:
 268                 return self
 269             else:
 270                 return A( (self.operator_matrix()**(n-1))*self.vector() )
 271
 272
 273         def characteristic_polynomial(self):
 274             """
 275             Return my characteristic polynomial (if I'm a regular
 276             element).
 277
 278             Eventually this should be implemented in terms of the parent
 279             algebra's characteristic polynomial that works for ALL
 280             elements.
 281             """
 282             if self.is_regular():
 283                 return self.minimal_polynomial()
 284             else:
 285                 raise NotImplementedError('irregular element')
 286
 287
 288         def inner_product(self, other):
 289             """
 290             Return the parent algebra's inner product of myself and ``other``.
 291
 292             EXAMPLES:
 293
 294             The inner product in the Jordan spin algebra is the usual
 295             inner product on `R^n` (this example only works because the
 296             basis for the Jordan algebra is the standard basis in `R^n`)::
 297
 298                 sage: J = JordanSpinSimpleEJA(3)
 299                 sage: x = vector(QQ,[1,2,3])
 300                 sage: y = vector(QQ,[4,5,6])
 301                 sage: x.inner_product(y)
 302                 32
 303                 sage: J(x).inner_product(J(y))
 304                 32
 305
 306             The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
 307             multiplication is the usual matrix multiplication in `S^n`,
 308             so the inner product of the identity matrix with itself
 309             should be the `n`::
 310
 311                 sage: J = RealSymmetricSimpleEJA(3)
 312                 sage: J.one().inner_product(J.one())
 313                 3
 314
 315             Likewise, the inner product on `C^n` is `<X,Y> =
 316             Re(trace(X*Y))`, where we must necessarily take the real
 317             part because the product of Hermitian matrices may not be
 318             Hermitian::
 319
 320                 sage: J = ComplexHermitianSimpleEJA(3)
 321                 sage: J.one().inner_product(J.one())
 322                 3
 323
 324             TESTS:
 325
 326             Ensure that we can always compute an inner product, and that
 327             it gives us back a real number::
 328
 329                 sage: set_random_seed()
 330                 sage: J = random_eja()
 331                 sage: x = J.random_element()
 332                 sage: y = J.random_element()
 333                 sage: x.inner_product(y) in RR
 334                 True
 335
 336             """
 337             P = self.parent()
 338             if not other in P:
 339                 raise TypeError("'other' must live in the same algebra")
 340
 341             return P.inner_product(self, other)
 342
 343
 344         def operator_commutes_with(self, other):
 345             """
 346             Return whether or not this element operator-commutes
 347             with ``other``.
 348
 349             EXAMPLES:
 350
 351             The definition of a Jordan algebra says that any element
 352             operator-commutes with its square::
 353
 354                 sage: set_random_seed()
 355                 sage: x = random_eja().random_element()
 356                 sage: x.operator_commutes_with(x^2)
 357                 True
 358
 359             TESTS:
 360
 361             Test Lemma 1 from Chapter III of Koecher::
 362
 363                 sage: set_random_seed()
 364                 sage: J = random_eja()
 365                 sage: u = J.random_element()
 366                 sage: v = J.random_element()
 367                 sage: lhs = u.operator_commutes_with(u*v)
 368                 sage: rhs = v.operator_commutes_with(u^2)
 369                 sage: lhs == rhs
 370                 True
 371
 372             """
 373             if not other in self.parent():
 374                 raise TypeError("'other' must live in the same algebra")
 375
 376             A = self.operator_matrix()
 377             B = other.operator_matrix()
 378             return (A*B == B*A)
 379
 380
 381         def det(self):
 382             """
 383             Return my determinant, the product of my eigenvalues.
 384
 385             EXAMPLES::
 386
 387                 sage: J = JordanSpinSimpleEJA(2)
 388                 sage: e0,e1 = J.gens()
 389                 sage: x = e0 + e1
 390                 sage: x.det()
 391                 0
 392                 sage: J = JordanSpinSimpleEJA(3)
 393                 sage: e0,e1,e2 = J.gens()
 394                 sage: x = e0 + e1 + e2
 395                 sage: x.det()
 396                 -1
 397
 398             """
 399             cs = self.characteristic_polynomial().coefficients(sparse=False)
 400             r = len(cs) - 1
 401             if r >= 0:
 402                 return cs[0] * (-1)**r
 403             else:
 404                 raise ValueError('charpoly had no coefficients')
 405
 406
 407         def inverse(self):
 408             """
 409             Return the Jordan-multiplicative inverse of this element.
 410
 411             We can't use the superclass method because it relies on the
 412             algebra being associative.
 413
 414             EXAMPLES:
 415
 416             The inverse in the spin factor algebra is given in Alizadeh's
 417             Example 11.11::
 418
 419                 sage: set_random_seed()
 420                 sage: n = ZZ.random_element(1,10)
 421                 sage: J = JordanSpinSimpleEJA(n)
 422                 sage: x = J.random_element()
 423                 sage: while x.is_zero():
 424                 ....:     x = J.random_element()
 425                 sage: x_vec = x.vector()
 426                 sage: x0 = x_vec[0]
 427                 sage: x_bar = x_vec[1:]
 428                 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
 429                 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
 430                 sage: x_inverse = coeff*inv_vec
 431                 sage: x.inverse() == J(x_inverse)
 432                 True
 433
 434             TESTS:
 435
 436             The identity element is its own inverse::
 437
 438                 sage: set_random_seed()
 439                 sage: J = random_eja()
 440                 sage: J.one().inverse() == J.one()
 441                 True
 442
 443             If an element has an inverse, it acts like one. TODO: this
 444             can be a lot less ugly once ``is_invertible`` doesn't crash
 445             on irregular elements::
 446
 447                 sage: set_random_seed()
 448                 sage: J = random_eja()
 449                 sage: x = J.random_element()
 450                 sage: try:
 451                 ....:     x.inverse()*x == J.one()
 452                 ....: except:
 453                 ....:     True
 454                 True
 455
 456             """
 457             if self.parent().is_associative():
 458                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 459                 return elt.inverse()
 460
 461             # TODO: we can do better once the call to is_invertible()
 462             # doesn't crash on irregular elements.
 463             #if not self.is_invertible():
 464             #    raise ValueError('element is not invertible')
 465
 466             # We do this a little different than the usual recursive
 467             # call to a finite-dimensional algebra element, because we
 468             # wind up with an inverse that lives in the subalgebra and
 469             # we need information about the parent to convert it back.
 470             V = self.span_of_powers()
 471             assoc_subalg = self.subalgebra_generated_by()
 472             # Mis-design warning: the basis used for span_of_powers()
 473             # and subalgebra_generated_by() must be the same, and in
 474             # the same order!
 475             elt = assoc_subalg(V.coordinates(self.vector()))
 476
 477             # This will be in the subalgebra's coordinates...
 478             fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
 479             subalg_inverse = fda_elt.inverse()
 480
 481             # So we have to convert back...
 482             basis = [ self.parent(v) for v in V.basis() ]
 483             pairs = zip(subalg_inverse.vector(), basis)
 484             return self.parent().linear_combination(pairs)
 485
 486
 487         def is_invertible(self):
 488             """
 489             Return whether or not this element is invertible.
 490
 491             We can't use the superclass method because it relies on
 492             the algebra being associative.
 493             """
 494             return not self.det().is_zero()
 495
 496
 497         def is_nilpotent(self):
 498             """
 499             Return whether or not some power of this element is zero.
 500
 501             The superclass method won't work unless we're in an
 502             associative algebra, and we aren't. However, we generate
 503             an assocoative subalgebra and we're nilpotent there if and
 504             only if we're nilpotent here (probably).
 505
 506             TESTS:
 507
 508             The identity element is never nilpotent::
 509
 510                 sage: set_random_seed()
 511                 sage: random_eja().one().is_nilpotent()
 512                 False
 513
 514             The additive identity is always nilpotent::
 515
 516                 sage: set_random_seed()
 517                 sage: random_eja().zero().is_nilpotent()
 518                 True
 519
 520             """
 521             # The element we're going to call "is_nilpotent()" on.
 522             # Either myself, interpreted as an element of a finite-
 523             # dimensional algebra, or an element of an associative
 524             # subalgebra.
 525             elt = None
 526
 527             if self.parent().is_associative():
 528                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 529             else:
 530                 V = self.span_of_powers()
 531                 assoc_subalg = self.subalgebra_generated_by()
 532                 # Mis-design warning: the basis used for span_of_powers()
 533                 # and subalgebra_generated_by() must be the same, and in
 534                 # the same order!
 535                 elt = assoc_subalg(V.coordinates(self.vector()))
 536
 537             # Recursive call, but should work since elt lives in an
 538             # associative algebra.
 539             return elt.is_nilpotent()
 540
 541
 542         def is_regular(self):
 543             """
 544             Return whether or not this is a regular element.
 545
 546             EXAMPLES:
 547
 548             The identity element always has degree one, but any element
 549             linearly-independent from it is regular::
 550
 551                 sage: J = JordanSpinSimpleEJA(5)
 552                 sage: J.one().is_regular()
 553                 False
 554                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 555                 sage: for x in J.gens():
 556                 ....:     (J.one() + x).is_regular()
 557                 False
 558                 True
 559                 True
 560                 True
 561                 True
 562
 563             """
 564             return self.degree() == self.parent().rank()
 565
 566
 567         def degree(self):
 568             """
 569             Compute the degree of this element the straightforward way
 570             according to the definition; by appending powers to a list
 571             and figuring out its dimension (that is, whether or not
 572             they're linearly dependent).
 573
 574             EXAMPLES::
 575
 576                 sage: J = JordanSpinSimpleEJA(4)
 577                 sage: J.one().degree()
 578                 1
 579                 sage: e0,e1,e2,e3 = J.gens()
 580                 sage: (e0 - e1).degree()
 581                 2
 582
 583             In the spin factor algebra (of rank two), all elements that
 584             aren't multiples of the identity are regular::
 585
 586                 sage: set_random_seed()
 587                 sage: n = ZZ.random_element(1,10)
 588                 sage: J = JordanSpinSimpleEJA(n)
 589                 sage: x = J.random_element()
 590                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 591                 True
 592
 593             """
 594             return self.span_of_powers().dimension()
 595
 596
 597         def minimal_polynomial(self):
 598             """
 599             EXAMPLES::
 600
 601                 sage: set_random_seed()
 602                 sage: x = random_eja().random_element()
 603                 sage: x.degree() == x.minimal_polynomial().degree()
 604                 True
 605
 606             ::
 607
 608                 sage: set_random_seed()
 609                 sage: x = random_eja().random_element()
 610                 sage: x.degree() == x.minimal_polynomial().degree()
 611                 True
 612
 613             The minimal polynomial and the characteristic polynomial coincide
 614             and are known (see Alizadeh, Example 11.11) for all elements of
 615             the spin factor algebra that aren't scalar multiples of the
 616             identity::
 617
 618                 sage: set_random_seed()
 619                 sage: n = ZZ.random_element(2,10)
 620                 sage: J = JordanSpinSimpleEJA(n)
 621                 sage: y = J.random_element()
 622                 sage: while y == y.coefficient(0)*J.one():
 623                 ....:     y = J.random_element()
 624                 sage: y0 = y.vector()[0]
 625                 sage: y_bar = y.vector()[1:]
 626                 sage: actual = y.minimal_polynomial()
 627                 sage: x = SR.symbol('x', domain='real')
 628                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 629                 sage: bool(actual == expected)
 630                 True
 631
 632             """
 633             # The element we're going to call "minimal_polynomial()" on.
 634             # Either myself, interpreted as an element of a finite-
 635             # dimensional algebra, or an element of an associative
 636             # subalgebra.
 637             elt = None
 638
 639             if self.parent().is_associative():
 640                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 641             else:
 642                 V = self.span_of_powers()
 643                 assoc_subalg = self.subalgebra_generated_by()
 644                 # Mis-design warning: the basis used for span_of_powers()
 645                 # and subalgebra_generated_by() must be the same, and in
 646                 # the same order!
 647                 elt = assoc_subalg(V.coordinates(self.vector()))
 648
 649             # Recursive call, but should work since elt lives in an
 650             # associative algebra.
 651             return elt.minimal_polynomial()
 652
 653
 654         def natural_representation(self):
 655             """
 656             Return a more-natural representation of this element.
 657
 658             Every finite-dimensional Euclidean Jordan Algebra is a
 659             direct sum of five simple algebras, four of which comprise
 660             Hermitian matrices. This method returns the original
 661             "natural" representation of this element as a Hermitian
 662             matrix, if it has one. If not, you get the usual representation.
 663
 664             EXAMPLES::
 665
 666                 sage: J = ComplexHermitianSimpleEJA(3)
 667                 sage: J.one()
 668                 e0 + e5 + e8
 669                 sage: J.one().natural_representation()
 670                 [1 0 0 0 0 0]
 671                 [0 1 0 0 0 0]
 672                 [0 0 1 0 0 0]
 673                 [0 0 0 1 0 0]
 674                 [0 0 0 0 1 0]
 675                 [0 0 0 0 0 1]
 676
 677             """
 678             B = self.parent().natural_basis()
 679             W = B[0].matrix_space()
 680             return W.linear_combination(zip(self.vector(), B))
 681
 682
 683         def operator_matrix(self):
 684             """
 685             Return the matrix that represents left- (or right-)
 686             multiplication by this element in the parent algebra.
 687
 688             We have to override this because the superclass method
 689             returns a matrix that acts on row vectors (that is, on
 690             the right).
 691
 692             EXAMPLES:
 693
 694             Test the first polarization identity from my notes, Koecher Chapter
 695             III, or from Baes (2.3)::
 696
 697                 sage: set_random_seed()
 698                 sage: J = random_eja()
 699                 sage: x = J.random_element()
 700                 sage: y = J.random_element()
 701                 sage: Lx = x.operator_matrix()
 702                 sage: Ly = y.operator_matrix()
 703                 sage: Lxx = (x*x).operator_matrix()
 704                 sage: Lxy = (x*y).operator_matrix()
 705                 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
 706                 True
 707
 708             Test the second polarization identity from my notes or from
 709             Baes (2.4)::
 710
 711                 sage: set_random_seed()
 712                 sage: J = random_eja()
 713                 sage: x = J.random_element()
 714                 sage: y = J.random_element()
 715                 sage: z = J.random_element()
 716                 sage: Lx = x.operator_matrix()
 717                 sage: Ly = y.operator_matrix()
 718                 sage: Lz = z.operator_matrix()
 719                 sage: Lzy = (z*y).operator_matrix()
 720                 sage: Lxy = (x*y).operator_matrix()
 721                 sage: Lxz = (x*z).operator_matrix()
 722                 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
 723                 True
 724
 725             Test the third polarization identity from my notes or from
 726             Baes (2.5)::
 727
 728                 sage: set_random_seed()
 729                 sage: J = random_eja()
 730                 sage: u = J.random_element()
 731                 sage: y = J.random_element()
 732                 sage: z = J.random_element()
 733                 sage: Lu = u.operator_matrix()
 734                 sage: Ly = y.operator_matrix()
 735                 sage: Lz = z.operator_matrix()
 736                 sage: Lzy = (z*y).operator_matrix()
 737                 sage: Luy = (u*y).operator_matrix()
 738                 sage: Luz = (u*z).operator_matrix()
 739                 sage: Luyz = (u*(y*z)).operator_matrix()
 740                 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
 741                 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
 742                 sage: bool(lhs == rhs)
 743                 True
 744
 745             """
 746             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 747             return fda_elt.matrix().transpose()
 748
 749
 750         def quadratic_representation(self, other=None):
 751             """
 752             Return the quadratic representation of this element.
 753
 754             EXAMPLES:
 755
 756             The explicit form in the spin factor algebra is given by
 757             Alizadeh's Example 11.12::
 758
 759                 sage: set_random_seed()
 760                 sage: n = ZZ.random_element(1,10)
 761                 sage: J = JordanSpinSimpleEJA(n)
 762                 sage: x = J.random_element()
 763                 sage: x_vec = x.vector()
 764                 sage: x0 = x_vec[0]
 765                 sage: x_bar = x_vec[1:]
 766                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 767                 sage: B = 2*x0*x_bar.row()
 768                 sage: C = 2*x0*x_bar.column()
 769                 sage: D = identity_matrix(QQ, n-1)
 770                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 771                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 772                 sage: Q = block_matrix(2,2,[A,B,C,D])
 773                 sage: Q == x.quadratic_representation()
 774                 True
 775
 776             Test all of the properties from Theorem 11.2 in Alizadeh::
 777
 778                 sage: set_random_seed()
 779                 sage: J = random_eja()
 780                 sage: x = J.random_element()
 781                 sage: y = J.random_element()
 782
 783             Property 1:
 784
 785                 sage: actual = x.quadratic_representation(y)
 786                 sage: expected = ( (x+y).quadratic_representation()
 787                 ....:              -x.quadratic_representation()
 788                 ....:              -y.quadratic_representation() ) / 2
 789                 sage: actual == expected
 790                 True
 791
 792             Property 2:
 793
 794                 sage: alpha = QQ.random_element()
 795                 sage: actual = (alpha*x).quadratic_representation()
 796                 sage: expected = (alpha^2)*x.quadratic_representation()
 797                 sage: actual == expected
 798                 True
 799
 800             Property 5:
 801
 802                 sage: Qy = y.quadratic_representation()
 803                 sage: actual = J(Qy*x.vector()).quadratic_representation()
 804                 sage: expected = Qy*x.quadratic_representation()*Qy
 805                 sage: actual == expected
 806                 True
 807
 808             Property 6:
 809
 810                 sage: k = ZZ.random_element(1,10)
 811                 sage: actual = (x^k).quadratic_representation()
 812                 sage: expected = (x.quadratic_representation())^k
 813                 sage: actual == expected
 814                 True
 815
 816             """
 817             if other is None:
 818                 other=self
 819             elif not other in self.parent():
 820                 raise TypeError("'other' must live in the same algebra")
 821
 822             L = self.operator_matrix()
 823             M = other.operator_matrix()
 824             return ( L*M + M*L - (self*other).operator_matrix() )
 825
 826
 827         def span_of_powers(self):
 828             """
 829             Return the vector space spanned by successive powers of
 830             this element.
 831             """
 832             # The dimension of the subalgebra can't be greater than
 833             # the big algebra, so just put everything into a list
 834             # and let span() get rid of the excess.
 835             V = self.vector().parent()
 836             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 837
 838
 839         def subalgebra_generated_by(self):
 840             """
 841             Return the associative subalgebra of the parent EJA generated
 842             by this element.
 843
 844             TESTS::
 845
 846                 sage: set_random_seed()
 847                 sage: x = random_eja().random_element()
 848                 sage: x.subalgebra_generated_by().is_associative()
 849                 True
 850
 851             Squaring in the subalgebra should be the same thing as
 852             squaring in the superalgebra::
 853
 854                 sage: set_random_seed()
 855                 sage: x = random_eja().random_element()
 856                 sage: u = x.subalgebra_generated_by().random_element()
 857                 sage: u.operator_matrix()*u.vector() == (u**2).vector()
 858                 True
 859
 860             """
 861             # First get the subspace spanned by the powers of myself...
 862             V = self.span_of_powers()
 863             F = self.base_ring()
 864
 865             # Now figure out the entries of the right-multiplication
 866             # matrix for the successive basis elements b0, b1,... of
 867             # that subspace.
 868             mats = []
 869             for b_right in V.basis():
 870                 eja_b_right = self.parent()(b_right)
 871                 b_right_rows = []
 872                 # The first row of the right-multiplication matrix by
 873                 # b1 is what we get if we apply that matrix to b1. The
 874                 # second row of the right multiplication matrix by b1
 875                 # is what we get when we apply that matrix to b2...
 876                 #
 877                 # IMPORTANT: this assumes that all vectors are COLUMN
 878                 # vectors, unlike our superclass (which uses row vectors).
 879                 for b_left in V.basis():
 880                     eja_b_left = self.parent()(b_left)
 881                     # Multiply in the original EJA, but then get the
 882                     # coordinates from the subalgebra in terms of its
 883                     # basis.
 884                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 885                     b_right_rows.append(this_row)
 886                 b_right_matrix = matrix(F, b_right_rows)
 887                 mats.append(b_right_matrix)
 888
 889             # It's an algebra of polynomials in one element, and EJAs
 890             # are power-associative.
 891             #
 892             # TODO: choose generator names intelligently.
 893             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 894
 895
 896         def subalgebra_idempotent(self):
 897             """
 898             Find an idempotent in the associative subalgebra I generate
 899             using Proposition 2.3.5 in Baes.
 900
 901             TESTS::
 902
 903                 sage: set_random_seed()
 904                 sage: J = eja_rn(5)
 905                 sage: c = J.random_element().subalgebra_idempotent()
 906                 sage: c^2 == c
 907                 True
 908                 sage: J = JordanSpinSimpleEJA(5)
 909                 sage: c = J.random_element().subalgebra_idempotent()
 910                 sage: c^2 == c
 911                 True
 912
 913             """
 914             if self.is_nilpotent():
 915                 raise ValueError("this only works with non-nilpotent elements!")
 916
 917             V = self.span_of_powers()
 918             J = self.subalgebra_generated_by()
 919             # Mis-design warning: the basis used for span_of_powers()
 920             # and subalgebra_generated_by() must be the same, and in
 921             # the same order!
 922             u = J(V.coordinates(self.vector()))
 923
 924             # The image of the matrix of left-u^m-multiplication
 925             # will be minimal for some natural number s...
 926             s = 0
 927             minimal_dim = V.dimension()
 928             for i in xrange(1, V.dimension()):
 929                 this_dim = (u**i).operator_matrix().image().dimension()
 930                 if this_dim < minimal_dim:
 931                     minimal_dim = this_dim
 932                     s = i
 933
 934             # Now minimal_matrix should correspond to the smallest
 935             # non-zero subspace in Baes's (or really, Koecher's)
 936             # proposition.
 937             #
 938             # However, we need to restrict the matrix to work on the
 939             # subspace... or do we? Can't we just solve, knowing that
 940             # A(c) = u^(s+1) should have a solution in the big space,
 941             # too?
 942             #
 943             # Beware, solve_right() means that we're using COLUMN vectors.
 944             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 945             u_next = u**(s+1)
 946             A = u_next.operator_matrix()
 947             c_coordinates = A.solve_right(u_next.vector())
 948
 949             # Now c_coordinates is the idempotent we want, but it's in
 950             # the coordinate system of the subalgebra.
 951             #
 952             # We need the basis for J, but as elements of the parent algebra.
 953             #
 954             basis = [self.parent(v) for v in V.basis()]
 955             return self.parent().linear_combination(zip(c_coordinates, basis))
 956
 957
 958         def trace(self):
 959             """
 960             Return my trace, the sum of my eigenvalues.
 961
 962             EXAMPLES::
 963
 964                 sage: J = JordanSpinSimpleEJA(3)
 965                 sage: e0,e1,e2 = J.gens()
 966                 sage: x = e0 + e1 + e2
 967                 sage: x.trace()
 968                 2
 969
 970             """
 971             cs = self.characteristic_polynomial().coefficients(sparse=False)
 972             if len(cs) >= 2:
 973                 return -1*cs[-2]
 974             else:
 975                 raise ValueError('charpoly had fewer than 2 coefficients')
 976
 977
 978         def trace_inner_product(self, other):
 979             """
 980             Return the trace inner product of myself and ``other``.
 981             """
 982             if not other in self.parent():
 983                 raise TypeError("'other' must live in the same algebra")
 984
 985             return (self*other).trace()
 986
 987
 988 def eja_rn(dimension, field=QQ):
 989     """
 990     Return the Euclidean Jordan Algebra corresponding to the set
 991     `R^n` under the Hadamard product.
 992
 993     EXAMPLES:
 994
 995     This multiplication table can be verified by hand::
 996
 997         sage: J = eja_rn(3)
 998         sage: e0,e1,e2 = J.gens()
 999         sage: e0*e0
1000         e0
1001         sage: e0*e1
1002         0
1003         sage: e0*e2
1004         0
1005         sage: e1*e1
1006         e1
1007         sage: e1*e2
1008         0
1009         sage: e2*e2
1010         e2
1011
1012     """
1013     # The FiniteDimensionalAlgebra constructor takes a list of
1014     # matrices, the ith representing right multiplication by the ith
1015     # basis element in the vector space. So if e_1 = (1,0,0), then
1016     # right (Hadamard) multiplication of x by e_1 picks out the first
1017     # component of x; and likewise for the ith basis element e_i.
1018     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
1019            for i in xrange(dimension) ]
1020
1021     return FiniteDimensionalEuclideanJordanAlgebra(field,
1022                                                    Qs,
1023                                                    rank=dimension,
1024                                                    inner_product=_usual_ip)
1025
1026
1027
1028 def random_eja():
1029     """
1030     Return a "random" finite-dimensional Euclidean Jordan Algebra.
1031
1032     ALGORITHM:
1033
1034     For now, we choose a random natural number ``n`` (greater than zero)
1035     and then give you back one of the following:
1036
1037       * The cartesian product of the rational numbers ``n`` times; this is
1038         ``QQ^n`` with the Hadamard product.
1039
1040       * The Jordan spin algebra on ``QQ^n``.
1041
1042       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1043         product.
1044
1045     Later this might be extended to return Cartesian products of the
1046     EJAs above.
1047
1048     TESTS::
1049
1050         sage: random_eja()
1051         Euclidean Jordan algebra of degree...
1052
1053     """
1054     n = ZZ.random_element(1,5)
1055     constructor = choice([eja_rn,
1056                           JordanSpinSimpleEJA,
1057                           RealSymmetricSimpleEJA,
1058                           ComplexHermitianSimpleEJA])
1059     return constructor(n, field=QQ)
1060
1061
1062
1063 def _real_symmetric_basis(n, field=QQ):
1064     """
1065     Return a basis for the space of real symmetric n-by-n matrices.
1066     """
1067     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1068     # coordinates.
1069     S = []
1070     for i in xrange(n):
1071         for j in xrange(i+1):
1072             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1073             if i == j:
1074                 Sij = Eij
1075             else:
1076                 # Beware, orthogonal but not normalized!
1077                 Sij = Eij + Eij.transpose()
1078             S.append(Sij)
1079     return tuple(S)
1080
1081
1082 def _complex_hermitian_basis(n, field=QQ):
1083     """
1084     Returns a basis for the space of complex Hermitian n-by-n matrices.
1085
1086     TESTS::
1087
1088         sage: set_random_seed()
1089         sage: n = ZZ.random_element(1,5)
1090         sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1091         True
1092
1093     """
1094     F = QuadraticField(-1, 'I')
1095     I = F.gen()
1096
1097     # This is like the symmetric case, but we need to be careful:
1098     #
1099     #   * We want conjugate-symmetry, not just symmetry.
1100     #   * The diagonal will (as a result) be real.
1101     #
1102     S = []
1103     for i in xrange(n):
1104         for j in xrange(i+1):
1105             Eij = matrix(field, n, lambda k,l: k==i and l==j)
1106             if i == j:
1107                 Sij = _embed_complex_matrix(Eij)
1108                 S.append(Sij)
1109             else:
1110                 # Beware, orthogonal but not normalized! The second one
1111                 # has a minus because it's conjugated.
1112                 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1113                 S.append(Sij_real)
1114                 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1115                 S.append(Sij_imag)
1116     return tuple(S)
1117
1118
1119 def _mat2vec(m):
1120         return vector(m.base_ring(), m.list())
1121
1122 def _vec2mat(v):
1123         return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1124
1125 def _multiplication_table_from_matrix_basis(basis):
1126     """
1127     At least three of the five simple Euclidean Jordan algebras have the
1128     symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1129     multiplication on the right is matrix multiplication. Given a basis
1130     for the underlying matrix space, this function returns a
1131     multiplication table (obtained by looping through the basis
1132     elements) for an algebra of those matrices. A reordered copy
1133     of the basis is also returned to work around the fact that
1134     the ``span()`` in this function will change the order of the basis
1135     from what we think it is, to... something else.
1136     """
1137     # In S^2, for example, we nominally have four coordinates even
1138     # though the space is of dimension three only. The vector space V
1139     # is supposed to hold the entire long vector, and the subspace W
1140     # of V will be spanned by the vectors that arise from symmetric
1141     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1142     field = basis[0].base_ring()
1143     dimension = basis[0].nrows()
1144
1145     V = VectorSpace(field, dimension**2)
1146     W = V.span( _mat2vec(s) for s in basis )
1147
1148     # Taking the span above reorders our basis (thanks, jerk!) so we
1149     # need to put our "matrix basis" in the same order as the
1150     # (reordered) vector basis.
1151     S = tuple( _vec2mat(b) for b in W.basis() )
1152
1153     Qs = []
1154     for s in S:
1155         # Brute force the multiplication-by-s matrix by looping
1156         # through all elements of the basis and doing the computation
1157         # to find out what the corresponding row should be. BEWARE:
1158         # these multiplication tables won't be symmetric! It therefore
1159         # becomes REALLY IMPORTANT that the underlying algebra
1160         # constructor uses ROW vectors and not COLUMN vectors. That's
1161         # why we're computing rows here and not columns.
1162         Q_rows = []
1163         for t in S:
1164             this_row = _mat2vec((s*t + t*s)/2)
1165             Q_rows.append(W.coordinates(this_row))
1166         Q = matrix(field, W.dimension(), Q_rows)
1167         Qs.append(Q)
1168
1169     return (Qs, S)
1170
1171
1172 def _embed_complex_matrix(M):
1173     """
1174     Embed the n-by-n complex matrix ``M`` into the space of real
1175     matrices of size 2n-by-2n via the map the sends each entry `z = a +
1176     bi` to the block matrix ``[[a,b],[-b,a]]``.
1177
1178     EXAMPLES::
1179
1180         sage: F = QuadraticField(-1,'i')
1181         sage: x1 = F(4 - 2*i)
1182         sage: x2 = F(1 + 2*i)
1183         sage: x3 = F(-i)
1184         sage: x4 = F(6)
1185         sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1186         sage: _embed_complex_matrix(M)
1187         [ 4 -2| 1  2]
1188         [ 2  4|-2  1]
1189         [-----+-----]
1190         [ 0 -1| 6  0]
1191         [ 1  0| 0  6]
1192
1193     """
1194     n = M.nrows()
1195     if M.ncols() != n:
1196         raise ValueError("the matrix 'M' must be square")
1197     field = M.base_ring()
1198     blocks = []
1199     for z in M.list():
1200         a = z.real()
1201         b = z.imag()
1202         blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1203
1204     # We can drop the imaginaries here.
1205     return block_matrix(field.base_ring(), n, blocks)
1206
1207
1208 def _unembed_complex_matrix(M):
1209     """
1210     The inverse of _embed_complex_matrix().
1211
1212     EXAMPLES::
1213
1214         sage: A = matrix(QQ,[ [ 1,  2,   3,  4],
1215         ....:                 [-2,  1,  -4,  3],
1216         ....:                 [ 9,  10, 11, 12],
1217         ....:                 [-10, 9, -12, 11] ])
1218         sage: _unembed_complex_matrix(A)
1219         [  -2*i + 1   -4*i + 3]
1220         [ -10*i + 9 -12*i + 11]
1221     """
1222     n = ZZ(M.nrows())
1223     if M.ncols() != n:
1224         raise ValueError("the matrix 'M' must be square")
1225     if not n.mod(2).is_zero():
1226         raise ValueError("the matrix 'M' must be a complex embedding")
1227
1228     F = QuadraticField(-1, 'i')
1229     i = F.gen()
1230
1231     # Go top-left to bottom-right (reading order), converting every
1232     # 2-by-2 block we see to a single complex element.
1233     elements = []
1234     for k in xrange(n/2):
1235         for j in xrange(n/2):
1236             submat = M[2*k:2*k+2,2*j:2*j+2]
1237             if submat[0,0] != submat[1,1]:
1238                 raise ValueError('bad real submatrix')
1239             if submat[0,1] != -submat[1,0]:
1240                 raise ValueError('bad imag submatrix')
1241             z = submat[0,0] + submat[1,0]*i
1242             elements.append(z)
1243
1244     return matrix(F, n/2, elements)
1245
1246 # The usual inner product on R^n.
1247 def _usual_ip(x,y):
1248     return x.vector().inner_product(y.vector())
1249
1250 # The inner product used for the real symmetric simple EJA.
1251 # We keep it as a separate function because e.g. the complex
1252 # algebra uses the same inner product, except divided by 2.
1253 def _matrix_ip(X,Y):
1254     X_mat = X.natural_representation()
1255     Y_mat = Y.natural_representation()
1256     return (X_mat*Y_mat).trace()
1257
1258
1259 def RealSymmetricSimpleEJA(n, field=QQ):
1260     """
1261     The rank-n simple EJA consisting of real symmetric n-by-n
1262     matrices, the usual symmetric Jordan product, and the trace inner
1263     product. It has dimension `(n^2 + n)/2` over the reals.
1264
1265     EXAMPLES::
1266
1267         sage: J = RealSymmetricSimpleEJA(2)
1268         sage: e0, e1, e2 = J.gens()
1269         sage: e0*e0
1270         e0
1271         sage: e1*e1
1272         e0 + e2
1273         sage: e2*e2
1274         e2
1275
1276     TESTS:
1277
1278     The degree of this algebra is `(n^2 + n) / 2`::
1279
1280         sage: set_random_seed()
1281         sage: n = ZZ.random_element(1,5)
1282         sage: J = RealSymmetricSimpleEJA(n)
1283         sage: J.degree() == (n^2 + n)/2
1284         True
1285
1286     """
1287     S = _real_symmetric_basis(n, field=field)
1288     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1289
1290     return FiniteDimensionalEuclideanJordanAlgebra(field,
1291                                                    Qs,
1292                                                    rank=n,
1293                                                    natural_basis=T,
1294                                                    inner_product=_matrix_ip)
1295
1296
1297 def ComplexHermitianSimpleEJA(n, field=QQ):
1298     """
1299     The rank-n simple EJA consisting of complex Hermitian n-by-n
1300     matrices over the real numbers, the usual symmetric Jordan product,
1301     and the real-part-of-trace inner product. It has dimension `n^2` over
1302     the reals.
1303
1304     TESTS:
1305
1306     The degree of this algebra is `n^2`::
1307
1308         sage: set_random_seed()
1309         sage: n = ZZ.random_element(1,5)
1310         sage: J = ComplexHermitianSimpleEJA(n)
1311         sage: J.degree() == n^2
1312         True
1313
1314     """
1315     S = _complex_hermitian_basis(n)
1316     (Qs, T) = _multiplication_table_from_matrix_basis(S)
1317
1318     # Since a+bi on the diagonal is represented as
1319     #
1320     #   a + bi  = [  a  b  ]
1321     #             [ -b  a  ],
1322     #
1323     # we'll double-count the "a" entries if we take the trace of
1324     # the embedding.
1325     ip = lambda X,Y: _matrix_ip(X,Y)/2
1326
1327     return FiniteDimensionalEuclideanJordanAlgebra(field,
1328                                                    Qs,
1329                                                    rank=n,
1330                                                    natural_basis=T,
1331                                                    inner_product=ip)
1332
1333
1334 def QuaternionHermitianSimpleEJA(n):
1335     """
1336     The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1337     matrices, the usual symmetric Jordan product, and the
1338     real-part-of-trace inner product. It has dimension `2n^2 - n` over
1339     the reals.
1340     """
1341     pass
1342
1343 def OctonionHermitianSimpleEJA(n):
1344     """
1345     This shit be crazy. It has dimension 27 over the reals.
1346     """
1347     n = 3
1348     pass
1349
1350 def JordanSpinSimpleEJA(n, field=QQ):
1351     """
1352     The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1353     with the usual inner product and jordan product ``x*y =
1354     (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1355     the reals.
1356
1357     EXAMPLES:
1358
1359     This multiplication table can be verified by hand::
1360
1361         sage: J = JordanSpinSimpleEJA(4)
1362         sage: e0,e1,e2,e3 = J.gens()
1363         sage: e0*e0
1364         e0
1365         sage: e0*e1
1366         e1
1367         sage: e0*e2
1368         e2
1369         sage: e0*e3
1370         e3
1371         sage: e1*e2
1372         0
1373         sage: e1*e3
1374         0
1375         sage: e2*e3
1376         0
1377
1378     In one dimension, this is the reals under multiplication::
1379
1380         sage: J1 = JordanSpinSimpleEJA(1)
1381         sage: J2 = eja_rn(1)
1382         sage: J1 == J2
1383         True
1384
1385     """
1386     Qs = []
1387     id_matrix = identity_matrix(field, n)
1388     for i in xrange(n):
1389         ei = id_matrix.column(i)
1390         Qi = zero_matrix(field, n)
1391         Qi.set_row(0, ei)
1392         Qi.set_column(0, ei)
1393         Qi += diagonal_matrix(n, [ei[0]]*n)
1394         # The addition of the diagonal matrix adds an extra ei[0] in the
1395         # upper-left corner of the matrix.
1396         Qi[0,0] = Qi[0,0] * ~field(2)
1397         Qs.append(Qi)
1398
1399     # The rank of the spin factor algebra is two, UNLESS we're in a
1400     # one-dimensional ambient space (the rank is bounded by the
1401     # ambient dimension).
1402     return FiniteDimensionalEuclideanJordanAlgebra(field,
1403                                                    Qs,
1404                                                    rank=min(n,2),
1405                                                    inner_product=_usual_ip)