]> gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/eja/euclidean_jordan_algebra.py
projects / sage.d.git / blob
? search:
summary | shortlog | log | commit | commitdiff | tree
history | raw | HEAD
eja: implement element trace in terms of charpoly coefficients.
[sage.d.git] / 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 natural_basis=None):
25 n = len(mult_table)
26 mult_table = [b.base_extend(field) for b in mult_table]
27 for b in mult_table:
28 b.set_immutable()
29 if not (is_Matrix(b) and b.dimensions() == (n, n)):
30 raise ValueError("input is not a multiplication table")
31 mult_table = tuple(mult_table)
32
33 cat = MagmaticAlgebras(field).FiniteDimensional().WithBasis()
34 cat.or_subcategory(category)
35 if assume_associative:
36 cat = cat.Associative()
37
38 names = normalize_names(n, names)
39
40 fda = super(FiniteDimensionalEuclideanJordanAlgebra, cls)
41 return fda.__classcall__(cls,
42 field,
43 mult_table,
44 assume_associative=assume_associative,
45 names=names,
46 category=cat,
47 rank=rank,
48 natural_basis=natural_basis)
49
50
51 def __init__(self,
52 field,
53 mult_table,
54 names='e',
55 assume_associative=False,
56 category=None,
57 rank=None,
58 natural_basis=None):
59 """
60 EXAMPLES:
61
62 By definition, Jordan multiplication commutes::
63
64 sage: set_random_seed()
65 sage: J = random_eja()
66 sage: x = J.random_element()
67 sage: y = J.random_element()
68 sage: x*y == y*x
69 True
70
71 """
72 self._rank = rank
73 self._natural_basis = natural_basis
74 self._multiplication_table = mult_table
75 fda = super(FiniteDimensionalEuclideanJordanAlgebra, self)
76 fda.__init__(field,
77 mult_table,
78 names=names,
79 category=category)
80
81
82 def _repr_(self):
83 """
84 Return a string representation of ``self``.
85 """
86 fmt = "Euclidean Jordan algebra of degree {} over {}"
87 return fmt.format(self.degree(), self.base_ring())
88
89
90
91 @cached_method
92 def _charpoly_coeff(self, i):
93 """
94 Return the coefficient polynomial "a_{i}" of this algebra's
95 general characteristic polynomial.
96
97 Having this be a separate cached method lets us compute and
98 store the trace/determinant (a_{r-1} and a_{0} respectively)
99 separate from the entire characteristic polynomial.
100 """
101 (A_of_x, x) = self._charpoly_matrix()
102 R = A_of_x.base_ring()
103 A_cols = A_of_x.columns()
104 A_cols[i] = (x**self.rank()).vector()
105 numerator = column_matrix(A_of_x.base_ring(), A_cols).det()
106 denominator = A_of_x.det()
107
108 # We're relying on the theory here to ensure that each a_i is
109 # indeed back in R, and the added negative signs are to make
110 # the whole charpoly expression sum to zero.
111 return R(-numerator/denominator)
112
113
114 @cached_method
115 def _charpoly_matrix(self):
116 """
117 Compute the matrix whose entries A_ij are polynomials in
118 X1,...,XN. This same matrix is used in more than one method and
119 it's not so fast to construct.
120 """
121 r = self.rank()
122 n = self.dimension()
123
124 # Construct a new algebra over a multivariate polynomial ring...
125 names = ['X' + str(i) for i in range(1,n+1)]
126 R = PolynomialRing(self.base_ring(), names)
127 J = FiniteDimensionalEuclideanJordanAlgebra(R,
128 self._multiplication_table,
129 rank=r)
130
131 idmat = identity_matrix(J.base_ring(), n)
132
133 x = J(vector(R, R.gens()))
134 l1 = [column_matrix((x**k).vector()) for k in range(r)]
135 l2 = [idmat.column(k-1).column() for k in range(r+1, n+1)]
136 A_of_x = block_matrix(R, 1, n, (l1 + l2))
137 return (A_of_x, x)
138
139
140 @cached_method
141 def characteristic_polynomial(self):
142 """
143 EXAMPLES:
144
145 The characteristic polynomial in the spin algebra is given in
146 Alizadeh, Example 11.11::
147
148 sage: J = JordanSpinEJA(3)
149 sage: p = J.characteristic_polynomial(); p
150 X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2
151 sage: xvec = J.one().vector()
152 sage: p(*xvec)
153 t^2 - 2*t + 1
154
155 """
156 r = self.rank()
157 n = self.dimension()
158
159 # The list of coefficient polynomials a_1, a_2, ..., a_n.
160 a = [ self._charpoly_coeff(i) for i in range(n) ]
161
162 # We go to a bit of trouble here to reorder the
163 # indeterminates, so that it's easier to evaluate the
164 # characteristic polynomial at x's coordinates and get back
165 # something in terms of t, which is what we want.
166 R = a[0].parent()
167 S = PolynomialRing(self.base_ring(),'t')
168 t = S.gen(0)
169 S = PolynomialRing(S, R.variable_names())
170 t = S(t)
171
172 # Note: all entries past the rth should be zero. The
173 # coefficient of the highest power (x^r) is 1, but it doesn't
174 # appear in the solution vector which contains coefficients
175 # for the other powers (to make them sum to x^r).
176 if (r < n):
177 a[r] = 1 # corresponds to x^r
178 else:
179 # When the rank is equal to the dimension, trying to
180 # assign a[r] goes out-of-bounds.
181 a.append(1) # corresponds to x^r
182
183 return sum( a[k]*(t**k) for k in range(len(a)) )
184
185
186 def inner_product(self, x, y):
187 """
188 The inner product associated with this Euclidean Jordan algebra.
189
190 Defaults to the trace inner product, but can be overridden by
191 subclasses if they are sure that the necessary properties are
192 satisfied.
193
194 EXAMPLES:
195
196 The inner product must satisfy its axiom for this algebra to truly
197 be a Euclidean Jordan Algebra::
198
199 sage: set_random_seed()
200 sage: J = random_eja()
201 sage: x = J.random_element()
202 sage: y = J.random_element()
203 sage: z = J.random_element()
204 sage: (x*y).inner_product(z) == y.inner_product(x*z)
205 True
206
207 """
208 if (not x in self) or (not y in self):
209 raise TypeError("arguments must live in this algebra")
210 return x.trace_inner_product(y)
211
212
213 def natural_basis(self):
214 """
215 Return a more-natural representation of this algebra's basis.
216
217 Every finite-dimensional Euclidean Jordan Algebra is a direct
218 sum of five simple algebras, four of which comprise Hermitian
219 matrices. This method returns the original "natural" basis
220 for our underlying vector space. (Typically, the natural basis
221 is used to construct the multiplication table in the first place.)
222
223 Note that this will always return a matrix. The standard basis
224 in `R^n` will be returned as `n`-by-`1` column matrices.
225
226 EXAMPLES::
227
228 sage: J = RealSymmetricEJA(2)
229 sage: J.basis()
230 Family (e0, e1, e2)
231 sage: J.natural_basis()
232 (
233 [1 0] [0 1] [0 0]
234 [0 0], [1 0], [0 1]
235 )
236
237 ::
238
239 sage: J = JordanSpinEJA(2)
240 sage: J.basis()
241 Family (e0, e1)
242 sage: J.natural_basis()
243 (
244 [1] [0]
245 [0], [1]
246 )
247
248 """
249 if self._natural_basis is None:
250 return tuple( b.vector().column() for b in self.basis() )
251 else:
252 return self._natural_basis
253
254
255 def rank(self):
256 """
257 Return the rank of this EJA.
258 """
259 if self._rank is None:
260 raise ValueError("no rank specified at genesis")
261 else:
262 return self._rank
263
264
265 class Element(FiniteDimensionalAlgebraElement):
266 """
267 An element of a Euclidean Jordan algebra.
268 """
269
270 def __init__(self, A, elt=None):
271 """
272 EXAMPLES:
273
274 The identity in `S^n` is converted to the identity in the EJA::
275
276 sage: J = RealSymmetricEJA(3)
277 sage: I = identity_matrix(QQ,3)
278 sage: J(I) == J.one()
279 True
280
281 This skew-symmetric matrix can't be represented in the EJA::
282
283 sage: J = RealSymmetricEJA(3)
284 sage: A = matrix(QQ,3, lambda i,j: i-j)
285 sage: J(A)
286 Traceback (most recent call last):
287 ...
288 ArithmeticError: vector is not in free module
289
290 """
291 # Goal: if we're given a matrix, and if it lives in our
292 # parent algebra's "natural ambient space," convert it
293 # into an algebra element.
294 #
295 # The catch is, we make a recursive call after converting
296 # the given matrix into a vector that lives in the algebra.
297 # This we need to try the parent class initializer first,
298 # to avoid recursing forever if we're given something that
299 # already fits into the algebra, but also happens to live
300 # in the parent's "natural ambient space" (this happens with
301 # vectors in R^n).
302 try:
303 FiniteDimensionalAlgebraElement.__init__(self, A, elt)
304 except ValueError:
305 natural_basis = A.natural_basis()
306 if elt in natural_basis[0].matrix_space():
307 # Thanks for nothing! Matrix spaces aren't vector
308 # spaces in Sage, so we have to figure out its
309 # natural-basis coordinates ourselves.
310 V = VectorSpace(elt.base_ring(), elt.nrows()**2)
311 W = V.span( _mat2vec(s) for s in natural_basis )
312 coords = W.coordinates(_mat2vec(elt))
313 FiniteDimensionalAlgebraElement.__init__(self, A, coords)
314
315 def __pow__(self, n):
316 """
317 Return ``self`` raised to the power ``n``.
318
319 Jordan algebras are always power-associative; see for
320 example Faraut and Koranyi, Proposition II.1.2 (ii).
321
322 .. WARNING:
323
324 We have to override this because our superclass uses row vectors
325 instead of column vectors! We, on the other hand, assume column
326 vectors everywhere.
327
328 EXAMPLES::
329
330 sage: set_random_seed()
331 sage: x = random_eja().random_element()
332 sage: x.operator_matrix()*x.vector() == (x^2).vector()
333 True
334
335 A few examples of power-associativity::
336
337 sage: set_random_seed()
338 sage: x = random_eja().random_element()
339 sage: x*(x*x)*(x*x) == x^5
340 True
341 sage: (x*x)*(x*x*x) == x^5
342 True
343
344 We also know that powers operator-commute (Koecher, Chapter
345 III, Corollary 1)::
346
347 sage: set_random_seed()
348 sage: x = random_eja().random_element()
349 sage: m = ZZ.random_element(0,10)
350 sage: n = ZZ.random_element(0,10)
351 sage: Lxm = (x^m).operator_matrix()
352 sage: Lxn = (x^n).operator_matrix()
353 sage: Lxm*Lxn == Lxn*Lxm
354 True
355
356 """
357 A = self.parent()
358 if n == 0:
359 return A.one()
360 elif n == 1:
361 return self
362 else:
363 return A( (self.operator_matrix()**(n-1))*self.vector() )
364
365
366 def apply_univariate_polynomial(self, p):
367 """
368 Apply the univariate polynomial ``p`` to this element.
369
370 A priori, SageMath won't allow us to apply a univariate
371 polynomial to an element of an EJA, because we don't know
372 that EJAs are rings (they are usually not associative). Of
373 course, we know that EJAs are power-associative, so the
374 operation is ultimately kosher. This function sidesteps
375 the CAS to get the answer we want and expect.
376
377 EXAMPLES::
378
379 sage: R = PolynomialRing(QQ, 't')
380 sage: t = R.gen(0)
381 sage: p = t^4 - t^3 + 5*t - 2
382 sage: J = RealCartesianProductEJA(5)
383 sage: J.one().apply_univariate_polynomial(p) == 3*J.one()
384 True
385
386 TESTS:
387
388 We should always get back an element of the algebra::
389
390 sage: set_random_seed()
391 sage: p = PolynomialRing(QQ, 't').random_element()
392 sage: J = random_eja()
393 sage: x = J.random_element()
394 sage: x.apply_univariate_polynomial(p) in J
395 True
396
397 """
398 if len(p.variables()) > 1:
399 raise ValueError("not a univariate polynomial")
400 P = self.parent()
401 R = P.base_ring()
402 # Convert the coeficcients to the parent's base ring,
403 # because a priori they might live in an (unnecessarily)
404 # larger ring for which P.sum() would fail below.
405 cs = [ R(c) for c in p.coefficients(sparse=False) ]
406 return P.sum( cs[k]*(self**k) for k in range(len(cs)) )
407
408
409 def characteristic_polynomial(self):
410 """
411 Return the characteristic polynomial of this element.
412
413 EXAMPLES:
414
415 The rank of `R^3` is three, and the minimal polynomial of
416 the identity element is `(t-1)` from which it follows that
417 the characteristic polynomial should be `(t-1)^3`::
418
419 sage: J = RealCartesianProductEJA(3)
420 sage: J.one().characteristic_polynomial()
421 t^3 - 3*t^2 + 3*t - 1
422
423 Likewise, the characteristic of the zero element in the
424 rank-three algebra `R^{n}` should be `t^{3}`::
425
426 sage: J = RealCartesianProductEJA(3)
427 sage: J.zero().characteristic_polynomial()
428 t^3
429
430 The characteristic polynomial of an element should evaluate
431 to zero on that element::
432
433 sage: set_random_seed()
434 sage: x = RealCartesianProductEJA(3).random_element()
435 sage: p = x.characteristic_polynomial()
436 sage: x.apply_univariate_polynomial(p)
437 0
438
439 """
440 p = self.parent().characteristic_polynomial()
441 return p(*self.vector())
442
443
444 def inner_product(self, other):
445 """
446 Return the parent algebra's inner product of myself and ``other``.
447
448 EXAMPLES:
449
450 The inner product in the Jordan spin algebra is the usual
451 inner product on `R^n` (this example only works because the
452 basis for the Jordan algebra is the standard basis in `R^n`)::
453
454 sage: J = JordanSpinEJA(3)
455 sage: x = vector(QQ,[1,2,3])
456 sage: y = vector(QQ,[4,5,6])
457 sage: x.inner_product(y)
458 32
459 sage: J(x).inner_product(J(y))
460 32
461
462 The inner product on `S^n` is `<X,Y> = trace(X*Y)`, where
463 multiplication is the usual matrix multiplication in `S^n`,
464 so the inner product of the identity matrix with itself
465 should be the `n`::
466
467 sage: J = RealSymmetricEJA(3)
468 sage: J.one().inner_product(J.one())
469 3
470
471 Likewise, the inner product on `C^n` is `<X,Y> =
472 Re(trace(X*Y))`, where we must necessarily take the real
473 part because the product of Hermitian matrices may not be
474 Hermitian::
475
476 sage: J = ComplexHermitianEJA(3)
477 sage: J.one().inner_product(J.one())
478 3
479
480 Ditto for the quaternions::
481
482 sage: J = QuaternionHermitianEJA(3)
483 sage: J.one().inner_product(J.one())
484 3
485
486 TESTS:
487
488 Ensure that we can always compute an inner product, and that
489 it gives us back a real number::
490
491 sage: set_random_seed()
492 sage: J = random_eja()
493 sage: x = J.random_element()
494 sage: y = J.random_element()
495 sage: x.inner_product(y) in RR
496 True
497
498 """
499 P = self.parent()
500 if not other in P:
501 raise TypeError("'other' must live in the same algebra")
502
503 return P.inner_product(self, other)
504
505
506 def operator_commutes_with(self, other):
507 """
508 Return whether or not this element operator-commutes
509 with ``other``.
510
511 EXAMPLES:
512
513 The definition of a Jordan algebra says that any element
514 operator-commutes with its square::
515
516 sage: set_random_seed()
517 sage: x = random_eja().random_element()
518 sage: x.operator_commutes_with(x^2)
519 True
520
521 TESTS:
522
523 Test Lemma 1 from Chapter III of Koecher::
524
525 sage: set_random_seed()
526 sage: J = random_eja()
527 sage: u = J.random_element()
528 sage: v = J.random_element()
529 sage: lhs = u.operator_commutes_with(u*v)
530 sage: rhs = v.operator_commutes_with(u^2)
531 sage: lhs == rhs
532 True
533
534 """
535 if not other in self.parent():
536 raise TypeError("'other' must live in the same algebra")
537
538 A = self.operator_matrix()
539 B = other.operator_matrix()
540 return (A*B == B*A)
541
542
543 def det(self):
544 """
545 Return my determinant, the product of my eigenvalues.
546
547 EXAMPLES::
548
549 sage: J = JordanSpinEJA(2)
550 sage: e0,e1 = J.gens()
551 sage: x = e0 + e1
552 sage: x.det()
553 0
554 sage: J = JordanSpinEJA(3)
555 sage: e0,e1,e2 = J.gens()
556 sage: x = e0 + e1 + e2
557 sage: x.det()
558 -1
559
560 """
561 cs = self.characteristic_polynomial().coefficients(sparse=False)
562 r = len(cs) - 1
563 if r >= 0:
564 return cs[0] * (-1)**r
565 else:
566 raise ValueError('charpoly had no coefficients')
567
568
569 def inverse(self):
570 """
571 Return the Jordan-multiplicative inverse of this element.
572
573 We can't use the superclass method because it relies on the
574 algebra being associative.
575
576 EXAMPLES:
577
578 The inverse in the spin factor algebra is given in Alizadeh's
579 Example 11.11::
580
581 sage: set_random_seed()
582 sage: n = ZZ.random_element(1,10)
583 sage: J = JordanSpinEJA(n)
584 sage: x = J.random_element()
585 sage: while x.is_zero():
586 ....: x = J.random_element()
587 sage: x_vec = x.vector()
588 sage: x0 = x_vec[0]
589 sage: x_bar = x_vec[1:]
590 sage: coeff = 1/(x0^2 - x_bar.inner_product(x_bar))
591 sage: inv_vec = x_vec.parent()([x0] + (-x_bar).list())
592 sage: x_inverse = coeff*inv_vec
593 sage: x.inverse() == J(x_inverse)
594 True
595
596 TESTS:
597
598 The identity element is its own inverse::
599
600 sage: set_random_seed()
601 sage: J = random_eja()
602 sage: J.one().inverse() == J.one()
603 True
604
605 If an element has an inverse, it acts like one. TODO: this
606 can be a lot less ugly once ``is_invertible`` doesn't crash
607 on irregular elements::
608
609 sage: set_random_seed()
610 sage: J = random_eja()
611 sage: x = J.random_element()
612 sage: try:
613 ....: x.inverse()*x == J.one()
614 ....: except:
615 ....: True
616 True
617
618 """
619 if self.parent().is_associative():
620 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
621 return elt.inverse()
622
623 # TODO: we can do better once the call to is_invertible()
624 # doesn't crash on irregular elements.
625 #if not self.is_invertible():
626 # raise ValueError('element is not invertible')
627
628 # We do this a little different than the usual recursive
629 # call to a finite-dimensional algebra element, because we
630 # wind up with an inverse that lives in the subalgebra and
631 # we need information about the parent to convert it back.
632 V = self.span_of_powers()
633 assoc_subalg = self.subalgebra_generated_by()
634 # Mis-design warning: the basis used for span_of_powers()
635 # and subalgebra_generated_by() must be the same, and in
636 # the same order!
637 elt = assoc_subalg(V.coordinates(self.vector()))
638
639 # This will be in the subalgebra's coordinates...
640 fda_elt = FiniteDimensionalAlgebraElement(assoc_subalg, elt)
641 subalg_inverse = fda_elt.inverse()
642
643 # So we have to convert back...
644 basis = [ self.parent(v) for v in V.basis() ]
645 pairs = zip(subalg_inverse.vector(), basis)
646 return self.parent().linear_combination(pairs)
647
648
649 def is_invertible(self):
650 """
651 Return whether or not this element is invertible.
652
653 We can't use the superclass method because it relies on
654 the algebra being associative.
655
656 ALGORITHM:
657
658 The usual way to do this is to check if the determinant is
659 zero, but we need the characteristic polynomial for the
660 determinant. The minimal polynomial is a lot easier to get,
661 so we use Corollary 2 in Chapter V of Koecher to check
662 whether or not the paren't algebra's zero element is a root
663 of this element's minimal polynomial.
664
665 TESTS:
666
667 The identity element is always invertible::
668
669 sage: set_random_seed()
670 sage: J = random_eja()
671 sage: J.one().is_invertible()
672 True
673
674 The zero element is never invertible::
675
676 sage: set_random_seed()
677 sage: J = random_eja()
678 sage: J.zero().is_invertible()
679 False
680
681 """
682 zero = self.parent().zero()
683 p = self.minimal_polynomial()
684 return not (p(zero) == zero)
685
686
687 def is_nilpotent(self):
688 """
689 Return whether or not some power of this element is zero.
690
691 The superclass method won't work unless we're in an
692 associative algebra, and we aren't. However, we generate
693 an assocoative subalgebra and we're nilpotent there if and
694 only if we're nilpotent here (probably).
695
696 TESTS:
697
698 The identity element is never nilpotent::
699
700 sage: set_random_seed()
701 sage: random_eja().one().is_nilpotent()
702 False
703
704 The additive identity is always nilpotent::
705
706 sage: set_random_seed()
707 sage: random_eja().zero().is_nilpotent()
708 True
709
710 """
711 # The element we're going to call "is_nilpotent()" on.
712 # Either myself, interpreted as an element of a finite-
713 # dimensional algebra, or an element of an associative
714 # subalgebra.
715 elt = None
716
717 if self.parent().is_associative():
718 elt = FiniteDimensionalAlgebraElement(self.parent(), self)
719 else:
720 V = self.span_of_powers()
721 assoc_subalg = self.subalgebra_generated_by()
722 # Mis-design warning: the basis used for span_of_powers()
723 # and subalgebra_generated_by() must be the same, and in
724 # the same order!
725 elt = assoc_subalg(V.coordinates(self.vector()))
726
727 # Recursive call, but should work since elt lives in an
728 # associative algebra.
729 return elt.is_nilpotent()
730
731
732 def is_regular(self):
733 """
734 Return whether or not this is a regular element.
735
736 EXAMPLES:
737
738 The identity element always has degree one, but any element
739 linearly-independent from it is regular::
740
741 sage: J = JordanSpinEJA(5)
742 sage: J.one().is_regular()
743 False
744 sage: e0, e1, e2, e3, e4 = J.gens() # e0 is the identity
745 sage: for x in J.gens():
746 ....: (J.one() + x).is_regular()
747 False
748 True
749 True
750 True
751 True
752
753 """
754 return self.degree() == self.parent().rank()
755
756
757 def degree(self):
758 """
759 Compute the degree of this element the straightforward way
760 according to the definition; by appending powers to a list
761 and figuring out its dimension (that is, whether or not
762 they're linearly dependent).
763
764 EXAMPLES::
765
766 sage: J = JordanSpinEJA(4)
767 sage: J.one().degree()
768 1
769 sage: e0,e1,e2,e3 = J.gens()
770 sage: (e0 - e1).degree()
771 2
772
773 In the spin factor algebra (of rank two), all elements that
774 aren't multiples of the identity are regular::
775
776 sage: set_random_seed()
777 sage: n = ZZ.random_element(1,10)
778 sage: J = JordanSpinEJA(n)
779 sage: x = J.random_element()
780 sage: x == x.coefficient(0)*J.one() or x.degree() == 2
781 True
782
783 """
784 return self.span_of_powers().dimension()
785
786
787 def minimal_polynomial(self):
788 """
789 Return the minimal polynomial of this element,
790 as a function of the variable `t`.
791
792 ALGORITHM:
793
794 We restrict ourselves to the associative subalgebra
795 generated by this element, and then return the minimal
796 polynomial of this element's operator matrix (in that
797 subalgebra). This works by Baes Proposition 2.3.16.
798
799 TESTS:
800
801 The minimal polynomial of the identity and zero elements are
802 always the same::
803
804 sage: set_random_seed()
805 sage: J = random_eja()
806 sage: J.one().minimal_polynomial()
807 t - 1
808 sage: J.zero().minimal_polynomial()
809 t
810
811 The degree of an element is (by one definition) the degree
812 of its minimal polynomial::
813
814 sage: set_random_seed()
815 sage: x = random_eja().random_element()
816 sage: x.degree() == x.minimal_polynomial().degree()
817 True
818
819 The minimal polynomial and the characteristic polynomial coincide
820 and are known (see Alizadeh, Example 11.11) for all elements of
821 the spin factor algebra that aren't scalar multiples of the
822 identity::
823
824 sage: set_random_seed()
825 sage: n = ZZ.random_element(2,10)
826 sage: J = JordanSpinEJA(n)
827 sage: y = J.random_element()
828 sage: while y == y.coefficient(0)*J.one():
829 ....: y = J.random_element()
830 sage: y0 = y.vector()[0]
831 sage: y_bar = y.vector()[1:]
832 sage: actual = y.minimal_polynomial()
833 sage: t = PolynomialRing(J.base_ring(),'t').gen(0)
834 sage: expected = t^2 - 2*y0*t + (y0^2 - norm(y_bar)^2)
835 sage: bool(actual == expected)
836 True
837
838 The minimal polynomial should always kill its element::
839
840 sage: set_random_seed()
841 sage: x = random_eja().random_element()
842 sage: p = x.minimal_polynomial()
843 sage: x.apply_univariate_polynomial(p)
844 0
845
846 """
847 V = self.span_of_powers()
848 assoc_subalg = self.subalgebra_generated_by()
849 # Mis-design warning: the basis used for span_of_powers()
850 # and subalgebra_generated_by() must be the same, and in
851 # the same order!
852 elt = assoc_subalg(V.coordinates(self.vector()))
853
854 # We get back a symbolic polynomial in 'x' but want a real
855 # polynomial in 't'.
856 p_of_x = elt.operator_matrix().minimal_polynomial()
857 return p_of_x.change_variable_name('t')
858
859
860 def natural_representation(self):
861 """
862 Return a more-natural representation of this element.
863
864 Every finite-dimensional Euclidean Jordan Algebra is a
865 direct sum of five simple algebras, four of which comprise
866 Hermitian matrices. This method returns the original
867 "natural" representation of this element as a Hermitian
868 matrix, if it has one. If not, you get the usual representation.
869
870 EXAMPLES::
871
872 sage: J = ComplexHermitianEJA(3)
873 sage: J.one()
874 e0 + e5 + e8
875 sage: J.one().natural_representation()
876 [1 0 0 0 0 0]
877 [0 1 0 0 0 0]
878 [0 0 1 0 0 0]
879 [0 0 0 1 0 0]
880 [0 0 0 0 1 0]
881 [0 0 0 0 0 1]
882
883 ::
884
885 sage: J = QuaternionHermitianEJA(3)
886 sage: J.one()
887 e0 + e9 + e14
888 sage: J.one().natural_representation()
889 [1 0 0 0 0 0 0 0 0 0 0 0]
890 [0 1 0 0 0 0 0 0 0 0 0 0]
891 [0 0 1 0 0 0 0 0 0 0 0 0]
892 [0 0 0 1 0 0 0 0 0 0 0 0]
893 [0 0 0 0 1 0 0 0 0 0 0 0]
894 [0 0 0 0 0 1 0 0 0 0 0 0]
895 [0 0 0 0 0 0 1 0 0 0 0 0]
896 [0 0 0 0 0 0 0 1 0 0 0 0]
897 [0 0 0 0 0 0 0 0 1 0 0 0]
898 [0 0 0 0 0 0 0 0 0 1 0 0]
899 [0 0 0 0 0 0 0 0 0 0 1 0]
900 [0 0 0 0 0 0 0 0 0 0 0 1]
901
902 """
903 B = self.parent().natural_basis()
904 W = B[0].matrix_space()
905 return W.linear_combination(zip(self.vector(), B))
906
907
908 def operator_matrix(self):
909 """
910 Return the matrix that represents left- (or right-)
911 multiplication by this element in the parent algebra.
912
913 We have to override this because the superclass method
914 returns a matrix that acts on row vectors (that is, on
915 the right).
916
917 EXAMPLES:
918
919 Test the first polarization identity from my notes, Koecher Chapter
920 III, or from Baes (2.3)::
921
922 sage: set_random_seed()
923 sage: J = random_eja()
924 sage: x = J.random_element()
925 sage: y = J.random_element()
926 sage: Lx = x.operator_matrix()
927 sage: Ly = y.operator_matrix()
928 sage: Lxx = (x*x).operator_matrix()
929 sage: Lxy = (x*y).operator_matrix()
930 sage: bool(2*Lx*Lxy + Ly*Lxx == 2*Lxy*Lx + Lxx*Ly)
931 True
932
933 Test the second polarization identity from my notes or from
934 Baes (2.4)::
935
936 sage: set_random_seed()
937 sage: J = random_eja()
938 sage: x = J.random_element()
939 sage: y = J.random_element()
940 sage: z = J.random_element()
941 sage: Lx = x.operator_matrix()
942 sage: Ly = y.operator_matrix()
943 sage: Lz = z.operator_matrix()
944 sage: Lzy = (z*y).operator_matrix()
945 sage: Lxy = (x*y).operator_matrix()
946 sage: Lxz = (x*z).operator_matrix()
947 sage: bool(Lx*Lzy + Lz*Lxy + Ly*Lxz == Lzy*Lx + Lxy*Lz + Lxz*Ly)
948 True
949
950 Test the third polarization identity from my notes or from
951 Baes (2.5)::
952
953 sage: set_random_seed()
954 sage: J = random_eja()
955 sage: u = J.random_element()
956 sage: y = J.random_element()
957 sage: z = J.random_element()
958 sage: Lu = u.operator_matrix()
959 sage: Ly = y.operator_matrix()
960 sage: Lz = z.operator_matrix()
961 sage: Lzy = (z*y).operator_matrix()
962 sage: Luy = (u*y).operator_matrix()
963 sage: Luz = (u*z).operator_matrix()
964 sage: Luyz = (u*(y*z)).operator_matrix()
965 sage: lhs = Lu*Lzy + Lz*Luy + Ly*Luz
966 sage: rhs = Luyz + Ly*Lu*Lz + Lz*Lu*Ly
967 sage: bool(lhs == rhs)
968 True
969
970 """
971 fda_elt = FiniteDimensionalAlgebraElement(self.parent(), self)
972 return fda_elt.matrix().transpose()
973
974
975 def quadratic_representation(self, other=None):
976 """
977 Return the quadratic representation of this element.
978
979 EXAMPLES:
980
981 The explicit form in the spin factor algebra is given by
982 Alizadeh's Example 11.12::
983
984 sage: set_random_seed()
985 sage: n = ZZ.random_element(1,10)
986 sage: J = JordanSpinEJA(n)
987 sage: x = J.random_element()
988 sage: x_vec = x.vector()
989 sage: x0 = x_vec[0]
990 sage: x_bar = x_vec[1:]
991 sage: A = matrix(QQ, 1, [x_vec.inner_product(x_vec)])
992 sage: B = 2*x0*x_bar.row()
993 sage: C = 2*x0*x_bar.column()
994 sage: D = identity_matrix(QQ, n-1)
995 sage: D = (x0^2 - x_bar.inner_product(x_bar))*D
996 sage: D = D + 2*x_bar.tensor_product(x_bar)
997 sage: Q = block_matrix(2,2,[A,B,C,D])
998 sage: Q == x.quadratic_representation()
999 True
1000
1001 Test all of the properties from Theorem 11.2 in Alizadeh::
1002
1003 sage: set_random_seed()
1004 sage: J = random_eja()
1005 sage: x = J.random_element()
1006 sage: y = J.random_element()
1007
1008 Property 1:
1009
1010 sage: actual = x.quadratic_representation(y)
1011 sage: expected = ( (x+y).quadratic_representation()
1012 ....: -x.quadratic_representation()
1013 ....: -y.quadratic_representation() ) / 2
1014 sage: actual == expected
1015 True
1016
1017 Property 2:
1018
1019 sage: alpha = QQ.random_element()
1020 sage: actual = (alpha*x).quadratic_representation()
1021 sage: expected = (alpha^2)*x.quadratic_representation()
1022 sage: actual == expected
1023 True
1024
1025 Property 5:
1026
1027 sage: Qy = y.quadratic_representation()
1028 sage: actual = J(Qy*x.vector()).quadratic_representation()
1029 sage: expected = Qy*x.quadratic_representation()*Qy
1030 sage: actual == expected
1031 True
1032
1033 Property 6:
1034
1035 sage: k = ZZ.random_element(1,10)
1036 sage: actual = (x^k).quadratic_representation()
1037 sage: expected = (x.quadratic_representation())^k
1038 sage: actual == expected
1039 True
1040
1041 """
1042 if other is None:
1043 other=self
1044 elif not other in self.parent():
1045 raise TypeError("'other' must live in the same algebra")
1046
1047 L = self.operator_matrix()
1048 M = other.operator_matrix()
1049 return ( L*M + M*L - (self*other).operator_matrix() )
1050
1051
1052 def span_of_powers(self):
1053 """
1054 Return the vector space spanned by successive powers of
1055 this element.
1056 """
1057 # The dimension of the subalgebra can't be greater than
1058 # the big algebra, so just put everything into a list
1059 # and let span() get rid of the excess.
1060 #
1061 # We do the extra ambient_vector_space() in case we're messing
1062 # with polynomials and the direct parent is a module.
1063 V = self.vector().parent().ambient_vector_space()
1064 return V.span( (self**d).vector() for d in xrange(V.dimension()) )
1065
1066
1067 def subalgebra_generated_by(self):
1068 """
1069 Return the associative subalgebra of the parent EJA generated
1070 by this element.
1071
1072 TESTS::
1073
1074 sage: set_random_seed()
1075 sage: x = random_eja().random_element()
1076 sage: x.subalgebra_generated_by().is_associative()
1077 True
1078
1079 Squaring in the subalgebra should be the same thing as
1080 squaring in the superalgebra::
1081
1082 sage: set_random_seed()
1083 sage: x = random_eja().random_element()
1084 sage: u = x.subalgebra_generated_by().random_element()
1085 sage: u.operator_matrix()*u.vector() == (u**2).vector()
1086 True
1087
1088 """
1089 # First get the subspace spanned by the powers of myself...
1090 V = self.span_of_powers()
1091 F = self.base_ring()
1092
1093 # Now figure out the entries of the right-multiplication
1094 # matrix for the successive basis elements b0, b1,... of
1095 # that subspace.
1096 mats = []
1097 for b_right in V.basis():
1098 eja_b_right = self.parent()(b_right)
1099 b_right_rows = []
1100 # The first row of the right-multiplication matrix by
1101 # b1 is what we get if we apply that matrix to b1. The
1102 # second row of the right multiplication matrix by b1
1103 # is what we get when we apply that matrix to b2...
1104 #
1105 # IMPORTANT: this assumes that all vectors are COLUMN
1106 # vectors, unlike our superclass (which uses row vectors).
1107 for b_left in V.basis():
1108 eja_b_left = self.parent()(b_left)
1109 # Multiply in the original EJA, but then get the
1110 # coordinates from the subalgebra in terms of its
1111 # basis.
1112 this_row = V.coordinates((eja_b_left*eja_b_right).vector())
1113 b_right_rows.append(this_row)
1114 b_right_matrix = matrix(F, b_right_rows)
1115 mats.append(b_right_matrix)
1116
1117 # It's an algebra of polynomials in one element, and EJAs
1118 # are power-associative.
1119 #
1120 # TODO: choose generator names intelligently.
1121 return FiniteDimensionalEuclideanJordanAlgebra(F, mats, assume_associative=True, names='f')
1122
1123
1124 def subalgebra_idempotent(self):
1125 """
1126 Find an idempotent in the associative subalgebra I generate
1127 using Proposition 2.3.5 in Baes.
1128
1129 TESTS::
1130
1131 sage: set_random_seed()
1132 sage: J = RealCartesianProductEJA(5)
1133 sage: c = J.random_element().subalgebra_idempotent()
1134 sage: c^2 == c
1135 True
1136 sage: J = JordanSpinEJA(5)
1137 sage: c = J.random_element().subalgebra_idempotent()
1138 sage: c^2 == c
1139 True
1140
1141 """
1142 if self.is_nilpotent():
1143 raise ValueError("this only works with non-nilpotent elements!")
1144
1145 V = self.span_of_powers()
1146 J = self.subalgebra_generated_by()
1147 # Mis-design warning: the basis used for span_of_powers()
1148 # and subalgebra_generated_by() must be the same, and in
1149 # the same order!
1150 u = J(V.coordinates(self.vector()))
1151
1152 # The image of the matrix of left-u^m-multiplication
1153 # will be minimal for some natural number s...
1154 s = 0
1155 minimal_dim = V.dimension()
1156 for i in xrange(1, V.dimension()):
1157 this_dim = (u**i).operator_matrix().image().dimension()
1158 if this_dim < minimal_dim:
1159 minimal_dim = this_dim
1160 s = i
1161
1162 # Now minimal_matrix should correspond to the smallest
1163 # non-zero subspace in Baes's (or really, Koecher's)
1164 # proposition.
1165 #
1166 # However, we need to restrict the matrix to work on the
1167 # subspace... or do we? Can't we just solve, knowing that
1168 # A(c) = u^(s+1) should have a solution in the big space,
1169 # too?
1170 #
1171 # Beware, solve_right() means that we're using COLUMN vectors.
1172 # Our FiniteDimensionalAlgebraElement superclass uses rows.
1173 u_next = u**(s+1)
1174 A = u_next.operator_matrix()
1175 c_coordinates = A.solve_right(u_next.vector())
1176
1177 # Now c_coordinates is the idempotent we want, but it's in
1178 # the coordinate system of the subalgebra.
1179 #
1180 # We need the basis for J, but as elements of the parent algebra.
1181 #
1182 basis = [self.parent(v) for v in V.basis()]
1183 return self.parent().linear_combination(zip(c_coordinates, basis))
1184
1185
1186 def trace(self):
1187 """
1188 Return my trace, the sum of my eigenvalues.
1189
1190 EXAMPLES::
1191
1192 sage: J = JordanSpinEJA(3)
1193 sage: x = sum(J.gens())
1194 sage: x.trace()
1195 2
1196
1197 ::
1198
1199 sage: J = RealCartesianProductEJA(5)
1200 sage: J.one().trace()
1201 5
1202
1203 TESTS:
1204
1205 The trace of an element is a real number::
1206
1207 sage: set_random_seed()
1208 sage: J = random_eja()
1209 sage: J.random_element().trace() in J.base_ring()
1210 True
1211
1212 """
1213 P = self.parent()
1214 r = P.rank()
1215 p = P._charpoly_coeff(r-1)
1216 # The _charpoly_coeff function already adds the factor of
1217 # -1 to ensure that _charpoly_coeff(r-1) is really what
1218 # appears in front of t^{r-1} in the charpoly. However,
1219 # we want the negative of THAT for the trace.
1220 return -p(*self.vector())
1221
1222
1223 def trace_inner_product(self, other):
1224 """
1225 Return the trace inner product of myself and ``other``.
1226 """
1227 if not other in self.parent():
1228 raise TypeError("'other' must live in the same algebra")
1229
1230 return (self*other).trace()
1231
1232
1233 class RealCartesianProductEJA(FiniteDimensionalEuclideanJordanAlgebra):
1234 """
1235 Return the Euclidean Jordan Algebra corresponding to the set
1236 `R^n` under the Hadamard product.
1237
1238 Note: this is nothing more than the Cartesian product of ``n``
1239 copies of the spin algebra. Once Cartesian product algebras
1240 are implemented, this can go.
1241
1242 EXAMPLES:
1243
1244 This multiplication table can be verified by hand::
1245
1246 sage: J = RealCartesianProductEJA(3)
1247 sage: e0,e1,e2 = J.gens()
1248 sage: e0*e0
1249 e0
1250 sage: e0*e1
1251 0
1252 sage: e0*e2
1253 0
1254 sage: e1*e1
1255 e1
1256 sage: e1*e2
1257 0
1258 sage: e2*e2
1259 e2
1260
1261 """
1262 @staticmethod
1263 def __classcall_private__(cls, n, field=QQ):
1264 # The FiniteDimensionalAlgebra constructor takes a list of
1265 # matrices, the ith representing right multiplication by the ith
1266 # basis element in the vector space. So if e_1 = (1,0,0), then
1267 # right (Hadamard) multiplication of x by e_1 picks out the first
1268 # component of x; and likewise for the ith basis element e_i.
1269 Qs = [ matrix(field, n, n, lambda k,j: 1*(k == j == i))
1270 for i in xrange(n) ]
1271
1272 fdeja = super(RealCartesianProductEJA, cls)
1273 return fdeja.__classcall_private__(cls, field, Qs, rank=n)
1274
1275 def inner_product(self, x, y):
1276 return _usual_ip(x,y)
1277
1278
1279 def random_eja():
1280 """
1281 Return a "random" finite-dimensional Euclidean Jordan Algebra.
1282
1283 ALGORITHM:
1284
1285 For now, we choose a random natural number ``n`` (greater than zero)
1286 and then give you back one of the following:
1287
1288 * The cartesian product of the rational numbers ``n`` times; this is
1289 ``QQ^n`` with the Hadamard product.
1290
1291 * The Jordan spin algebra on ``QQ^n``.
1292
1293 * The ``n``-by-``n`` rational symmetric matrices with the symmetric
1294 product.
1295
1296 * The ``n``-by-``n`` complex-rational Hermitian matrices embedded
1297 in the space of ``2n``-by-``2n`` real symmetric matrices.
1298
1299 * The ``n``-by-``n`` quaternion-rational Hermitian matrices embedded
1300 in the space of ``4n``-by-``4n`` real symmetric matrices.
1301
1302 Later this might be extended to return Cartesian products of the
1303 EJAs above.
1304
1305 TESTS::
1306
1307 sage: random_eja()
1308 Euclidean Jordan algebra of degree...
1309
1310 """
1311
1312 # The max_n component lets us choose different upper bounds on the
1313 # value "n" that gets passed to the constructor. This is needed
1314 # because e.g. R^{10} is reasonable to test, while the Hermitian
1315 # 10-by-10 quaternion matrices are not.
1316 (constructor, max_n) = choice([(RealCartesianProductEJA, 6),
1317 (JordanSpinEJA, 6),
1318 (RealSymmetricEJA, 5),
1319 (ComplexHermitianEJA, 4),
1320 (QuaternionHermitianEJA, 3)])
1321 n = ZZ.random_element(1, max_n)
1322 return constructor(n, field=QQ)
1323
1324
1325
1326 def _real_symmetric_basis(n, field=QQ):
1327 """
1328 Return a basis for the space of real symmetric n-by-n matrices.
1329 """
1330 # The basis of symmetric matrices, as matrices, in their R^(n-by-n)
1331 # coordinates.
1332 S = []
1333 for i in xrange(n):
1334 for j in xrange(i+1):
1335 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1336 if i == j:
1337 Sij = Eij
1338 else:
1339 # Beware, orthogonal but not normalized!
1340 Sij = Eij + Eij.transpose()
1341 S.append(Sij)
1342 return tuple(S)
1343
1344
1345 def _complex_hermitian_basis(n, field=QQ):
1346 """
1347 Returns a basis for the space of complex Hermitian n-by-n matrices.
1348
1349 TESTS::
1350
1351 sage: set_random_seed()
1352 sage: n = ZZ.random_element(1,5)
1353 sage: all( M.is_symmetric() for M in _complex_hermitian_basis(n) )
1354 True
1355
1356 """
1357 F = QuadraticField(-1, 'I')
1358 I = F.gen()
1359
1360 # This is like the symmetric case, but we need to be careful:
1361 #
1362 # * We want conjugate-symmetry, not just symmetry.
1363 # * The diagonal will (as a result) be real.
1364 #
1365 S = []
1366 for i in xrange(n):
1367 for j in xrange(i+1):
1368 Eij = matrix(field, n, lambda k,l: k==i and l==j)
1369 if i == j:
1370 Sij = _embed_complex_matrix(Eij)
1371 S.append(Sij)
1372 else:
1373 # Beware, orthogonal but not normalized! The second one
1374 # has a minus because it's conjugated.
1375 Sij_real = _embed_complex_matrix(Eij + Eij.transpose())
1376 S.append(Sij_real)
1377 Sij_imag = _embed_complex_matrix(I*Eij - I*Eij.transpose())
1378 S.append(Sij_imag)
1379 return tuple(S)
1380
1381
1382 def _quaternion_hermitian_basis(n, field=QQ):
1383 """
1384 Returns a basis for the space of quaternion Hermitian n-by-n matrices.
1385
1386 TESTS::
1387
1388 sage: set_random_seed()
1389 sage: n = ZZ.random_element(1,5)
1390 sage: all( M.is_symmetric() for M in _quaternion_hermitian_basis(n) )
1391 True
1392
1393 """
1394 Q = QuaternionAlgebra(QQ,-1,-1)
1395 I,J,K = Q.gens()
1396
1397 # This is like the symmetric case, but we need to be careful:
1398 #
1399 # * We want conjugate-symmetry, not just symmetry.
1400 # * The diagonal will (as a result) be real.
1401 #
1402 S = []
1403 for i in xrange(n):
1404 for j in xrange(i+1):
1405 Eij = matrix(Q, n, lambda k,l: k==i and l==j)
1406 if i == j:
1407 Sij = _embed_quaternion_matrix(Eij)
1408 S.append(Sij)
1409 else:
1410 # Beware, orthogonal but not normalized! The second,
1411 # third, and fourth ones have a minus because they're
1412 # conjugated.
1413 Sij_real = _embed_quaternion_matrix(Eij + Eij.transpose())
1414 S.append(Sij_real)
1415 Sij_I = _embed_quaternion_matrix(I*Eij - I*Eij.transpose())
1416 S.append(Sij_I)
1417 Sij_J = _embed_quaternion_matrix(J*Eij - J*Eij.transpose())
1418 S.append(Sij_J)
1419 Sij_K = _embed_quaternion_matrix(K*Eij - K*Eij.transpose())
1420 S.append(Sij_K)
1421 return tuple(S)
1422
1423
1424 def _mat2vec(m):
1425 return vector(m.base_ring(), m.list())
1426
1427 def _vec2mat(v):
1428 return matrix(v.base_ring(), sqrt(v.degree()), v.list())
1429
1430 def _multiplication_table_from_matrix_basis(basis):
1431 """
1432 At least three of the five simple Euclidean Jordan algebras have the
1433 symmetric multiplication (A,B) |-> (AB + BA)/2, where the
1434 multiplication on the right is matrix multiplication. Given a basis
1435 for the underlying matrix space, this function returns a
1436 multiplication table (obtained by looping through the basis
1437 elements) for an algebra of those matrices. A reordered copy
1438 of the basis is also returned to work around the fact that
1439 the ``span()`` in this function will change the order of the basis
1440 from what we think it is, to... something else.
1441 """
1442 # In S^2, for example, we nominally have four coordinates even
1443 # though the space is of dimension three only. The vector space V
1444 # is supposed to hold the entire long vector, and the subspace W
1445 # of V will be spanned by the vectors that arise from symmetric
1446 # matrices. Thus for S^2, dim(V) == 4 and dim(W) == 3.
1447 field = basis[0].base_ring()
1448 dimension = basis[0].nrows()
1449
1450 V = VectorSpace(field, dimension**2)
1451 W = V.span( _mat2vec(s) for s in basis )
1452
1453 # Taking the span above reorders our basis (thanks, jerk!) so we
1454 # need to put our "matrix basis" in the same order as the
1455 # (reordered) vector basis.
1456 S = tuple( _vec2mat(b) for b in W.basis() )
1457
1458 Qs = []
1459 for s in S:
1460 # Brute force the multiplication-by-s matrix by looping
1461 # through all elements of the basis and doing the computation
1462 # to find out what the corresponding row should be. BEWARE:
1463 # these multiplication tables won't be symmetric! It therefore
1464 # becomes REALLY IMPORTANT that the underlying algebra
1465 # constructor uses ROW vectors and not COLUMN vectors. That's
1466 # why we're computing rows here and not columns.
1467 Q_rows = []
1468 for t in S:
1469 this_row = _mat2vec((s*t + t*s)/2)
1470 Q_rows.append(W.coordinates(this_row))
1471 Q = matrix(field, W.dimension(), Q_rows)
1472 Qs.append(Q)
1473
1474 return (Qs, S)
1475
1476
1477 def _embed_complex_matrix(M):
1478 """
1479 Embed the n-by-n complex matrix ``M`` into the space of real
1480 matrices of size 2n-by-2n via the map the sends each entry `z = a +
1481 bi` to the block matrix ``[[a,b],[-b,a]]``.
1482
1483 EXAMPLES::
1484
1485 sage: F = QuadraticField(-1,'i')
1486 sage: x1 = F(4 - 2*i)
1487 sage: x2 = F(1 + 2*i)
1488 sage: x3 = F(-i)
1489 sage: x4 = F(6)
1490 sage: M = matrix(F,2,[[x1,x2],[x3,x4]])
1491 sage: _embed_complex_matrix(M)
1492 [ 4 -2| 1 2]
1493 [ 2 4|-2 1]
1494 [-----+-----]
1495 [ 0 -1| 6 0]
1496 [ 1 0| 0 6]
1497
1498 TESTS:
1499
1500 Embedding is a homomorphism (isomorphism, in fact)::
1501
1502 sage: set_random_seed()
1503 sage: n = ZZ.random_element(5)
1504 sage: F = QuadraticField(-1, 'i')
1505 sage: X = random_matrix(F, n)
1506 sage: Y = random_matrix(F, n)
1507 sage: actual = _embed_complex_matrix(X) * _embed_complex_matrix(Y)
1508 sage: expected = _embed_complex_matrix(X*Y)
1509 sage: actual == expected
1510 True
1511
1512 """
1513 n = M.nrows()
1514 if M.ncols() != n:
1515 raise ValueError("the matrix 'M' must be square")
1516 field = M.base_ring()
1517 blocks = []
1518 for z in M.list():
1519 a = z.real()
1520 b = z.imag()
1521 blocks.append(matrix(field, 2, [[a,b],[-b,a]]))
1522
1523 # We can drop the imaginaries here.
1524 return block_matrix(field.base_ring(), n, blocks)
1525
1526
1527 def _unembed_complex_matrix(M):
1528 """
1529 The inverse of _embed_complex_matrix().
1530
1531 EXAMPLES::
1532
1533 sage: A = matrix(QQ,[ [ 1, 2, 3, 4],
1534 ....: [-2, 1, -4, 3],
1535 ....: [ 9, 10, 11, 12],
1536 ....: [-10, 9, -12, 11] ])
1537 sage: _unembed_complex_matrix(A)
1538 [ 2*i + 1 4*i + 3]
1539 [ 10*i + 9 12*i + 11]
1540
1541 TESTS:
1542
1543 Unembedding is the inverse of embedding::
1544
1545 sage: set_random_seed()
1546 sage: F = QuadraticField(-1, 'i')
1547 sage: M = random_matrix(F, 3)
1548 sage: _unembed_complex_matrix(_embed_complex_matrix(M)) == M
1549 True
1550
1551 """
1552 n = ZZ(M.nrows())
1553 if M.ncols() != n:
1554 raise ValueError("the matrix 'M' must be square")
1555 if not n.mod(2).is_zero():
1556 raise ValueError("the matrix 'M' must be a complex embedding")
1557
1558 F = QuadraticField(-1, 'i')
1559 i = F.gen()
1560
1561 # Go top-left to bottom-right (reading order), converting every
1562 # 2-by-2 block we see to a single complex element.
1563 elements = []
1564 for k in xrange(n/2):
1565 for j in xrange(n/2):
1566 submat = M[2*k:2*k+2,2*j:2*j+2]
1567 if submat[0,0] != submat[1,1]:
1568 raise ValueError('bad on-diagonal submatrix')
1569 if submat[0,1] != -submat[1,0]:
1570 raise ValueError('bad off-diagonal submatrix')
1571 z = submat[0,0] + submat[0,1]*i
1572 elements.append(z)
1573
1574 return matrix(F, n/2, elements)
1575
1576
1577 def _embed_quaternion_matrix(M):
1578 """
1579 Embed the n-by-n quaternion matrix ``M`` into the space of real
1580 matrices of size 4n-by-4n by first sending each quaternion entry
1581 `z = a + bi + cj + dk` to the block-complex matrix
1582 ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into
1583 a real matrix.
1584
1585 EXAMPLES::
1586
1587 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1588 sage: i,j,k = Q.gens()
1589 sage: x = 1 + 2*i + 3*j + 4*k
1590 sage: M = matrix(Q, 1, [[x]])
1591 sage: _embed_quaternion_matrix(M)
1592 [ 1 2 3 4]
1593 [-2 1 -4 3]
1594 [-3 4 1 -2]
1595 [-4 -3 2 1]
1596
1597 Embedding is a homomorphism (isomorphism, in fact)::
1598
1599 sage: set_random_seed()
1600 sage: n = ZZ.random_element(5)
1601 sage: Q = QuaternionAlgebra(QQ,-1,-1)
1602 sage: X = random_matrix(Q, n)
1603 sage: Y = random_matrix(Q, n)
1604 sage: actual = _embed_quaternion_matrix(X)*_embed_quaternion_matrix(Y)
1605 sage: expected = _embed_quaternion_matrix(X*Y)
1606 sage: actual == expected
1607 True
1608
1609 """
1610 quaternions = M.base_ring()
1611 n = M.nrows()
1612 if M.ncols() != n:
1613 raise ValueError("the matrix 'M' must be square")
1614
1615 F = QuadraticField(-1, 'i')
1616 i = F.gen()
1617
1618 blocks = []
1619 for z in M.list():
1620 t = z.coefficient_tuple()
1621 a = t[0]
1622 b = t[1]
1623 c = t[2]
1624 d = t[3]
1625 cplx_matrix = matrix(F, 2, [[ a + b*i, c + d*i],
1626 [-c + d*i, a - b*i]])
1627 blocks.append(_embed_complex_matrix(cplx_matrix))
1628
1629 # We should have real entries by now, so use the realest field
1630 # we've got for the return value.
1631 return block_matrix(quaternions.base_ring(), n, blocks)
1632
1633
1634 def _unembed_quaternion_matrix(M):
1635 """
1636 The inverse of _embed_quaternion_matrix().
1637
1638 EXAMPLES::
1639
1640 sage: M = matrix(QQ, [[ 1, 2, 3, 4],
1641 ....: [-2, 1, -4, 3],
1642 ....: [-3, 4, 1, -2],
1643 ....: [-4, -3, 2, 1]])
1644 sage: _unembed_quaternion_matrix(M)
1645 [1 + 2*i + 3*j + 4*k]
1646
1647 TESTS:
1648
1649 Unembedding is the inverse of embedding::
1650
1651 sage: set_random_seed()
1652 sage: Q = QuaternionAlgebra(QQ, -1, -1)
1653 sage: M = random_matrix(Q, 3)
1654 sage: _unembed_quaternion_matrix(_embed_quaternion_matrix(M)) == M
1655 True
1656
1657 """
1658 n = ZZ(M.nrows())
1659 if M.ncols() != n:
1660 raise ValueError("the matrix 'M' must be square")
1661 if not n.mod(4).is_zero():
1662 raise ValueError("the matrix 'M' must be a complex embedding")
1663
1664 Q = QuaternionAlgebra(QQ,-1,-1)
1665 i,j,k = Q.gens()
1666
1667 # Go top-left to bottom-right (reading order), converting every
1668 # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1
1669 # quaternion block.
1670 elements = []
1671 for l in xrange(n/4):
1672 for m in xrange(n/4):
1673 submat = _unembed_complex_matrix(M[4*l:4*l+4,4*m:4*m+4])
1674 if submat[0,0] != submat[1,1].conjugate():
1675 raise ValueError('bad on-diagonal submatrix')
1676 if submat[0,1] != -submat[1,0].conjugate():
1677 raise ValueError('bad off-diagonal submatrix')
1678 z = submat[0,0].real() + submat[0,0].imag()*i
1679 z += submat[0,1].real()*j + submat[0,1].imag()*k
1680 elements.append(z)
1681
1682 return matrix(Q, n/4, elements)
1683
1684
1685 # The usual inner product on R^n.
1686 def _usual_ip(x,y):
1687 return x.vector().inner_product(y.vector())
1688
1689 # The inner product used for the real symmetric simple EJA.
1690 # We keep it as a separate function because e.g. the complex
1691 # algebra uses the same inner product, except divided by 2.
1692 def _matrix_ip(X,Y):
1693 X_mat = X.natural_representation()
1694 Y_mat = Y.natural_representation()
1695 return (X_mat*Y_mat).trace()
1696
1697
1698 class RealSymmetricEJA(FiniteDimensionalEuclideanJordanAlgebra):
1699 """
1700 The rank-n simple EJA consisting of real symmetric n-by-n
1701 matrices, the usual symmetric Jordan product, and the trace inner
1702 product. It has dimension `(n^2 + n)/2` over the reals.
1703
1704 EXAMPLES::
1705
1706 sage: J = RealSymmetricEJA(2)
1707 sage: e0, e1, e2 = J.gens()
1708 sage: e0*e0
1709 e0
1710 sage: e1*e1
1711 e0 + e2
1712 sage: e2*e2
1713 e2
1714
1715 TESTS:
1716
1717 The degree of this algebra is `(n^2 + n) / 2`::
1718
1719 sage: set_random_seed()
1720 sage: n = ZZ.random_element(1,5)
1721 sage: J = RealSymmetricEJA(n)
1722 sage: J.degree() == (n^2 + n)/2
1723 True
1724
1725 The Jordan multiplication is what we think it is::
1726
1727 sage: set_random_seed()
1728 sage: n = ZZ.random_element(1,5)
1729 sage: J = RealSymmetricEJA(n)
1730 sage: x = J.random_element()
1731 sage: y = J.random_element()
1732 sage: actual = (x*y).natural_representation()
1733 sage: X = x.natural_representation()
1734 sage: Y = y.natural_representation()
1735 sage: expected = (X*Y + Y*X)/2
1736 sage: actual == expected
1737 True
1738 sage: J(expected) == x*y
1739 True
1740
1741 """
1742 @staticmethod
1743 def __classcall_private__(cls, n, field=QQ):
1744 S = _real_symmetric_basis(n, field=field)
1745 (Qs, T) = _multiplication_table_from_matrix_basis(S)
1746
1747 fdeja = super(RealSymmetricEJA, cls)
1748 return fdeja.__classcall_private__(cls,
1749 field,
1750 Qs,
1751 rank=n,
1752 natural_basis=T)
1753
1754 def inner_product(self, x, y):
1755 return _matrix_ip(x,y)
1756
1757
1758 class ComplexHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1759 """
1760 The rank-n simple EJA consisting of complex Hermitian n-by-n
1761 matrices over the real numbers, the usual symmetric Jordan product,
1762 and the real-part-of-trace inner product. It has dimension `n^2` over
1763 the reals.
1764
1765 TESTS:
1766
1767 The degree of this algebra is `n^2`::
1768
1769 sage: set_random_seed()
1770 sage: n = ZZ.random_element(1,5)
1771 sage: J = ComplexHermitianEJA(n)
1772 sage: J.degree() == n^2
1773 True
1774
1775 The Jordan multiplication is what we think it is::
1776
1777 sage: set_random_seed()
1778 sage: n = ZZ.random_element(1,5)
1779 sage: J = ComplexHermitianEJA(n)
1780 sage: x = J.random_element()
1781 sage: y = J.random_element()
1782 sage: actual = (x*y).natural_representation()
1783 sage: X = x.natural_representation()
1784 sage: Y = y.natural_representation()
1785 sage: expected = (X*Y + Y*X)/2
1786 sage: actual == expected
1787 True
1788 sage: J(expected) == x*y
1789 True
1790
1791 """
1792 @staticmethod
1793 def __classcall_private__(cls, n, field=QQ):
1794 S = _complex_hermitian_basis(n)
1795 (Qs, T) = _multiplication_table_from_matrix_basis(S)
1796
1797 fdeja = super(ComplexHermitianEJA, cls)
1798 return fdeja.__classcall_private__(cls,
1799 field,
1800 Qs,
1801 rank=n,
1802 natural_basis=T)
1803
1804 def inner_product(self, x, y):
1805 # Since a+bi on the diagonal is represented as
1806 #
1807 # a + bi = [ a b ]
1808 # [ -b a ],
1809 #
1810 # we'll double-count the "a" entries if we take the trace of
1811 # the embedding.
1812 return _matrix_ip(x,y)/2
1813
1814
1815 class QuaternionHermitianEJA(FiniteDimensionalEuclideanJordanAlgebra):
1816 """
1817 The rank-n simple EJA consisting of self-adjoint n-by-n quaternion
1818 matrices, the usual symmetric Jordan product, and the
1819 real-part-of-trace inner product. It has dimension `2n^2 - n` over
1820 the reals.
1821
1822 TESTS:
1823
1824 The degree of this algebra is `n^2`::
1825
1826 sage: set_random_seed()
1827 sage: n = ZZ.random_element(1,5)
1828 sage: J = QuaternionHermitianEJA(n)
1829 sage: J.degree() == 2*(n^2) - n
1830 True
1831
1832 The Jordan multiplication is what we think it is::
1833
1834 sage: set_random_seed()
1835 sage: n = ZZ.random_element(1,5)
1836 sage: J = QuaternionHermitianEJA(n)
1837 sage: x = J.random_element()
1838 sage: y = J.random_element()
1839 sage: actual = (x*y).natural_representation()
1840 sage: X = x.natural_representation()
1841 sage: Y = y.natural_representation()
1842 sage: expected = (X*Y + Y*X)/2
1843 sage: actual == expected
1844 True
1845 sage: J(expected) == x*y
1846 True
1847
1848 """
1849 @staticmethod
1850 def __classcall_private__(cls, n, field=QQ):
1851 S = _quaternion_hermitian_basis(n)
1852 (Qs, T) = _multiplication_table_from_matrix_basis(S)
1853
1854 fdeja = super(QuaternionHermitianEJA, cls)
1855 return fdeja.__classcall_private__(cls,
1856 field,
1857 Qs,
1858 rank=n,
1859 natural_basis=T)
1860
1861 def inner_product(self, x, y):
1862 # Since a+bi+cj+dk on the diagonal is represented as
1863 #
1864 # a + bi +cj + dk = [ a b c d]
1865 # [ -b a -d c]
1866 # [ -c d a -b]
1867 # [ -d -c b a],
1868 #
1869 # we'll quadruple-count the "a" entries if we take the trace of
1870 # the embedding.
1871 return _matrix_ip(x,y)/4
1872
1873
1874 class JordanSpinEJA(FiniteDimensionalEuclideanJordanAlgebra):
1875 """
1876 The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)``
1877 with the usual inner product and jordan product ``x*y =
1878 (<x_bar,y_bar>, x0*y_bar + y0*x_bar)``. It has dimension `n` over
1879 the reals.
1880
1881 EXAMPLES:
1882
1883 This multiplication table can be verified by hand::
1884
1885 sage: J = JordanSpinEJA(4)
1886 sage: e0,e1,e2,e3 = J.gens()
1887 sage: e0*e0
1888 e0
1889 sage: e0*e1
1890 e1
1891 sage: e0*e2
1892 e2
1893 sage: e0*e3
1894 e3
1895 sage: e1*e2
1896 0
1897 sage: e1*e3
1898 0
1899 sage: e2*e3
1900 0
1901
1902 """
1903 @staticmethod
1904 def __classcall_private__(cls, n, field=QQ):
1905 Qs = []
1906 id_matrix = identity_matrix(field, n)
1907 for i in xrange(n):
1908 ei = id_matrix.column(i)
1909 Qi = zero_matrix(field, n)
1910 Qi.set_row(0, ei)
1911 Qi.set_column(0, ei)
1912 Qi += diagonal_matrix(n, [ei[0]]*n)
1913 # The addition of the diagonal matrix adds an extra ei[0] in the
1914 # upper-left corner of the matrix.
1915 Qi[0,0] = Qi[0,0] * ~field(2)
1916 Qs.append(Qi)
1917
1918 # The rank of the spin algebra is two, unless we're in a
1919 # one-dimensional ambient space (because the rank is bounded by
1920 # the ambient dimension).
1921 fdeja = super(JordanSpinEJA, cls)
1922 return fdeja.__classcall_private__(cls, field, Qs, rank=min(n,2))
1923
1924 def inner_product(self, x, y):
1925 return _usual_ip(x,y)
My personal library of Sage code.