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