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