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 is_regular(self):
 130             """
 131             Return whether or not this is a regular element.
 132
 133             EXAMPLES:
 134
 135             The identity element always has degree one, but any element
 136             linearly-independent from it is regular::
 137
 138                 sage: J = eja_ln(5)
 139                 sage: J.one().is_regular()
 140                 False
 141                 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
 142                 sage: for x in J.gens():
 143                 ....:     (J.one() + x).is_regular()
 144                 False
 145                 True
 146                 True
 147                 True
 148                 True
 149
 150             """
 151             return self.degree() == self.parent().rank()
 152
 153         def span_of_powers(self):
 154             """
 155             Return the vector space spanned by successive powers of
 156             this element.
 157             """
 158             # The dimension of the subalgebra can't be greater than
 159             # the big algebra, so just put everything into a list
 160             # and let span() get rid of the excess.
 161             V = self.vector().parent()
 162             return V.span( (self**d).vector() for d in xrange(V.dimension()) )
 163
 164
 165         def degree(self):
 166             """
 167             Compute the degree of this element the straightforward way
 168             according to the definition; by appending powers to a list
 169             and figuring out its dimension (that is, whether or not
 170             they're linearly dependent).
 171
 172             EXAMPLES::
 173
 174                 sage: J = eja_ln(4)
 175                 sage: J.one().degree()
 176                 1
 177                 sage: e0,e1,e2,e3 = J.gens()
 178                 sage: (e0 - e1).degree()
 179                 2
 180
 181             In the spin factor algebra (of rank two), all elements that
 182             aren't multiples of the identity are regular::
 183
 184                 sage: set_random_seed()
 185                 sage: n = ZZ.random_element(1,10).abs()
 186                 sage: J = eja_ln(n)
 187                 sage: x = J.random_element()
 188                 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
 189                 True
 190
 191             """
 192             return self.span_of_powers().dimension()
 193
 194
 195         def matrix(self):
 196             """
 197             Return the matrix that represents left- (or right-)
 198             multiplication by this element in the parent algebra.
 199
 200             We have to override this because the superclass method
 201             returns a matrix that acts on row vectors (that is, on
 202             the right).
 203             """
 204             fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 205             return fda_elt.matrix().transpose()
 206
 207
 208         def subalgebra_generated_by(self):
 209             """
 210             Return the associative subalgebra of the parent EJA generated
 211             by this element.
 212
 213             TESTS::
 214
 215                 sage: set_random_seed()
 216                 sage: n = ZZ.random_element(1,10).abs()
 217                 sage: J = eja_rn(n)
 218                 sage: x = J.random_element()
 219                 sage: x.subalgebra_generated_by().is_associative()
 220                 True
 221                 sage: J = eja_ln(n)
 222                 sage: x = J.random_element()
 223                 sage: x.subalgebra_generated_by().is_associative()
 224                 True
 225
 226             Squaring in the subalgebra should be the same thing as
 227             squaring in the superalgebra::
 228
 229                 sage: J = eja_ln(5)
 230                 sage: x = J.random_element()
 231                 sage: u = x.subalgebra_generated_by().random_element()
 232                 sage: u.matrix()*u.vector() == (u**2).vector()
 233                 True
 234
 235             """
 236             # First get the subspace spanned by the powers of myself...
 237             V = self.span_of_powers()
 238             F = self.base_ring()
 239
 240             # Now figure out the entries of the right-multiplication
 241             # matrix for the successive basis elements b0, b1,... of
 242             # that subspace.
 243             mats = []
 244             for b_right in V.basis():
 245                 eja_b_right = self.parent()(b_right)
 246                 b_right_rows = []
 247                 # The first row of the right-multiplication matrix by
 248                 # b1 is what we get if we apply that matrix to b1. The
 249                 # second row of the right multiplication matrix by b1
 250                 # is what we get when we apply that matrix to b2...
 251                 #
 252                 # IMPORTANT: this assumes that all vectors are COLUMN
 253                 # vectors, unlike our superclass (which uses row vectors).
 254                 for b_left in V.basis():
 255                     eja_b_left = self.parent()(b_left)
 256                     # Multiply in the original EJA, but then get the
 257                     # coordinates from the subalgebra in terms of its
 258                     # basis.
 259                     this_row = V.coordinates((eja_b_left*eja_b_right).vector())
 260                     b_right_rows.append(this_row)
 261                 b_right_matrix = matrix(F, b_right_rows)
 262                 mats.append(b_right_matrix)
 263
 264             # It's an algebra of polynomials in one element, and EJAs
 265             # are power-associative.
 266             #
 267             # TODO: choose generator names intelligently.
 268             return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
 269
 270
 271         def minimal_polynomial(self):
 272             """
 273             EXAMPLES::
 274
 275                 sage: set_random_seed()
 276                 sage: n = ZZ.random_element(1,10).abs()
 277                 sage: J = eja_rn(n)
 278                 sage: x = J.random_element()
 279                 sage: x.degree() == x.minimal_polynomial().degree()
 280                 True
 281
 282             ::
 283
 284                 sage: set_random_seed()
 285                 sage: n = ZZ.random_element(1,10).abs()
 286                 sage: J = eja_ln(n)
 287                 sage: x = J.random_element()
 288                 sage: x.degree() == x.minimal_polynomial().degree()
 289                 True
 290
 291             The minimal polynomial and the characteristic polynomial coincide
 292             and are known (see Alizadeh, Example 11.11) for all elements of
 293             the spin factor algebra that aren't scalar multiples of the
 294             identity::
 295
 296                 sage: set_random_seed()
 297                 sage: n = ZZ.random_element(2,10).abs()
 298                 sage: J = eja_ln(n)
 299                 sage: y = J.random_element()
 300                 sage: while y == y.coefficient(0)*J.one():
 301                 ....:     y = J.random_element()
 302                 sage: y0 = y.vector()[0]
 303                 sage: y_bar = y.vector()[1:]
 304                 sage: actual = y.minimal_polynomial()
 305                 sage: x = SR.symbol('x', domain='real')
 306                 sage: expected = x^2 - 2*y0*x + (y0^2 - norm(y_bar)^2)
 307                 sage: bool(actual == expected)
 308                 True
 309
 310             """
 311             # The element we're going to call "minimal_polynomial()" on.
 312             # Either myself, interpreted as an element of a finite-
 313             # dimensional algebra, or an element of an associative
 314             # subalgebra.
 315             elt = None
 316
 317             if self.parent().is_associative():
 318                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 319             else:
 320                 V = self.span_of_powers()
 321                 assoc_subalg = self.subalgebra_generated_by()
 322                 # Mis-design warning: the basis used for span_of_powers()
 323                 # and subalgebra_generated_by() must be the same, and in
 324                 # the same order!
 325                 elt = assoc_subalg(V.coordinates(self.vector()))
 326
 327             # Recursive call, but should work since elt lives in an
 328             # associative algebra.
 329             return elt.minimal_polynomial()
 330
 331
 332         def is_nilpotent(self):
 333             """
 334             Return whether or not some power of this element is zero.
 335
 336             The superclass method won't work unless we're in an
 337             associative algebra, and we aren't. However, we generate
 338             an assocoative subalgebra and we're nilpotent there if and
 339             only if we're nilpotent here (probably).
 340
 341             TESTS:
 342
 343             The identity element is never nilpotent::
 344
 345                 sage: set_random_seed()
 346                 sage: n = ZZ.random_element(2,10).abs()
 347                 sage: J = eja_rn(n)
 348                 sage: J.one().is_nilpotent()
 349                 False
 350                 sage: J = eja_ln(n)
 351                 sage: J.one().is_nilpotent()
 352                 False
 353
 354             The additive identity is always nilpotent::
 355
 356                 sage: set_random_seed()
 357                 sage: n = ZZ.random_element(2,10).abs()
 358                 sage: J = eja_rn(n)
 359                 sage: J.zero().is_nilpotent()
 360                 True
 361                 sage: J = eja_ln(n)
 362                 sage: J.zero().is_nilpotent()
 363                 True
 364
 365             """
 366             # The element we're going to call "is_nilpotent()" on.
 367             # Either myself, interpreted as an element of a finite-
 368             # dimensional algebra, or an element of an associative
 369             # subalgebra.
 370             elt = None
 371
 372             if self.parent().is_associative():
 373                 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
 374             else:
 375                 V = self.span_of_powers()
 376                 assoc_subalg = self.subalgebra_generated_by()
 377                 # Mis-design warning: the basis used for span_of_powers()
 378                 # and subalgebra_generated_by() must be the same, and in
 379                 # the same order!
 380                 elt = assoc_subalg(V.coordinates(self.vector()))
 381
 382             # Recursive call, but should work since elt lives in an
 383             # associative algebra.
 384             return elt.is_nilpotent()
 385
 386
 387         def subalgebra_idempotent(self):
 388             """
 389             Find an idempotent in the associative subalgebra I generate
 390             using Proposition 2.3.5 in Baes.
 391
 392             TESTS::
 393
 394                 sage: set_random_seed()
 395                 sage: J = eja_rn(5)
 396                 sage: c = J.random_element().subalgebra_idempotent()
 397                 sage: c^2 == c
 398                 True
 399                 sage: J = eja_ln(5)
 400                 sage: c = J.random_element().subalgebra_idempotent()
 401                 sage: c^2 == c
 402                 True
 403
 404             """
 405             if self.is_nilpotent():
 406                 raise ValueError("this only works with non-nilpotent elements!")
 407
 408             V = self.span_of_powers()
 409             J = self.subalgebra_generated_by()
 410             # Mis-design warning: the basis used for span_of_powers()
 411             # and subalgebra_generated_by() must be the same, and in
 412             # the same order!
 413             u = J(V.coordinates(self.vector()))
 414
 415             # The image of the matrix of left-u^m-multiplication
 416             # will be minimal for some natural number s...
 417             s = 0
 418             minimal_dim = V.dimension()
 419             for i in xrange(1, V.dimension()):
 420                 this_dim = (u**i).matrix().image().dimension()
 421                 if this_dim < minimal_dim:
 422                     minimal_dim = this_dim
 423                     s = i
 424
 425             # Now minimal_matrix should correspond to the smallest
 426             # non-zero subspace in Baes's (or really, Koecher's)
 427             # proposition.
 428             #
 429             # However, we need to restrict the matrix to work on the
 430             # subspace... or do we? Can't we just solve, knowing that
 431             # A(c) = u^(s+1) should have a solution in the big space,
 432             # too?
 433             #
 434             # Beware, solve_right() means that we're using COLUMN vectors.
 435             # Our FiniteDimensionalAlgebraElement superclass uses rows.
 436             u_next = u**(s+1)
 437             A = u_next.matrix()
 438             c_coordinates = A.solve_right(u_next.vector())
 439
 440             # Now c_coordinates is the idempotent we want, but it's in
 441             # the coordinate system of the subalgebra.
 442             #
 443             # We need the basis for J, but as elements of the parent algebra.
 444             #
 445             basis = [self.parent(v) for v in V.basis()]
 446             return self.parent().linear_combination(zip(c_coordinates, basis))
 447
 448
 449
 450         def characteristic_polynomial(self):
 451             return self.matrix().characteristic_polynomial()
 452
 453
 454 def eja_rn(dimension, field=QQ):
 455     """
 456     Return the Euclidean Jordan Algebra corresponding to the set
 457     `R^n` under the Hadamard product.
 458
 459     EXAMPLES:
 460
 461     This multiplication table can be verified by hand::
 462
 463         sage: J = eja_rn(3)
 464         sage: e0,e1,e2 = J.gens()
 465         sage: e0*e0
 466         e0
 467         sage: e0*e1
 468         0
 469         sage: e0*e2
 470         0
 471         sage: e1*e1
 472         e1
 473         sage: e1*e2
 474         0
 475         sage: e2*e2
 476         e2
 477
 478     """
 479     # The FiniteDimensionalAlgebra constructor takes a list of
 480     # matrices, the ith representing right multiplication by the ith
 481     # basis element in the vector space. So if e_1 = (1,0,0), then
 482     # right (Hadamard) multiplication of x by e_1 picks out the first
 483     # component of x; and likewise for the ith basis element e_i.
 484     Qs = [ matrix(field, dimension, dimension, lambda k,j: 1*(k == j == i))
 485            for i in xrange(dimension) ]
 486
 487     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=dimension)
 488
 489
 490 def eja_ln(dimension, field=QQ):
 491     """
 492     Return the Jordan algebra corresponding to the Lorentz "ice cream"
 493     cone of the given ``dimension``.
 494
 495     EXAMPLES:
 496
 497     This multiplication table can be verified by hand::
 498
 499         sage: J = eja_ln(4)
 500         sage: e0,e1,e2,e3 = J.gens()
 501         sage: e0*e0
 502         e0
 503         sage: e0*e1
 504         e1
 505         sage: e0*e2
 506         e2
 507         sage: e0*e3
 508         e3
 509         sage: e1*e2
 510         0
 511         sage: e1*e3
 512         0
 513         sage: e2*e3
 514         0
 515
 516     In one dimension, this is the reals under multiplication::
 517
 518       sage: J1 = eja_ln(1)
 519       sage: J2 = eja_rn(1)
 520       sage: J1 == J2
 521       True
 522
 523     """
 524     Qs = []
 525     id_matrix = identity_matrix(field,dimension)
 526     for i in xrange(dimension):
 527         ei = id_matrix.column(i)
 528         Qi = zero_matrix(field,dimension)
 529         Qi.set_row(0, ei)
 530         Qi.set_column(0, ei)
 531         Qi += diagonal_matrix(dimension, [ei[0]]*dimension)
 532         # The addition of the diagonal matrix adds an extra ei[0] in the
 533         # upper-left corner of the matrix.
 534         Qi[0,0] = Qi[0,0] * ~field(2)
 535         Qs.append(Qi)
 536
 537     # The rank of the spin factor algebra is two, UNLESS we're in a
 538     # one-dimensional ambient space (the rank is bounded by the
 539     # ambient dimension).
 540     rank = min(dimension,2)
 541     return FiniteDimensionalEuclideanJordanAlgebra(field,Qs,rank=rank)