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