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: J = eja_ln(5)
 388                 sage: x = J.random_element()
 389                 sage: u = x.subalgebra_generated_by().random_element()
 390                 sage: u.matrix()*u.vector() == (u**2).vector()
 391                 True
 392
 393             """
 394             # First get the subspace spanned by the powers of myself...
 395             V = self.span_of_powers()
 396             F = self.base_ring()
 397
 398             # Now figure out the entries of the right-multiplication
 399             # matrix for the successive basis elements b0, b1,... of
 400             # that subspace.
 401             mats = []
 402             for b_right in V.basis():
 403                 eja_b_right = self.parent()(b_right)
 404                 b_right_rows = []
 405                 # The first row of the right-multiplication matrix by
 406                 # b1 is what we get if we apply that matrix to b1. The
 407                 # second row of the right multiplication matrix by b1
 408                 # is what we get when we apply that matrix to b2...
 409                 #
 410                 # IMPORTANT: this assumes that all vectors are COLUMN
 411                 # vectors, unlike our superclass (which uses row vectors).
 412                 for b_left in V.basis():
 413                     eja_b_left = self.parent()(b_left)
 414                     # Multiply in the original EJA, but then get the
 415                     # coordinates from the subalgebra in terms of its
 416                     # basis.
 417                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 418                     b_right_rows.append(this_row)
 419                 b_right_matrix = matrix(F, b_right_rows)
 420                 mats.append(b_right_matrix)
 421
 422             # It's an algebra of polynomials in one element, and EJAs
 423             # are power-associative.
 424             #
 425             # TODO: choose generator names intelligently.
 426             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 427
 428
 429         def subalgebra_idempotent(self):
 430             """
 431             Find an idempotent in the associative subalgebra I generate
 432             using Proposition 2.3.5 in Baes.
 433
 434             TESTS::
 435
 436                 sage: set_random_seed()
 437                 sage: J = eja_rn(5)
 438                 sage: c = J.random_element().subalgebra_idempotent()
 439                 sage: c^2 == c
 440                 True
 441                 sage: J = eja_ln(5)
 442                 sage: c = J.random_element().subalgebra_idempotent()
 443                 sage: c^2 == c
 444                 True
 445
 446             """
 447             if self.is_nilpotent():
 448                 raise ValueError("this only works with non-nilpotent elements!")
 449
 450             V = self.span_of_powers()
 451             J = self.subalgebra_generated_by()
 452             # Mis-design warning: the basis used for span_of_powers()
 453             # and subalgebra_generated_by() must be the same, and in
 454             # the same order!
 455             u = J(V.coordinates(self.vector()))
 456
 457             # The image of the matrix of left-u^m-multiplication
 458             # will be minimal for some natural number s...
 459             s = 0
 460             minimal_dim = V.dimension()
 461             for i in xrange(1, V.dimension()):
 462                 this_dim = (u**i).matrix().image().dimension()
 463                 if this_dim < minimal_dim:
 464                     minimal_dim = this_dim
 465                     s = i
 466
 467             # Now minimal_matrix should correspond to the smallest
 468             # non-zero subspace in Baes's (or really, Koecher's)
 469             # proposition.
 470             #
 471             # However, we need to restrict the matrix to work on the
 472             # subspace... or do we? Can't we just solve, knowing that
 473             # A(c) = u^(s+1) should have a solution in the big space,
 474             # too?
 475             #
 476             # Beware, solve_right() means that we're using COLUMN vectors.
 477             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 478             u_next = u**(s+1)
 479             A = u_next.matrix()
 480             c_coordinates = A.solve_right(u_next.vector())
 481
 482             # Now c_coordinates is the idempotent we want, but it's in
 483             # the coordinate system of the subalgebra.
 484             #
 485             # We need the basis for J, but as elements of the parent algebra.
 486             #
 487             basis = [self.parent(v) for v in V.basis()]
 488             return self.parent().linear_combination(zip(c_coordinates, basis))
 489
 490
 491         def trace(self):
 492             """
 493             Return my trace, the sum of my eigenvalues.
 494
 495             EXAMPLES::
 496
 497                 sage: J = eja_ln(3)
 498                 sage: e0,e1,e2 = J.gens()
 499                 sage: x = e0 + e1 + e2
 500                 sage: x.trace()
 501                 2
 502
 503             """
 504             cs = self.characteristic_polynomial().coefficients(sparse=False)
 505             if len(cs) >= 2:
 506                 return -1*cs[-2]
 507             else:
 508                 raise ValueError('charpoly had fewer than 2 coefficients')
 509
 510
 511 def eja_rn(dimension, field=QQ):
 512     """
 513     Return the Euclidean Jordan Algebra corresponding to the set
 514     `R^n` under the Hadamard product.
 515
 516     EXAMPLES:
 517
 518     This multiplication table can be verified by hand::
 519
 520         sage: J = eja_rn(3)
 521         sage: e0,e1,e2 = J.gens()
 522         sage: e0*e0
 523         e0
 524         sage: e0*e1
 525         0
 526         sage: e0*e2
 527         0
 528         sage: e1*e1
 529         e1
 530         sage: e1*e2
 531         0
 532         sage: e2*e2
 533         e2
 534
 535     """
 536     # The FiniteDimensionalAlgebra constructor takes a list of
 537     # matrices, the ith representing right multiplication by the ith
 538     # basis element in the vector space. So if e_1 = (1,0,0), then
 539     # right (Hadamard) multiplication of x by e_1 picks out the first
 540     # component of x; and likewise for the ith basis element e_i.
 541     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 542            for i in xrange(dimension) ]
 543
 544     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 545
 546
 547 def eja_ln(dimension, field=QQ):
 548     """
 549     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 550     cone of the given ``dimension``.
 551
 552     EXAMPLES:
 553
 554     This multiplication table can be verified by hand::
 555
 556         sage: J = eja_ln(4)
 557         sage: e0,e1,e2,e3 = J.gens()
 558         sage: e0*e0
 559         e0
 560         sage: e0*e1
 561         e1
 562         sage: e0*e2
 563         e2
 564         sage: e0*e3
 565         e3
 566         sage: e1*e2
 567         0
 568         sage: e1*e3
 569         0
 570         sage: e2*e3
 571         0
 572
 573     In one dimension, this is the reals under multiplication::
 574
 575       sage: J1 = eja_ln(1)
 576       sage: J2 = eja_rn(1)
 577       sage: J1 == J2
 578       True
 579
 580     """
 581     Qs = []
 582     id_matrix = identity_matrix(field,dimension)
 583     for i in xrange(dimension):
 584         ei = id_matrix.column(i)
 585         Qi = zero_matrix(field,dimension)
 586         Qi.set_row(0, ei)
 587         Qi.set_column(0, ei)
 588         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 589         # The addition of the diagonal matrix adds an extra ei[0] in the
 590         # upper-left corner of the matrix.
 591         Qi[0,0] = Qi[0,0] * ~field(2)
 592         Qs.append(Qi)
 593
 594     # The rank of the spin factor algebra is two, UNLESS we're in a
 595     # one-dimensional ambient space (the rank is bounded by the
 596     # ambient dimension).
 597     rank = min(dimension,2)
 598     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)