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 quadratic_representation(self):
 327             """
 328             Return the quadratic representation of this element.
 329
 330             EXAMPLES:
 331
 332             The explicit form in the spin factor algebra is given by
 333             Alizadeh's Example 11.12::
 334
 335                 sage: n = ZZ.random_element(1,10).abs()
 336                 sage: J = eja_ln(n)
 337                 sage: x = J.random_element()
 338                 sage: x_vec = x.vector()
 339                 sage: x0 = x_vec[0]
 340                 sage: x_bar = x_vec[1:]
 341                 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
 342                 sage: B = 2*x0*x_bar.row()
 343                 sage: C = 2*x0*x_bar.column()
 344                 sage: D = identity_matrix(QQ, n-1)
 345                 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
 346                 sage: D = D + 2*x_bar.tensor_product(x_bar)
 347                 sage: Q = block_matrix(2,2,[A,B,C,D])
 348                 sage: Q == x.quadratic_representation()
 349                 True
 350
 351             """
 352             return 2*(self.matrix()**2) - (self**2).matrix()
 353
 354
 355         def span_of_powers(self):
 356             """
 357             Return the vector space spanned by successive powers of
 358             this element.
 359             """
 360             # The dimension of the subalgebra can't be greater than
 361             # the big algebra, so just put everything into a list
 362             # and let span() get rid of the excess.
 363             V = self.vector().parent()
 364             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 365
 366
 367         def subalgebra_generated_by(self):
 368             """
 369             Return the associative subalgebra of the parent EJA generated
 370             by this element.
 371
 372             TESTS::
 373
 374                 sage: set_random_seed()
 375                 sage: x = random_eja().random_element()
 376                 sage: x.subalgebra_generated_by().is_associative()
 377                 True
 378
 379             Squaring in the subalgebra should be the same thing as
 380             squaring in the superalgebra::
 381
 382                 sage: set_random_seed()
 383                 sage: x = random_eja().random_element()
 384                 sage: u = x.subalgebra_generated_by().random_element()
 385                 sage: u.matrix()*u.vector() == (u**2).vector()
 386                 True
 387
 388             """
 389             # First get the subspace spanned by the powers of myself...
 390             V = self.span_of_powers()
 391             F = self.base_ring()
 392
 393             # Now figure out the entries of the right-multiplication
 394             # matrix for the successive basis elements b0, b1,... of
 395             # that subspace.
 396             mats = []
 397             for b_right in V.basis():
 398                 eja_b_right = self.parent()(b_right)
 399                 b_right_rows = []
 400                 # The first row of the right-multiplication matrix by
 401                 # b1 is what we get if we apply that matrix to b1. The
 402                 # second row of the right multiplication matrix by b1
 403                 # is what we get when we apply that matrix to b2...
 404                 #
 405                 # IMPORTANT: this assumes that all vectors are COLUMN
 406                 # vectors, unlike our superclass (which uses row vectors).
 407                 for b_left in V.basis():
 408                     eja_b_left = self.parent()(b_left)
 409                     # Multiply in the original EJA, but then get the
 410                     # coordinates from the subalgebra in terms of its
 411                     # basis.
 412                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 413                     b_right_rows.append(this_row)
 414                 b_right_matrix = matrix(F, b_right_rows)
 415                 mats.append(b_right_matrix)
 416
 417             # It's an algebra of polynomials in one element, and EJAs
 418             # are power-associative.
 419             #
 420             # TODO: choose generator names intelligently.
 421             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 422
 423
 424         def subalgebra_idempotent(self):
 425             """
 426             Find an idempotent in the associative subalgebra I generate
 427             using Proposition 2.3.5 in Baes.
 428
 429             TESTS::
 430
 431                 sage: set_random_seed()
 432                 sage: J = eja_rn(5)
 433                 sage: c = J.random_element().subalgebra_idempotent()
 434                 sage: c^2 == c
 435                 True
 436                 sage: J = eja_ln(5)
 437                 sage: c = J.random_element().subalgebra_idempotent()
 438                 sage: c^2 == c
 439                 True
 440
 441             """
 442             if self.is_nilpotent():
 443                 raise ValueError("this only works with non-nilpotent elements!")
 444
 445             V = self.span_of_powers()
 446             J = self.subalgebra_generated_by()
 447             # Mis-design warning: the basis used for span_of_powers()
 448             # and subalgebra_generated_by() must be the same, and in
 449             # the same order!
 450             u = J(V.coordinates(self.vector()))
 451
 452             # The image of the matrix of left-u^m-multiplication
 453             # will be minimal for some natural number s...
 454             s = 0
 455             minimal_dim = V.dimension()
 456             for i in xrange(1, V.dimension()):
 457                 this_dim = (u**i).matrix().image().dimension()
 458                 if this_dim < minimal_dim:
 459                     minimal_dim = this_dim
 460                     s = i
 461
 462             # Now minimal_matrix should correspond to the smallest
 463             # non-zero subspace in Baes's (or really, Koecher's)
 464             # proposition.
 465             #
 466             # However, we need to restrict the matrix to work on the
 467             # subspace... or do we? Can't we just solve, knowing that
 468             # A(c) = u^(s+1) should have a solution in the big space,
 469             # too?
 470             #
 471             # Beware, solve_right() means that we're using COLUMN vectors.
 472             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 473             u_next = u**(s+1)
 474             A = u_next.matrix()
 475             c_coordinates = A.solve_right(u_next.vector())
 476
 477             # Now c_coordinates is the idempotent we want, but it's in
 478             # the coordinate system of the subalgebra.
 479             #
 480             # We need the basis for J, but as elements of the parent algebra.
 481             #
 482             basis = [self.parent(v) for v in V.basis()]
 483             return self.parent().linear_combination(zip(c_coordinates, basis))
 484
 485
 486         def trace(self):
 487             """
 488             Return my trace, the sum of my eigenvalues.
 489
 490             EXAMPLES::
 491
 492                 sage: J = eja_ln(3)
 493                 sage: e0,e1,e2 = J.gens()
 494                 sage: x = e0 + e1 + e2
 495                 sage: x.trace()
 496                 2
 497
 498             """
 499             cs = self.characteristic_polynomial().coefficients(sparse=False)
 500             if len(cs) >= 2:
 501                 return -1*cs[-2]
 502             else:
 503                 raise ValueError('charpoly had fewer than 2 coefficients')
 504
 505
 506 def eja_rn(dimension, field=QQ):
 507     """
 508     Return the Euclidean Jordan Algebra corresponding to the set
 509     `R^n` under the Hadamard product.
 510
 511     EXAMPLES:
 512
 513     This multiplication table can be verified by hand::
 514
 515         sage: J = eja_rn(3)
 516         sage: e0,e1,e2 = J.gens()
 517         sage: e0*e0
 518         e0
 519         sage: e0*e1
 520         0
 521         sage: e0*e2
 522         0
 523         sage: e1*e1
 524         e1
 525         sage: e1*e2
 526         0
 527         sage: e2*e2
 528         e2
 529
 530     """
 531     # The FiniteDimensionalAlgebra constructor takes a list of
 532     # matrices, the ith representing right multiplication by the ith
 533     # basis element in the vector space. So if e_1 = (1,0,0), then
 534     # right (Hadamard) multiplication of x by e_1 picks out the first
 535     # component of x; and likewise for the ith basis element e_i.
 536     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 537            for i in xrange(dimension) ]
 538
 539     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 540
 541
 542 def eja_ln(dimension, field=QQ):
 543     """
 544     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 545     cone of the given ``dimension``.
 546
 547     EXAMPLES:
 548
 549     This multiplication table can be verified by hand::
 550
 551         sage: J = eja_ln(4)
 552         sage: e0,e1,e2,e3 = J.gens()
 553         sage: e0*e0
 554         e0
 555         sage: e0*e1
 556         e1
 557         sage: e0*e2
 558         e2
 559         sage: e0*e3
 560         e3
 561         sage: e1*e2
 562         0
 563         sage: e1*e3
 564         0
 565         sage: e2*e3
 566         0
 567
 568     In one dimension, this is the reals under multiplication::
 569
 570       sage: J1 = eja_ln(1)
 571       sage: J2 = eja_rn(1)
 572       sage: J1 == J2
 573       True
 574
 575     """
 576     Qs = []
 577     id_matrix = identity_matrix(field,dimension)
 578     for i in xrange(dimension):
 579         ei = id_matrix.column(i)
 580         Qi = zero_matrix(field,dimension)
 581         Qi.set_row(0, ei)
 582         Qi.set_column(0, ei)
 583         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 584         # The addition of the diagonal matrix adds an extra ei[0] in the
 585         # upper-left corner of the matrix.
 586         Qi[0,0] = Qi[0,0] * ~field(2)
 587         Qs.append(Qi)
 588
 589     # The rank of the spin factor algebra is two, UNLESS we're in a
 590     # one-dimensional ambient space (the rank is bounded by the
 591     # ambient dimension).
 592     rank = min(dimension,2)
 593     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)
 594
 595
 596 def eja_sn(dimension, field=QQ):
 597     """
 598     Return the simple Jordan algebra of ``dimension``-by-``dimension``
 599     symmetric matrices over ``field``.
 600
 601     EXAMPLES::
 602
 603         sage: J = eja_sn(2)
 604         sage: e0, e1, e2 = J.gens()
 605         sage: e0*e0
 606         e0
 607         sage: e1*e1
 608         e0 + e2
 609         sage: e2*e2
 610         e2
 611
 612     """
 613     Qs = []
 614
 615     # In S^2, for example, we nominally have four coordinates even
 616     # though the space is of dimension three only. The vector space V
 617     # is supposed to hold the entire long vector, and the subspace W
 618     # of V will be spanned by the vectors that arise from symmetric
 619     # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
 620     V = VectorSpace(field, dimension**2)
 621
 622     # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
 623     # coordinates.
 624     S = []
 625
 626     for i in xrange(dimension):
 627         for j in xrange(i+1):
 628             Eij = matrix(field, dimension, lambda k,l: k==i and l==j)
 629             if i == j:
 630                 Sij = Eij
 631             else:
 632                 Sij = Eij + Eij.transpose()
 633             S.append(Sij)
 634
 635     def mat2vec(m):
 636         return vector(field, m.list())
 637
 638     def vec2mat(v):
 639         return matrix(field, dimension, v.list())
 640
 641     W = V.span( mat2vec(s) for s in S )
 642
 643     # Taking the span above reorders our basis (thanks, jerk!) so we
 644     # need to put our "matrix basis" in the same order as the
 645     # (reordered) vector basis.
 646     S = [ vec2mat(b) for b in W.basis() ]
 647
 648     for s in S:
 649         # Brute force the multiplication-by-s matrix by looping
 650         # through all elements of the basis and doing the computation
 651         # to find out what the corresponding row should be. BEWARE:
 652         # these multiplication tables won't be symmetric! It therefore
 653         # becomes REALLY IMPORTANT that the underlying algebra
 654         # constructor uses ROW vectors and not COLUMN vectors. That's
 655         # why we're computing rows here and not columns.
 656         Q_rows = []
 657         for t in S:
 658             this_row = mat2vec((s*t + t*s)/2)
 659             Q_rows.append(W.coordinates(this_row))
 660         Q = matrix(field,Q_rows)
 661         Qs.append(Q)
 662
 663     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 664
 665
 666 def random_eja():
 667     """
 668     Return a "random" finite-dimensional Euclidean Jordan Algebra.
 669
 670     ALGORITHM:
 671
 672     For now, we choose a random natural number ``n`` (greater than zero)
 673     and then give you back one of the following:
 674
 675       * The cartesian product of the rational numbers ``n`` times; this is
 676         ``QQ^n`` with the Hadamard product.
 677
 678       * The Jordan spin algebra on ``QQ^n``.
 679
 680       * The ``n``-by-``n`` rational symmetric matrices with the symmetric
 681         product.
 682
 683     Later this might be extended to return Cartesian products of the
 684     EJAs above.
 685
 686     TESTS::
 687
 688         sage: random_eja()
 689         Euclidean Jordan algebra of degree...
 690
 691     """
 692     n = ZZ.random_element(1,10).abs()
 693     constructor = choice([eja_rn, eja_ln, eja_sn])
 694     return constructor(dimension=n, field=QQ)