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