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