mjo/eja/euclidean_jordan_algebra.py

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